Hello, i use fenicsx 0.7.2 and i am facing issue in convergence for non linear system with dirichlet and neumann boundary condition.
One of these neumann bc is dependant on uh the current solution hence non linearity.
the system is the one from Gate-induced carrier density modulation in bulk graphene: theories and electrostatic simulation using Matlab pdetool | Journal of Computational Electronics figure 4 :
The system a big rectangle with 4 Neumann BC, 3 of them with g=0, and the bottom one with g = - constant * sign(u) * u^2. Inside the big rectangle there is a smaller one, that doesn’t touch the bigger one, it has Dirichlet BC u=1.5. Objectiv is to get the potential in the space (epsilon =1) between the two rectangles.
I recreated this system in matlab pdetool, and got the same result as the article, but i need to do it in python since i have other script in this language.
- I created via gmsh the system.
- Using tutoriels for dirichlet/neumann/robin i created a function for Neumann BC. I also used the tutoriel for non linear poisson for newton solver.
- The system converge after 6 iterations if i put rtol = 1e-2, and the potentiel map looks smooth and okay even if the value of u is not exactly what it should be at y=0, see image :
- Now if i go beyond that, rtol = 1e-3, it diverge/converge and so on, never reaching 1e-3, and it hit the limit of 50 iterations.
Can you please tell me if the implementation of Neumann BC, the creation of F and the call of Newton solver looks good?
if needed i can provide the whole code and gmsh code, it’s just too long for posting as a first question (lot’s of function defined), i suppose.
Neumann Fct :
def set_neumann_bc(self, markers, value, test_fct):
if type(markers)==int :
markers=[markers]
facet_indices, facet_markers = [], []
for marker in markers:
facets = self.mesh_boundaries_tags.find(marker)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(self.mesh_geom, self.mesh_geom.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
dss = ufl.Measure("ds", domain=self.mesh_geom, subdomain_id = facet_markers[0], subdomain_data = facet_tag)
Formul = inner(value, test_fct) * dss
return Formul
adding Neumann to F :
F= self.dielectric_cons * inner(grad(self.uh), grad(self.v)) * dx # f=0 in air
###Apply Dirichlet boundary conditions where needed using markers
self.set_dirichlet_bc(self.markers_1d["square"], bias)
####### Neumann 2dmat
g2d = -13.295*pow(self.uh,2)*ufl.sign(self.uh)
F+=self.set_neumann_bc(self.markers_1d["mat2d"],g2d,self.v)
self.problem = NonlinearProblem(F,self.uh, bcs=self.bcs)
Non linear Newton solver call :
solver = NewtonSolver(MPI.COMM_WORLD, to_solve.problem)
# solver.convergence_criterion = "residual"
solver.rtol = 1e-2
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}setInitialGuessNonzero"] = "True"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(to_solve.uh)
assert (converged)
print(f"Number of iterations: {n:d}")
Thanks for reading