Advice on: Newton solver does not converge

variational formulation
1

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)[0]
    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.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?
3

I am having a similar issue as I posted on

Have you been able to figure out what the issue with your code is?

4

I found out some things. It depends a lot on the initial guess. Depending on the initial guess I have seen the following things:

  • Does not converge
  • Converges to the physical solution
  • Converges to a non physical solution

So yeah it is pretty fragile and a better problem discretization would be nice.

5

I’ve posted some tips and tricks for nonlinear problems with a Newton solver here.

1 Like