I would like to solve the following system:
-∇ϕ = E div(E) = ρ 0 = ∂ₜρ = div(μEρ) + s
It can be interpreted as the steady state of charged particle density ρ, generated by a process s
and moving around under under an external potential ϕ as well as its own electric field.
I tried the following:
from fenics import * mesh = UnitIntervalMesh(5) class ManufacturedSolution(): x = SpatialCoordinate(mesh) phi = x**2 E = -grad(phi) rho = div(E) mu = 2 s = div(E*mu*rho) truth = ManufacturedSolution() cell = mesh.ufl_cell() deg = 1 fe_rho = FiniteElement("CG", cell, deg) fe_E = VectorElement("CG", cell, deg) fe_phi = FiniteElement("CG", cell, deg) fe = MixedElement([fe_rho, fe_E, fe_phi]) V = FunctionSpace(mesh, fe) v_rho, v_E, v_phi = TestFunctions(V) u = Function(V) rho, E, phi = split(u) V_rho = V.sub(0) V_E = V.sub(1) V_phi = V.sub(2) class LeftBoundary(SubDomain): def inside(self, x, on_boundary): return on_boundary and near(x, 0.0) bc_rho = DirichletBC(V_rho, project(truth.rho, V_rho.collapse()), LeftBoundary()) bc_phi = DirichletBC(V_phi, project(truth.phi, V_phi.collapse()), "on_boundary") bcs=[bc_rho, bc_phi] F_rho = inner(truth.s - div(truth.mu*E*rho), v_rho) * dx F_E = inner(div(E) - rho, v_phi)*dx F_phi = inner(-grad(phi) - E, v_E) * dx F = F_rho + F_phi + F_E def relerror(phi_exact, phi): scale = norm(phi_exact, "L2") return errornorm(phi_exact, phi, "L2")/scale solve(F == 0, u, bcs=bcs) rho,E, phi = u.split() print(relerror(project(truth.rho, V_rho.collapse()), rho)) print(relerror(project(truth.E , V_E .collapse()), E )) print(relerror(project(truth.phi, V_phi.collapse()), phi)) plot(rho)
But I get
*** Error: Unable to solve nonlinear system with NewtonSolver. *** Reason: Newton solver did not converge because maximum number of iterations reached. *** Where: This error was encountered inside NewtonSolver.cpp. *** Process: 0 No Jacobian form specified for nonlinear variational problem. Differentiating residual form F to obtain Jacobian J = F'. Solving nonlinear variational problem. Newton iteration 0: r (abs) = 4.045e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09) Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09) ...
- I am looking for both being able to solve the issue for the specific example at hand and some general advice.
- Why does the solver not converge?
- Is the system instable? In the linear case, I would look at the assembled matrix or the condition number. Here I don’t really know how to start investigating the problem.
- How to reformulate the problem in a way that can be solved?