Second-derivative incorrect on boundary

Question

Hi all,

I thought this should be a simple problem, but I am having some trouble so hoping someone can help. I have a parabolic function defined on a unit square (like a fully-developed steady 2D fluid flow). I want to calculate the second derivative of this function. Since ${\rm vel} = 4y - 4y^2$, I know the answer should be -8 everywhere. My code is given below. What I get is +8 everywhere, except at the (top/bottom) boundaries where vel is zero, where I get large noise-like results. The thickness of the noisy "boundary layer" is wider for a lower mesh density and gets worse if I increase the element order. I tried including the boundary term which arises after integration by parts, but this introduces error at the left/right boundaries. If I make vel a function of FS rather than VFS, then the ()*ds term does help, but I don't understand why.

Any help would be greatly appreciated.

mesh = UnitSquareMesh(100, 100)
VFS = VectorFunctionSpace(mesh, "CG", 1)
FS = FunctionSpace(mesh, "CG", 1)
velexp = Expression(('4.0*x[1]*(1 - x[1])', '0'), element=VFS.ufl_element())
vel = Function(VFS)
vel.interpolate(velexp)
du2 = Function(VFS)
u = TestFunction(VFS)
v = TrialFunction(VFS)

a0 = assemble(inner(v, u)*dx)
L0 = assemble(inner(grad(v), grad(vel))*dx )

# L0 = assemble(inner(grad(v), -grad(vel))*dx )  # Adding negative sign does nothing
# n  = FacetNormal(mesh)
# L0 = assemble(inner(grad(v), grad(vel))*dx - inner(v, grad(vel).T*n)*ds )

solve(a0, du2.vector(), L0,  'gmres', 'ilu')
File("du2.pvd") << du2

Comments

Hi, you have a better chance for an answer if you use math (like $\LaTeX$) and code formatting (the {} in the editor toolbar) for your question.

---

Thanks, I've put the code in {}, but I don't know how to use math. I tried $ ... $ but that didn't work.

---

It should though. add $ $ to x^2 to get $x^2$. Maybe it is just the preview that doesn't work.

---

Yep, it was just the preview that was misleading. Thanks!

Answer 1

This is correct term

L0 = assemble(inner(grad(v), grad(vel))*dx - inner(v, grad(vel)*n)*ds )

resulting from integration by parts.

Then the values in the bulk are correct. Regarding the oscillations on the boundary, this is tricky thing I haven't satisfactorily solved also. The nature of the problem is that vel is a function on the space which is not $H^2$-conforming and you want to compute a second derivative.

You could simply get rid of boundary values. For another (incomplete) suggestions check this thread on scicomp.

Answer 2

The problem with the sign being wrong was just an embarrassing mistake (I was only viewing the magnitude plot in Paraview!).

Regarding the boundary problems, I noticed that with cubic or higher degree polynomials, i.e. ${\rm vel} = 4y^k - y^{k+1}$ (where k > 1), the solution is correct at the $y=0$ boundary if I do not include the boundary flux term... however, not including this term causes the solution at the $y=1$ boundary to be worse.