Assembly of integrals over interal boundaries messes up DirichletBCs

Joseph Weston

I believe that I have encountered a similar problem.

The MWE that I have is:

Solve the Poisson equation on a unit square, with homogeneous Dirichlet BCs applied on the top/bottom of the square.

Let the RHS of the Poisson equation be given by an integral over internal facets, on a line bisecting the square parallel to the x axis.

Note that this case is slightly different to the OP’s: I apply Dirichlet BCs on external facets, but the RHS consists of an integral over internal facets.

Nevertheless, I believe the underlying cause to be the same, as different results are obtained when using SystemAssembler (called byassemble_system) or Assembler (called by assemble.

I observed this issue on version 2019.1.0, but have not tested against other versions.

What follows is a script that implements the MWE (it seems Bitbucket does not allow commenters to attach files in a way that do not modify the OP):

# Electrostatics surface charge capacitor test

# We are going to solve the Poisson equation on a unit square domain $\Omega$, where the charge density is zero everywhere except on a line $\Gamma$ that bisects the square at $y=0.7$.
# $$
# \int_\Omega \epsilon_r(x)\nabla u \cdot \nabla v \; dx = \int_\Gamma \rho \; dS
# $$


from matplotlib import pyplot as plt
import numpy as np
import fenics


class YSubdomain(fenics.SubDomain):

    def __init__(self, y):
        self.y = float(y)
        super().__init__()


class Interface(YSubdomain):
    def inside(self, x, on_boundary):
        return fenics.near(x[1], self.y, 1e-6)


class Region(YSubdomain):
    def inside(self, x, on_boundary):
        return x[1] <= self.y


mesh = fenics.UnitSquareMesh(10, 10)

top = Interface(y=1.0)
bottom = Interface(y=0.0)
interface = Interface(y=0.7)

lower = Region(y=0.7)

# NOTE: the bug illustrated in this example does _not_ appear if a "crossed"
#       mesh is not used. It is, however, present when an unstructured mesh
#       (e.g. as generated by GMesh) is used.
mesh = fenics.UnitSquareMesh(10, 10, "crossed")

subdomains = fenics.MeshFunction("size_t", mesh, 2)
subdomains.set_all(1)
lower.mark(subdomains, 0)

facet_subdomains = fenics.MeshFunction("size_t", mesh, 1)
facet_subdomains.set_all(0)
bottom.mark(facet_subdomains, 1)
interface.mark(facet_subdomains, 2)
top.mark(facet_subdomains, 3)


fenics.plot(mesh)


# ## Set up problem

# We are solving the Poisson equation where the only charge is on the line dividing the 2 subdomains.


# Relative permittivities of regions (unitless)
eps_a = 15
eps_b = 25
# Voltage boundary condition, units: V
U = 0.0
# signed number density, units: um^-2
rho = -1000
# conversion factor q_e / eps0, units 'V um'
qe_over_eps0 = 0.01809512817972783  


V = fenics.FunctionSpace(mesh, "CG", 1)
u, v = fenics.TrialFunction(V), fenics.TestFunction(V)
dx = fenics.Measure("cell", domain=mesh, subdomain_data=subdomains)
ds = fenics.Measure("interior_facet", domain=mesh, subdomain_data=facet_subdomains)


# Units 'V^2 um^-2 * um^3' = V^2 um
a = (
    + eps_a * fenics.dot(fenics.grad(u), fenics.grad(v)) * dx(0)
    + eps_b * fenics.dot(fenics.grad(u), fenics.grad(v)) * dx(1)
)
# Units 'um^-2 * V um * V * um^2' = V^2 um
L = rho * qe_over_eps0 * fenics.avg(v) * ds(2)

equation = (a == L)


bcs = [
    fenics.DirichletBC(V, U, facet_subdomains, 1),
    fenics.DirichletBC(V, U, facet_subdomains, 3),
]


# ## Define solvers

# We will assemble the system in 3 different ways to compare:
# 1. `fenics.assemble_system`
# 2. `fenics.assemble` followed by `bc.apply` (non-symmetric BC application)
# 3. `fenics.assemble` followed by `bc.zero_columns` (symmetric BC application, AFAIK)


def assemble_system(equation, bcs):
    return fenics.assemble_system(equation.lhs, equation.rhs, bcs)


def assemble_and_apply(equation, bcs):
    A = fenics.assemble(equation.lhs)
    b = fenics.assemble(equation.rhs)
    for bc in bcs:
        bc.apply(A, b)
    return A, b


def assemble_and_zero_columns(equation, bcs):
    A = fenics.assemble(equation.lhs)
    b = fenics.assemble(equation.rhs)
    # AFAIK this applies the boundary conditions in such a manner that
    # the symmetry of A is preserved.
    for bc in bcs:
        bc.zero_columns(A, b, 1.0)
    return A, b


def solve(V, equation, bcs, assemble):
    A, b = assemble(equation, bcs)
    solution = fenics.Function(V)
    fenics.la.solver.solve(A, solution.vector(), b, "gmres", "ilu")
    return solution


# ## Solve using different assembly methods


fenics.parameters["krylov_solver"]["relative_tolerance"] = 1e-9


sol_assemble_system = solve(V, equation, bcs, assemble_system)

sol_apply = solve(V, equation, bcs, assemble_and_apply)

sol_zero_columns = solve(V, equation, bcs, assemble_and_zero_columns)


# ## Compare solutions


def evaluate_vectorized(f, xs):
    return [f(x) for x in xs]



ys = np.hstack([np.linspace(0, 0.7, 700, endpoint=False), np.linspace(0.7, 1.0, 300)])
xs = np.zeros_like(ys)
coords = np.array([xs, ys]).T

plt.plot(ys, evaluate_vectorized(sol_assemble_system, coords), label="assemble_system")
plt.plot(ys, evaluate_vectorized(sol_apply, coords), "--", label="assemble_and_apply")
plt.plot(ys, evaluate_vectorized(sol_zero_columns, coords), ":", label="assemble_and_zero_columns")

plt.xlabel("Y position [um]")
plt.ylabel("potential [V]")

plt.legend()


# We see that `assemble_and_apply` and `assemble_and_zero_columns` agree with each other, but `assemble_system` produces a different result.
# 
# The results look even wilder if we set the boundary conditions to be different from 0.

‌

2021-04-06T00:29:43+00:00

Comments (2)