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)
