Discrepancy, with assemble and assmbled operator

[tucker]

I have created a mesh in gmsh for the Haut Glacier d’Arolla (TC - Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP–HOM)) and converted to .xml format for use in FEniCS using dolfin-convert.

I decided to test my mesh by solving the following Poisson problem.

After solving the discretized problem, checking that the solution is reasonable by visualizing in ParaView I then check the H1 semi-norm of the solution, defined as

in two ways. First by using the assemble command in dolfin, the numerical value returned is nan. There is an alternative way to compute the H1 semi-norm. That is

Where A is the finite element stiffness matrix. u^TA u evaluates to a reasonable nonnegative finite number.

I have included the full code below., but of course you will not be able to test without the mesh.

I am hoping to get insight into why this would be the case, what goes on internally in FEniCS when one assembles a finite element system matrix, such as A, with Dirichlet conditions that causes the discrepancy that I am seeing.

import dolfin as dl
import numpy as np

# mesh of PDE domain
meshdir = "../"
meshname = "Arolla2"
mesh = dl.Mesh(meshdir+meshname+".xml")
boundary_markers = dl.MeshFunction("size_t", mesh,\
                   meshdir+meshname+"_facet_region.xml")

# discrete function spaces
P1 = dl.FiniteElement("CG", mesh.ufl_cell(), 1)
Vh  = dl.FunctionSpace(mesh, P1)

# Homogeneous Dirichlet conditions
bc = [dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 1),\
      dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 2)]

# Right hand side forcing term
f = dl.Constant(1.0)

vh = dl.TestFunction(Vh)
uh = dl.TrialFunction(Vh)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx
A, b = dl.assemble_system(a, L, bc)

Asolver = dl.PETScKrylovSolver()
dl.PETScOptions.set("ksp_monitor")
dl.PETScOptions.set("ksp_monitor_true_residual")
dl.PETScOptions.set("ksp_type", "cg")
dl.PETScOptions.set("ksp_rtol", "1.e-12")
dl.PETScOptions.set("ksp_atol", "1.e-12")
Asolver.set_from_options()
Asolver.set_operator(A)

u = dl.Function(Vh)
Asolver.solve(u.vector(), b)

# for the determined solution,
# compute sqrt(int_(Omega) ||grad(u)||^2 dV) 
u_H1norm_1 = np.sqrt(dl.assemble(dl.inner(dl.grad(u), dl.grad(u))*dl.dx(mesh))) 
print("||u||_H^1(Omega)^2 = %1.3e (directly assembled)" %(u_H1norm_1))

# vector to hold the product of the solution with
# the finite element system matrix A
Au = dl.Vector()
A.init_vector(Au, 0)

# Au = A*u
A.mult(u.vector(), Au)

# sqrt(u^T Au)
u_H1norm_2 = np.sqrt(Au.inner(u.vector()))
print("||u||_H^1(Omega) = %1.3e (from assembled operators)" \
             %(u_H1norm_2))

[dokken]

I cannot reproduce this with a built in mesh:

import dolfin as dl
import numpy as np

mesh = dl.UnitSquareMesh(30,30)
# discrete function spaces
P1 = dl.FiniteElement("CG", mesh.ufl_cell(), 1)
Vh  = dl.FunctionSpace(mesh, P1)

# Homogeneous Dirichlet conditions
bc = [dl.DirichletBC(Vh, dl.Constant(0.0), "on_boundary")]

# Right hand side forcing term
f = dl.Constant(1.0)

vh = dl.TestFunction(Vh)
uh = dl.TrialFunction(Vh)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx
A, b = dl.assemble_system(a, L, bc)

Asolver = dl.PETScKrylovSolver()
dl.PETScOptions.set("ksp_monitor")
dl.PETScOptions.set("ksp_monitor_true_residual")
dl.PETScOptions.set("ksp_type", "cg")
dl.PETScOptions.set("ksp_rtol", "1.e-12")
dl.PETScOptions.set("ksp_atol", "1.e-12")
Asolver.set_from_options()
Asolver.set_operator(A)

u = dl.Function(Vh)
Asolver.solve(u.vector(), b)

# for the determined solution,
# compute sqrt(int_(Omega) ||grad(u)||^2 dV) 
u_H1norm_1 = np.sqrt(dl.assemble(dl.inner(dl.grad(u), dl.grad(u))*dl.dx(mesh))) 
print("||u||_H^1(Omega)^2 = %1.3e (directly assembled)" %(u_H1norm_1))

# vector to hold the product of the solution with
# the finite element system matrix A
Au = dl.Vector()
A.init_vector(Au, 0)

# Au = A*u
A.mult(u.vector(), Au)

# sqrt(u^T Au)
u_H1norm_2 = np.sqrt(Au.inner(u.vector()))
print("||u||_H^1(Omega) = %1.3e (from assembled operators)" \
             %(u_H1norm_2))

which returns:

||u||_H^1(Omega)^2 = 1.871e-01 (directly assembled)
||u||_H^1(Omega) = 1.871e-01 (from assembled operators)

(using docker and: (docker run -it -v $(pwd):/home/fenics/shared -w /home/fenics/shared --rm quay.io/fenicsproject/dev:latest)
Therefore, it is most likely your mesh that has issues.
I also reproduced it with an unstructured mesh, using mshr:

cylinder = mshr.Cylinder(dl.Point(0, 0, 0), dl.Point(0, 0, -9), 8.0, 8.0)
domain = cylinder
mesh = mshr.generate_mesh(domain, 20)

with result:

||u||_H^1(Omega)^2 = 7.148e+01 (directly assembled)
||u||_H^1(Omega) = 7.148e+01 (from assembled operators)

[tucker]

The issue appears to have been due to an issue with how I generated the mesh in gmsh. I was still hoping to get insight into the specifics of how a matrix is altered when Dirichlet conditions are applied, as I expected the results above to be the same.

[dokken]

To assemble a matrix system symmetrically with Dirichlet conditions, the most common thing to do is to apply lifting of the boundary conditions. This is done by ignoring all contributions from cells containing Dirichlet-dofs during assembly, and subtract the matrix contribution for the Dirichlet dofs from the RHS.
Each row or column corresponding to a Dirichlet condition will simply have one on the diagonal, and be zero everywhere else.

[tucker]

Thank you so much, this is helpful.