Unexpected Non-Uniform Flux with Homogeneous Neumann BC

General
1

Hi, I’m learning fenicsx-0.9.0 to work on some 3D problems. I found that when I apply a homogeneous Neumann boundary condition to the 3D domain, the actual Neumann flux ends up pretty uneven, depending on the mesh cell type. Here’s a minimal example to show what I mean:

from pathlib import Path
from mpi4py import MPI
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem

import ufl
import numpy as np

msh = mesh.create_box(
    comm = MPI.COMM_WORLD,
    points = ((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)),
    n = (10, 10, 10),
)

V = fem.functionspace(msh, ('CG', 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(msh,0.1)

# Define ds(1)
boundary = lambda x : np.isclose(x[2],10.0)
facets = mesh.locate_entities_boundary(msh, 2, boundary)
marks = np.ones((len(facets),),dtype=np.int32)
meshtag = mesh.meshtags(msh,2,facets,marks)
ds = ufl.Measure("ds", domain = msh, subdomain_data = meshtag)

# Define and solve the problem
a = u * v * ufl.dx
L = f * v * ds(1)
problem = LinearProblem(
    a,L,petsc_options={"ksp_type":"gmres","pc_type": "ilu"})
uh = problem.solve()

# Output xdmf results
out_folder = Path("output")
out_folder.mkdir(parents=True, exist_ok=True)
with io.XDMFFile(msh.comm, out_folder / "results.xdmf", "w") as file:
    file.write_mesh(msh)
    file.write_function(uh)

This minumum working code simply solves:

uh*v*dx = f*v*ds(1) with spatially constant f

and in thoery, the value of uh on the boundary ds(1) should be uniform. However, when I use the built-in structured tetrahedral mesh, the values for uh on the boundary ds(1) turn out uneven (see the first figure). If I switch to an unstructured tetra mesh, the results get worse (see the second figure). On the contrary, if I set the cell type as hexahedron, the results would be perfect (see the thrid figure).

structured tetra
structured tetra1807×1021 277 KB

unstructured tetra
unstructured tetra1812×1033 300 KB

quad
quad1809×1024 213 KB

So how can I properly apply the Neumann boundary condition when using a tetrahedral mesh? Did I miss something in my setup?

2

It is not quite so simple. FEM operates in weak forms; equations are satisfied only in terms of their actions on testfunctions (often states as “in a weighted average sense”).

Your simple-looking problem is balancing u^h = 0 inside of the domain with u^h = f on domain boundary ds(1). Both cannot hold at the same time (specifically on the element adjacent ds(1)), and hence the solution that satisfies this in a “weighted average sense” is one that satisfies neither exactly and wiggles a bit.

(BTW. this has little to do with Neumann BC’s at this point)

3

Hi, Stein, thanks for your reply. Sorry for my oversimplification of my problem. I used to think that
uh*v*dx = f*v*ds(1)
is equivalent to, for exampel, the weak form of Richards equation
C*(uh-uh0)*v*dx = -K*dot(grad(uh+z),grad(v))*dx + f*v*ds(1) with C=Constant(1), uh0=Constant(0), K=Constant(0).
I have corrected the weak form in my code to a more meaningful one:

from pathlib import Path
from mpi4py import MPI
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
import ufl
import numpy as np

msh = mesh.create_box(
    comm = MPI.COMM_WORLD,
    points = ((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)),
    n = (10, 10, 10),
)

V = fem.functionspace(msh, ('CG', 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(msh,0.1)

# Define ds(1)
boundary = lambda x : np.isclose(x[2],10.0)
facets = mesh.locate_entities_boundary(msh, 2, boundary)
marks = np.ones((len(facets),),dtype=np.int32)
meshtag = mesh.meshtags(msh,2,facets,marks)
ds = ufl.Measure("ds", domain = msh, subdomain_data = meshtag)

# Define and solve the problem
u0 = fem.Constant(msh, 0.0)
K = fem.Constant(msh, 1.0)
a = u * v * ufl.dx + K * ufl.dot(ufl.grad(u),ufl.grad(v)) * ufl.dx
L = u0 * v * ufl.dx + f * v * ds(1)

problem = LinearProblem(
a,L,petsc_options={"ksp_type":"gmres","pc_type": "ilu"})
uh = problem.solve()

# Output xdmf results
out_folder = Path("output")
out_folder.mkdir(parents=True, exist_ok=True)
with io.XDMFFile(msh.comm, out_folder / "results.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(uh)

Basing on the new example, I found that the “unexpected non-uniform flux” is simply a issue of mesh resoultion. If I modify

msh = mesh.create_box(
    comm = MPI.COMM_WORLD,
    points = ((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)),
    n = (10, 10, 10))

to

msh = mesh.create_box(
    comm = MPI.COMM_WORLD,
    points = ((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)),
    n = (10, 10, 100))

the outcome would be much better.

4

Yes indeed, which essentially comes down to the same point;

In the case of PDEs and natural/essential boundary conditions, you’re also trying to satisfy both the PDE in the interior and the BC on the boundary.

Of course, the more you mesh refine, the more the discrete solution can satisfy both at the same time and the smaller the missmatches become.