Error in Result at boundary for a Poisson Subdomains/Interior boundary problem

Hello there,

I am currently trying to solve a simple Poisson example with subdomains, following the answer provided to me in this post [Electrostatics Problem with Surface Charge]. He solves it with superposition in that example, which i cannot do in my application.

So, i attempted to solve it in the way shown in my MWI below: (It is a simple example with a square domain with some of it being a subdomain with different properties)

The code solves fine, however, I am getting a strange error on the top boundary, see these images:

The graph on the right is taken down the vertical central axis. Everything seems to be correct aside from the obvious kink at the top boundary (You can see it quite clearly oscillating if you zoom into the colour plot). It acts the same even if I set \sigma = 0 (No surface charge). If I delete the entire surface integral \int \sigma \ dS (ignoring surface charge altogether), then the oscillations at the top also disappear, and the solution (of a perfectly linear decrease in potential) is correct, but I cant seem to find out why the interior facet integral is affecting the solution in such a way.

If anybody could test my code and see where I have went wrong/what is wrong that would be really appreciated, the code is below:
Thanks,

from fenics import *
from dolfin import *
from matplotlib import pyplot as plt
plt.style.use('classic')

# Classes for boundaries and subdomains
# Create classes defining the geometry
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], len_x)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], len_y)
    
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= wall_thickness
    
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], wall_thickness)

left = Left()
right = Right()
top = Top()
bottom = Bottom()
wall = Wall()
interface = Interface()
    
# Constants
qe = Constant(1.602e-19)
e0 = Constant(8.854e-12)

# Domain dimensions
len_x = 4e-2
len_y = 4e-2
wall_thickness = 1e-2

# Material property in domains
e_wall = Constant(3)
e_air = Constant(1)

# Value on the interface
sigma = Constant(1e14)

# Dirichlet boundary values
u_top = 100e+3
u_bottom = 0

# Generate a mesh
mesh = RectangleMesh(Point(0,0), Point(len_x,len_y), 100, 100, "right/left")

# Meshfunctions for domains and boundaries
domains = MeshFunction('size_t',mesh,mesh.topology().dim())
boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
domains.set_all(0)
boundaries.set_all(0)

# Set values
wall.mark(domains,1)
top.mark(boundaries,1)
bottom.mark(boundaries,2)
left.mark(boundaries,3)
right.mark(boundaries,4)
interface.mark(boundaries,5)

# Measures
dx = Measure('dx',domain=mesh,subdomain_data=domains)
ds = Measure('ds',domain=mesh,subdomain_data=boundaries)
dS = Measure('dS',domain=mesh,subdomain_data=boundaries)

# Set up functions and functionspace
V = FunctionSpace(mesh,'Lagrange',1)
u = Function(V)
w = TestFunction(V)

# Boundary conditions
bc1 = DirichletBC(V,u_top,boundaries,1)
bc2 = DirichletBC(V,u_bottom,boundaries,2)
bcs = [bc1,bc2]

# Variational form
F1 = (1/qe)*dot((e0*e_air)*grad(u),grad(w))*dx(0) - sigma*w('-')*dS(5)
F2 = (1/qe)*dot((e0*e_wall)*grad(u),grad(w))*dx(1) - sigma*w('-')*dS(5)
F = F1 + F2

# Set up solver and solve
J = derivative(F,u)
problem = NonlinearVariationalProblem(F,u,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['maximum_iterations'] = 15
prm['snes_solver']['method'] = 'newtonls'
prm['snes_solver']['linear_solver'] = 'gmres'
prm['snes_solver']['preconditioner'] = 'hypre_amg'
prm['snes_solver']['relative_tolerance'] = 1e-6
solver.solve()

# Plot u
plot(u)
# Save filed as XDMF
u.rename('potential','potential')
file_u = XDMFFile('surface_charge_test/solution.xdmf')
file_u.write(u)

Hi @timwong,

As far as I can tell, you haven’t gone wrong anywhere in your code. To me, it seems like a bug in FEniCS. This was exactly the reason that I used superposition in the post you linked. When I test your code, I get the same issue that you do.

I can confirm that I get exactly this behavior:

Hello @conpierce8,

Thanks for confirming this behaviour.

I still haven’t found a fix for the problem above, but strangely enough, it appears to have something to do with using the NonlinearVariationalProblem and solver, because if you attempt to solve it just as a linear problem using:

a = (1/qe)*(e0*epsilon)*dot(grad(u),grad(w))*dx
L = sigma*w('-')*dS(5)

solve(a == L, u, bcs)

It gives the correct results and the oscillations disappear.

However, for my application, I am solving the Poisson equation coupled to a bunch of other non-linear equations, which is why I choose the form of F=0 instead.

Any ideas what could be done?

Thanks,

It looks like the calculation is confused about the boundaries. Try making your boundaries strictly unique. At the moment a subset of your top boundary ( near(x[1], len_y) also satisfies your wall (x[0] <= wall_thickness). Try making them strictly exclusive. My guess is you’d want top boundary to be “inclusive”, in which case I’d try tightening the definition of Wall (and Interface).

1 Like

This is interesting. I didn’t realize the linear problem could be solved correctly (I thought this failed for linear and nonlinear alike since the erroneous solution I encountered previously was for a linear problem). It looks like the error might be related to this 2-year-old issue reported on Bitbucket: https://bitbucket.org/fenics-project/dolfin/issues/1075/assembly-of-integrals-over-interal. In that bug report, it is suggested that SystemAssembler (used by NonlinearVariationalSolver) is the culprit. Indeed, I have created a code example here demonstrating this–it solves correctly when using solve(a == L, ...) or when using assemble followed by bc.apply(), but gives an incorrect solution when using assemble_system. @dparsons: This example also shows that this issue persists even when the charged surface is strictly exclusive of the boundaries with DirichletBCs.

Happily, inspired by this discussion, I think I can provide a solution/workaround. It seems that calling bc.apply(u) before creating the NonlinearVariationalProblem solves the issue, although I’m not sure why.

from fenics import *
from dolfin import *
from matplotlib import pyplot as plt
plt.style.use('classic')

# Classes for boundaries and subdomains
# Create classes defining the geometry
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], len_x)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], len_y)
    
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= wall_thickness
    
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], wall_thickness)

left = Left()
right = Right()
top = Top()
bottom = Bottom()
wall = Wall()
interface = Interface()
    
# Constants
qe = Constant(1.602e-19)
e0 = Constant(8.854e-12)

# Domain dimensions
len_x = 4e-2
len_y = 4e-2
wall_thickness = 1e-2

# Material property in domains
e_wall = Constant(3)
e_air = Constant(1)

# Value on the interface
sigma = Constant(1e14)

# Dirichlet boundary values
u_top = 100e+3
u_bottom = 0

# Generate a mesh
mesh = RectangleMesh(Point(0,0), Point(len_x,len_y), 100, 100, "right/left")

# Meshfunctions for domains and boundaries
domains = MeshFunction('size_t',mesh,mesh.topology().dim())
boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
domains.set_all(0)
boundaries.set_all(0)

# Set values
wall.mark(domains,1)
top.mark(boundaries,1)
bottom.mark(boundaries,2)
left.mark(boundaries,3)
right.mark(boundaries,4)
interface.mark(boundaries,5)

# Measures
dx = Measure('dx',domain=mesh,subdomain_data=domains)
ds = Measure('ds',domain=mesh,subdomain_data=boundaries)
dS = Measure('dS',domain=mesh,subdomain_data=boundaries)

# Set up functions and functionspace
V = FunctionSpace(mesh,'Lagrange',1)
u = Function(V)
w = TestFunction(V)

# Boundary conditions
bc1 = DirichletBC(V,u_top,boundaries,1)
bc2 = DirichletBC(V,u_bottom,boundaries,2)
bcs = [bc1,bc2]

for b in bcs:
    b.apply(u.vector())

# Variational form
F1 = (1/qe)*dot((e0*e_air)*grad(u),grad(w))*dx(0) - sigma*w('-')*dS(5)
F2 = (1/qe)*dot((e0*e_wall)*grad(u),grad(w))*dx(1) - sigma*w('-')*dS(5)
F = F1 + F2

# Set up solver and solve
J = derivative(F,u)
problem = NonlinearVariationalProblem(F,u,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['maximum_iterations'] = 15
prm['snes_solver']['method'] = 'newtonls'
prm['snes_solver']['linear_solver'] = 'gmres'
prm['snes_solver']['preconditioner'] = 'hypre_amg'
prm['snes_solver']['relative_tolerance'] = 1e-6
solver.solve()

# Plot u
plot(u)
# Save filed as XDMF
u.rename('potential','potential')
file_u = XDMFFile('surface_charge_test/solution.xdmf')
file_u.write(u)

@timwong

On a side note, your variational formulation may be incorrect. I’m not 100% sure, but I think that including - sigma*w('-')*dS(5) in both F1 and F2 causes you to double-count the surface charge.

@conpierce8

That’s interesting - and you’re right, I tried what @dparsons suggested and it yielded the same error.

However, your workaround appears to work for my application, though, another problem seems to have arisen. It may be a separate issue from the above, so let me know if I should make a new post. If by chance it is related though, its regarding running the above code in parallel.

So with some modifications, the quantity which I am interested in is the electric field, E, which i calculate in the modified code below as E = -grad(u) and project it onto a Vector space to save to paraview. However, when I run this code with mpirun I get the following:

MWI_example3_ZOOM

The second image is the first, but zoomed in to the bottom corner. The boundaries of the mesh which are distributed over each process are clearly visible, and the field experiences some fluctuations and artefacts. I tried to implement the same fix in my larger code, which is a transient problem and the same issue drove the entire problem to nonconvergence due to these unphysical boundary lines.

The modified code of the above MWE is below, any ideas? (Some things like the solver I have also changed just to isolate variables, it appears the error is independent of the solver (SNES or newton), the ghost mode (shared_vertex or shared_facet) or the partitioner (SCOTCH or ParMETIS))

Update: Further testing shows that this is also exclusive to using NonlinearVariationalProblem and solver - these boundary lines do not appear when simply solving a == L.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jul  9 12:55:01 2021

@author: timothy
"""
from fenics import *
from dolfin import *
from matplotlib import pyplot as plt
plt.style.use('classic')

parameters["ghost_mode"] = "shared_facet"
# Classes for boundaries and subdomains
# Create classes defining the geometry
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], len_x)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], len_y)
    
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= wall_thickness
    
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], wall_thickness)

left = Left()
right = Right()
top = Top()
bottom = Bottom()
wall = Wall()
interface = Interface()
    
# Constants
qe = Constant(1.602e-19)
e0 = Constant(8.854e-12)

# Domain dimensions
len_x = 4e-2
len_y = 4e-2
wall_thickness = 1e-2

# Material property in domains
e_wall = Constant(3)
e_air = Constant(1)

# Value on the interface
sigma = Constant(2)

# Dirichlet boundary values
u_top = 100e+3
u_bottom = 0

# Generate a mesh
mesh = RectangleMesh(Point(0,0), Point(len_x,len_y), 500, 500, "right/left")

# Meshfunctions for domains and boundaries
domains = MeshFunction('size_t',mesh,mesh.topology().dim())
boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
domains.set_all(0)
boundaries.set_all(0)

# Set values
wall.mark(domains,1)
top.mark(boundaries,1)
bottom.mark(boundaries,2)
left.mark(boundaries,3)
right.mark(boundaries,4)
interface.mark(boundaries,5)

# Measures
dx = Measure('dx',domain=mesh,subdomain_data=domains)
ds = Measure('ds',domain=mesh,subdomain_data=boundaries)
dS = Measure('dS',domain=mesh,subdomain_data=boundaries)

# Set up functions and functionspace
V = FunctionSpace(mesh,'Lagrange',1)
u = Function(V)
w = TestFunction(V)

# Boundary conditions
bc1 = DirichletBC(V,u_top,boundaries,1)
bc2 = DirichletBC(V,u_bottom,boundaries,2)
bcs = [bc1,bc2]

for b in bcs:
    b.apply(u.vector())

# Variational form
F1 = dot((e0*e_air)*grad(u),grad(w))*dx(0) - sigma*w('-')*dS(5)
F2 = dot((e0*e_wall)*grad(u),grad(w))*dx(1)
F = F1 + F2

# Set up solver and solve
J = derivative(F,u)
problem = NonlinearVariationalProblem(F,u,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['maximum_iterations'] = 15
prm['snes_solver']['method'] = 'newtonls'
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'hypre_amg'
prm['newton_solver']['relative_tolerance'] = 1e-6
solver.solve()

V_e = VectorFunctionSpace(mesh,'Lagrange',1)
E = project(-grad(u),V_e, solver_type = 'gmres', preconditioner_type = 'hypre_amg')

# Save files as XDMF
u.rename('potential','potential')
E.rename('field','field')
file_u = XDMFFile(MPI.comm_world,'surface_charge_test/solution.xdmf')
file_E = XDMFFile(MPI.comm_world,'surface_charge_test/field.xdmf')
file_u.parameters["flush_output"] = True
file_E.parameters["flush_output"] = True
file_u.write(u)
file_E.write(E)

P.S. @conpierce8 Thanks for the heads up about the double surface charge, I have changed it now :slight_smile:

Hi @timwong,

Disclaimer: I’m a novice at nonlinear problems and parallel solving. :slightly_smiling_face:

Having said that, I played around with your problem for a bit. First, when I visualized the solution (from “solution.xdmf”) in Paraview with a “Warp By Scalar” filter, I noticed that the potential had irregularities along the mesh partition boundaries. (Hence the irregularities in the field.) The Warp By Scalar filter helped to amplify irregularities, which were subtle. See below—mesh partitioning and solution for the potential.

solution

These irregularities, as you can see, always occur on the left of the internal boundary, and never occur on the right. This asymmetry was unaffected by changing the size of the wall and by changing the relative values of e_wall and e_air (I tried simulating with e_wall = 1 and e_air = 3). But this got me wondering if the “asymmetry” in the definition of the variational form was causing the issue, since the surface integral dS appears in F1 but not in F2.

So I defined an Expression called epsilon to evaluate the permittivity anywhere in the domain, which allows me to reduced the variational form to a single volume integral over dx (rather than two separate integrals over dx(0) and dx(1)), and changed the variational form to

F = dot((e0*epsilon)*grad(u),grad(w))*dx - sigma*w('-')*dS(5)

This seems to have solved the issue on my machine when using 8 processors; don’t ask me to explain why :slightly_smiling_face:. See solution below for the potential and electric field.
potential
field

Here’s the code:

from fenics import *
from dolfin import *
from matplotlib import pyplot as plt
plt.style.use('classic')

parameters["ghost_mode"] = "shared_facet"
# Classes for boundaries and subdomains
# Create classes defining the geometry
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0)
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], len_x)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], len_y)
    
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= wall_thickness
    
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], wall_thickness)

left = Left()
right = Right()
top = Top()
bottom = Bottom()
wall = Wall()
interface = Interface()
    
# Constants
qe = Constant(1.602e-19)
e0 = Constant(8.854e-12)

# Domain dimensions
len_x = 4e-2
len_y = 4e-2
wall_thickness = 1e-2

# Material property in domains
e_wall = Constant(3)
e_air = Constant(1)
epsilon = Expression('(x[0]<wall_x)?e_wall:e_air',
                      degree=2,
                      wall_x=wall_thickness+DOLFIN_EPS,
                      e_wall=e_wall, e_air=e_air)

# Value on the interface
sigma = Constant(2)

# Dirichlet boundary values
u_top = 100e+3
u_bottom = 0

# Generate a mesh
mesh = RectangleMesh(Point(0,0), Point(len_x,len_y), 500, 500, "right/left")

# Meshfunctions for domains and boundaries
domains = MeshFunction('size_t',mesh,mesh.topology().dim())
boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
domains.set_all(0)
boundaries.set_all(0)

# Set values
wall.mark(domains,1)
top.mark(boundaries,1)
bottom.mark(boundaries,2)
left.mark(boundaries,3)
right.mark(boundaries,4)
interface.mark(boundaries,5)

# Measures
dx = Measure('dx',domain=mesh,subdomain_data=domains)
ds = Measure('ds',domain=mesh,subdomain_data=boundaries)
dS = Measure('dS',domain=mesh,subdomain_data=boundaries)

# Set up functions and functionspace
V = FunctionSpace(mesh,'Lagrange',1)
u = Function(V)
w = TestFunction(V)

# Boundary conditions
bc1 = DirichletBC(V,u_top,boundaries,1)
bc2 = DirichletBC(V,u_bottom,boundaries,2)
bcs = [bc1,bc2]

for b in bcs:
    b.apply(u.vector())

# Variational form
F = dot((e0*epsilon)*grad(u),grad(w))*dx - sigma*w('-')*dS(5)

# Set up solver and solve
J = derivative(F,u)
problem = NonlinearVariationalProblem(F,u,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['maximum_iterations'] = 15
prm['snes_solver']['method'] = 'newtonls'
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'hypre_amg'
prm['newton_solver']['relative_tolerance'] = 1e-6
solver.solve()

V_e = VectorFunctionSpace(mesh,'Lagrange',1)
E = project(-grad(u),V_e, solver_type = 'gmres', preconditioner_type = 'hypre_amg')

# Save files as XDMF
u.rename('potential','potential')
E.rename('field','field')
file_u = XDMFFile(MPI.comm_world,'surface_charge_test/solution.xdmf')
file_E = XDMFFile(MPI.comm_world,'surface_charge_test/field.xdmf')
file_u.parameters["flush_output"] = True
file_E.parameters["flush_output"] = True
file_u.write(u)
file_E.write(E)
1 Like

Hello @conpierce8,

Thank you so much! That seems to have fixed it for me too even for >8 MPI processes, even if the reasons why remain a mystery. Thank you also for going through your debugging thought process, you have been incredibly helpful :slight_smile:

Thanks again,

1 Like