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:
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) class Right(SubDomain): def inside(self, x, on_boundary): return on_boundary and near(x, len_x) class Top(SubDomain): def inside(self, x, on_boundary): return on_boundary and near(x, len_y) class Bottom(SubDomain): def inside(self, x, on_boundary): return on_boundary and near(x, 0.0) class Wall(SubDomain): def inside(self, x, on_boundary): return x <= wall_thickness class Interface(SubDomain): def inside(self, x, on_boundary): return near(x, 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)