Neumann BC not fulfilled? (dolfinx)

I had some indications that in one of my dolfinx-scripts the Neumann BC is not fulfilled, so I barely modified one of the tutorials from @dokken to create a MWE and was surprised to find the same indication there as well: The Neumann BC does not seem to be fulfilled.

MWE (see Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial):

import numpy as np
import pyvista

from dolfinx.fem import (Constant,  Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
u_ex = lambda x: 1 + x[0]**2 + 2*x[1]**2
x = SpatialCoordinate(mesh)
# Define physical parameters and boundary condtions
s = u_ex(x)
f = -div(grad(u_ex(x)))
n = FacetNormal(mesh)
g = -dot(n, grad(u_ex(x)))
# g = 0.0

kappa = Constant(mesh, ScalarType(1))
r = Constant(mesh, ScalarType(1000))
# Define function space and standard part of variational form
V = FunctionSpace(mesh, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
F = kappa * inner(grad(u), grad(v)) * dx - inner(f, v) * dx

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(mesh.comm, "output/facet_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tag)

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
                       BoundaryCondition("Dirichlet", 2, u_ex),
                       BoundaryCondition("Robin", 3, (r, s)),
                       BoundaryCondition("Neumann", 4, g)]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem 
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

solution_trace_norm = assemble_scalar(form((dot(n, grad(uh))-g)**2*ds(4)))**0.5
print(f"Neumann BC trace: {solution_trace_norm}")

The output I receive should be close to zero, but is:

Neumann BC trace: 7.800058502092114

Even if I set the Neumann BC to g=0, the output is still too high with

Neumann BC trace: 1.366936098344996

I am using the version fenics-dolfinx 0.6.0 via conda-forge on Ubuntu via WSL2.

Is this expected behaviour?
Did I do any mistakes evaluating (last line of the code) the Neumann BC?

1 Like

There is a sign mistake, your g is pointing inwards, not outwards.
As you can see from the following example, the trace goes closer and closer to the actual trace.

import numpy as np

from dolfinx.fem import (Constant,  Function, functionspace, assemble_scalar,
                         dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile

for N in [10, 20, 40, 80, 100]:
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    def u_ex(x): return 1 + x[0]**2 + 2*x[1]**2

    x = SpatialCoordinate(mesh)
    # Define physical parameters and boundary condtions
    s = u_ex(x)
    f = -div(grad(u_ex(x)))
    n = FacetNormal(mesh)
    g = -dot(n, grad(u_ex(x)))
    # g = 0.0

    kappa = Constant(mesh, ScalarType(1))
    r = Constant(mesh, ScalarType(1000))
    # Define function space and standard part of variational form
    V = functionspace(mesh, ("Lagrange", 1))
    u, v = TrialFunction(V), TestFunction(V)
    F = kappa * inner(grad(u), grad(v)) * dx - inner(f, v) * dx

    boundaries = [(1, lambda x: np.isclose(x[0], 0)),
                  (2, lambda x: np.isclose(x[0], 1)),
                  (3, lambda x: np.isclose(x[1], 0)),
                  (4, lambda x: np.isclose(x[1], 1))]

    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    for (marker, locator) in boundaries:
        facets = locate_entities(mesh, fdim, locator)
        facet_indices.append(facets)
        facet_markers.append(np.full_like(facets, marker))
    facet_indices = np.hstack(facet_indices).astype(np.int32)
    facet_markers = np.hstack(facet_markers).astype(np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = meshtags(
        mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

    mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

    ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

    class BoundaryCondition():
        def __init__(self, type, marker, values):
            self._type = type
            if type == "Dirichlet":
                u_D = Function(V)
                u_D.interpolate(values)
                facets = facet_tag.find(marker)
                dofs = locate_dofs_topological(V, fdim, facets)
                self._bc = dirichletbc(u_D, dofs)
            elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
            elif type == "Robin":
                self._bc = values[0] * inner(u-values[1], v) * ds(marker)
            else:
                raise TypeError(
                    "Unknown boundary condition: {0:s}".format(type))

        @property
        def bc(self):
            return self._bc

        @property
        def type(self):
            return self._type

    # Define the Dirichlet condition
    boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
                           BoundaryCondition("Dirichlet", 2, u_ex),
                           BoundaryCondition("Robin", 3, (r, s)),
                           BoundaryCondition("Neumann", 4, g)]

    bcs = []
    for condition in boundary_conditions:
        if condition.type == "Dirichlet":
            bcs.append(condition.bc)
        else:
            F += condition.bc

    # Solve linear variational problem
    a = lhs(F)
    L = rhs(F)
    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
                            )
    uh = problem.solve()
    print(
        f"Integral over boundary {assemble_scalar(form(dot(n, grad(uh))*ds(4)))},{assemble_scalar(form(-g*ds(4)))}")
    solution_trace_norm = assemble_scalar(
        form((dot(n, grad(uh))+g)**2*ds(4)))**0.5
    print(f"Neumann BC trace: {solution_trace_norm}")

    print(f"L2 norm: {assemble_scalar(form(dot(uh-u_ex(x), uh-u_ex(x))*dx))}")
    with XDMFFile(mesh.comm, "output/facet_tags.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_meshtags(facet_tag, mesh.geometry)
        xdmf.write_function(uh)

(here Im using the main branch of DOLFINx), so a slight change in syntax:

Integral over boundary 3.8000585020380435,4.0
Neumann BC trace: 0.19994150007125447
L2 norm: 2.3597310792740984e-05
Integral over boundary 3.900007213125378,4.000000000000002
Neumann BC trace: 0.09999278693537393
L2 norm: 1.4778471852531955e-06
Integral over boundary 3.9500008985183026,4.000000000000001
Neumann BC trace: 0.049999101483557185
L2 norm: 9.242779112544334e-08
Integral over boundary 3.9750001122171454,3.999999999999993
Neumann BC trace: 0.024999887782911453
L2 norm: 5.777939472233128e-09
Integral over boundary 3.9800000574491663,4.000000000000002
Neumann BC trace: 0.019999942550853547
L2 norm: 2.366711199034677e-09