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

5

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)
1 Like