Dear NGSolve community,
recently, I got reason to numerically approximate high order derivatives (like the bilaplacian) of a phase field function.
To my surprise, I got into some trouble doing so and hope for being pointed out
a major fault in my method of using NGSolve or directed to a better method. The code I have been using to debug the issue is attached.
I noticed that the laplacian - and accordingly the bilaplacian - approximates show oscillatory behaviour. Sure, the phase field parameter I employ may be small and so I played around with making it larger or, on the other hand, refining the mesh while keeping it small.
Both approaches did not lead to a satisfactory result: The oscillatory behaviour still remains outside the
interfacial layer. To my understanding, computing gradients and divergences with NGSolve like in the attached code is based on differences of the phase field on a reference element. Am I mistaken here?
I am not very experienced in numerical issues so any recommendations will be appreciated.
All the best,
Edit: I just substituted the phase field with the constant function 1 and even then the laplacian is not approximated correctly (Image attached), so there must be some error in using NGSolve I suppose.