Compute H^2 norm

10 months ago by
Hi everyone,

I would like to compute the H^2 norm of a field. Let say $w$w   is the numerical solution of a Poisson problem with homogenous Dirichlet boundary condition, a constant source term  $f$ƒ   and a piecewise constant coefficient  $\mu$μ .

  $\mu\Delta w=f$μΔw=ƒ  

This solution is computed with Lagrange finite element of order 2.

N = 128
mesh = UnitSquareMesh(N,N)
V2 = FunctionSpace(mesh,'CG',2)
def no_slip_bound(x,on_boundary):
    return on_boundary
bc = DirichletBC(V2,Constant(0.),no_slip_bound)
# Parameters
f = Constant(1.)
mu1 = 1.
mu2 = 1.
mu = Expression("mu1 + mu2*(x[0]>=0.5)",mu1=mu1,mu2=mu2,degree=0)
# Weak formulation
u = TrialFunction(V2)
v = TestFunction(V2)
w = Function(V2)
a = mu*inner(grad(u),grad(v))*dx
L = dot(f,v)*dx

I found two possible ways to compute its H^2 norm, but they do not give the same result :

First solution :
V1 = FunctionSpace(mesh,"Lagrange",1)
V0 = FunctionSpace(mesh,"DG",0)

# Compute derivatives
dxw = w.dx(0); dxw = project(dxw,V1)
dyw = w.dx(1); dyw = project(dyw,V1)
dxxw = dxw.dx(0); dxxw = project(dxxw,V0)
dxyw = dxw.dx(1); dxyw = project(dxyw,V0)
dyxw = dyw.dx(0); dyxw = project(dyxw,V0)
dyyw = dyw.dx(1); dyyw = project(dyyw,V0)
# Compute L2, H1 and H2 norms
normL2 = sqrt(assemble(dot(w,w)*dx))
normH1 = sqrt(assemble(dot(w,w)*dx) + assemble(dot(dxw,dxw)*dx + dot(dyw,dyw)*dx))
normH2 = sqrt(assemble((dot(w,w) + dot(dxw,dxw) + dot(dyw,dyw) + dot(dxxw,dxxw) + dot(dxyw,dxyw) + dot(dyxw,dyxw) + dot(dyyw,dyyw))*dx))

Second solution :
normL2 = sqrt(assemble(dot(w,w)*dx))
normH1 = sqrt(assemble(dot(w,w)*dx) + assemble((dot(w.dx(0),w.dx(0)) + dot(w.dx(1),w.dx(1)))*dx))
normH2 = sqrt(assemble((dot(w,w) + dot(w.dx(0),w.dx(0)) + dot(w.dx(1),w.dx(1)) + dot(w.dx(0).dx(0),w.dx(0).dx(0)) + dot(w.dx(0).dx(1),w.dx(0).dx(1)) + dot(w.dx(1).dx(0),w.dx(1).dx(0)) + dot(w.dx(1).dx(1),w.dx(1).dx(1)))*dx))​

In the first case, the H^2 norm is not bounded, which is what we should expect, since the coefficient in the Poisson equation is not continuous.

However, in the second case, the H^2 norm is bounded.

Questions :
Why the two computed H^2 norms are not the same?
2. Can we really compute the H^2 norm with Lagrange elements, since they are not C^1 ? ( I am a bit confused about that, sorry...)

Thank you for your help.

Community: FEniCS Project

1 Answer

10 months ago by
It happens because you are projecting the functions to a smoother space, which in result gives a different function. If your field u is in \( H^2 \), then of course you would have \( D u \in H^1 \) and \( D^2 u \in L^2 \), but instead your are projecting the second derivatives in a second degree space. Please try using 

V0 = FonctionSpace(mesh,"DG",0)​

instead of the one you are using, surely the norms will look similar. There is also the chance that you are having numerical artifacts, so maybe compare this methods to calculating the norm with

normH2 = norm(u, 'H1') + norm(grad(u), 'H10', mesh = mesh)​

where the 0 indicates that you are not getting the \( L^2 \) term ( and the second one is receiving the mesh because grad(u) is no longer a Function object.... but maybe it won't run anyway. If it doesn't, try projecting first the gradient as you do with the partial derivatives to \( H^1 \). Let me know how it goes. 

Best regards

Thank you for your response Nicolas.

Concerning your first remark, I made a mistake in writing my question... (sorry about that ! I edited my question). In my code I have

V0 = FunctionSpace(mesh,"DG",0)​​

Concerning your second remark, calculting the norm with

normH2 = norm(w, 'H1') + norm(grad(w), 'H10', mesh = mesh)

do not run. However, projecting first the gradient with the partial derivatives to  $H^1$H1 does,

normH2 = norm(w, 'H1') + norm(dxw, 'H10') + norm(dyw, 'H10')

The bad news is that it gives a different results from the other two ways of computing the  $H^2$H2  norm and I do not know which one is correct.


written 10 months ago by Fabien Vergnet  
Haha it is good news! I guess they are different because projections change the solution a bit (not too much). Maybe you could post the numeric values you are getting to see how much it changes, but all of them should be usable. I tried a simple case myself:
from fenics import *

N = 100
mesh = UnitSquareMesh(N,N)
f = Expression("x[0] + x[1]", degree = 1, domain = mesh)
dx = Measure('dx', domain = mesh)
Df = grad(f)
DDf = grad(grad(f))
norm1 = sqrt(assemble( (f*f + dot(Df, Df) + inner(DDf, DDf))*dx ))

V1 = FunctionSpace(mesh, 'CG', 1)
V0 = FunctionSpace(mesh, 'DG', 0)
dxf = f.dx(0); dxf = project(dxf,V1)
dyf = f.dx(1); dyf = project(dyf,V1)
dxxf = dxf.dx(0); dxxf = project(dxxf,V0)
dxyf = dxf.dx(1); dxyf = project(dxyf,V0)
dyxf = dyf.dx(0); dyxf = project(dyxf,V0)
dyyf = dyf.dx(1); dyyf = project(dyyf,V0)

norm2 = sqrt(assemble( (f*f +  dxf*dxf + dyf*dyf + dxxf*dxxf + \
            dxyf*dxyf + dyxf*dyxf + dyyf*dyyf)*dx))
gradf = as_vector([dxf, dyf])
gradf = project(gradf, VectorFunctionSpace(mesh, 'CG', 1))
norm3 = norm(f, 'H1') + norm(gradf, 'H10')

print norm1, norm2, norm3, sqrt(2./3 + 1./2 + 2.)​

where the last term is the exact solution, and I get the same value for all three. In any case, I recommend using the first one I did here, of maybe doing 

norm(f, 'H1') + sqrt(assemble(inner(DDf, DDf)*dx))

because it is cleaner. Also, it avoids projecting, which will change a little bit your numbers (but shouldn't be wrong!), as fenics will (behiind curtains) use the exact values per element, because it just takes derivatives of the local polynomial expression. Try creating functions and comparing the projections vs interpolations to see this effect. 

Best regards!

written 10 months ago by Nicolás Barnafi  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »