About NewtonSolver error in Dolfin

AOlz · October 29, 2025, 7:03pm

Please can you help resolve this error in my code:
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"

Code:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

==================== Nondimensional parameters ====================

alpha = 1.5
beta = 0.75
mu_phi = -np.pi
gamma_vals = [0.0, 0.1] # line tension (γ̄)
f_v = 1e-3 # traction force
kbar = 1.0 # bending rigidity
rho = 10.0 # domain radius
dt = 1e-3 # nondimensional time increment
times = [0.0, 0.91, 1.45] # times for profiles
Ns = 800
epsr = 5e-4
delta = 1e-4 # asymptotic origin shift

==================== Mesh and FunctionSpace ====================

mesh1d = IntervalMesh(Ns-1, delta, rho)

2nd-order Lagrange element and 7-field mixed element

P2 = FiniteElement(“CG”, “interval”, 2)
mixed_elem = MixedElement([P2] * 7)
W = FunctionSpace(mesh1d, mixed_elem)
print(“ Mixed space created — DOFs:”, W.dim())

Unknowns [r, z, ψ, M, L, a, λ] and test functions

Y = Function(W)
v = TestFunction(W)
r, z, psi, M, L, a, lam = split(Y)
vr, vz, vpsi, vM, vL, va, vlam = split(v)
x = SpatialCoordinate(mesh1d)[0]

==================== Protein distribution (σ̄ = exp(σ̃)) ====================

Vsig = FunctionSpace(mesh1d, “CG”, 1)
sigma_tilde = Function(Vsig)
sigma_tilde.vector()[:] = 0.0
sigma_bar = exp(sigma_tilde)
C = mu_phi * sigma_bar
r_safe = sqrt(r2 + epsr2)

==================== Governing equations====================

dr = r.dx(0)
dz = z.dx(0)
dpsi = 2*(M/kbar) - sin(psi)/r_safe + 2C
dM = M.dx(0)
L_def = L - rdM

Fv = (alphasigma_bar - beta)*2 + kbar(M/kbar - C)**2 + lam
- (2M/kbar - sin(psi)/r_safe + 2*C)*M
gamma_sym = Constant(0.0)

dL = r*((Fv - gamma_symcos(psi)/r_safe) *
(2M/kbar - sin(psi)/r_safe + 2C)
+ (Fv - gamma_symcos(psi)/r_safe)sin(psi)/r_safe
+ 2M*((M/kbar) - sin(psi)/r_safe + C)*sin(psi)/r_safe)

da = 2np.pir
dsigma = sigma_bar.dx(0)
dlam = 2*(mu_phiM - alpha(alpha*sigma_bar - beta))*dsigma

Weak form (residual)

F = (
(dr - cos(psi))*vr

(dz - sin(psi))*vz
(psi.dx(0) - dpsi)*vpsi
L_def*vM
(L.dx(0) - dL)*vL
(a.dx(0) - da)*va
(lam.dx(0) - dlam)*vlam
) * dx

==================== Boundary Condition====================

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], delta)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], rho)

left = Left()
right = Right()

Apply BCs

bc_r0 = DirichletBC(W.sub(0), Constant(0.0), left) # r(0) = 0
bc_psi0 = DirichletBC(W.sub(2), Constant(0.0), left) # ψ(0) = 0
bc_L0 = DirichletBC(W.sub(4), Constant(0.0), left) # L(0) = 0
bc_a0 = DirichletBC(W.sub(5), Constant(0.0), left) # a(0) = 0

bc_zN = DirichletBC(W.sub(1), Constant(0.0), right) # z(ρ) = 0
bc_psiN = DirichletBC(W.sub(2), Constant(0.0), right) # ψ(ρ) = 0
bc_LN = DirichletBC(W.sub(4), Constant(0.0), right) # L(ρ) = 0
bc_aN = DirichletBC(W.sub(5), Constant(np.pi*rho**2), right) # a(ρ) = πρ²

bcs = [bc_r0, bc_psi0, bc_L0, bc_a0, bc_zN, bc_psiN, bc_LN, bc_aN]

==================== Jacobian and Nonlinear solver ====================

Automatic symbolic Jacobian of the residual

J = derivative(F, Y)

problem = NonlinearVariationalProblem(F, Y, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm[“newton_solver”][“relative_tolerance”] = 1e-6
prm[“newton_solver”][“absolute_tolerance”] = 1e-8
prm[“newton_solver”][“maximum_iterations”] = 30

==================== Time stepping and results ====================

results = {}
for gamma in gamma_vals:
gamma_sym.assign(gamma)
print(f"\n=== Solving for γ̄ = {gamma} ===“)
results[gamma] =
for t in times:
print(f” t̄ = {t}")
solver.solve()
r_vals = Y.sub(0).vector().get_local()
z_vals = Y.sub(1).vector().get_local()
results[gamma].append((r_vals.copy(), z_vals.copy()))

    # ---- Protein diffusion update (forward Euler) ----
    s_nodes = mesh1d.coordinates().flatten()
    sigma_old = np.exp(sigma_tilde.vector().get_local())
    dsigma_num = np.gradient(sigma_old, s_nodes)
    M_arr = Y.sub(3).vector().get_local()
    dM_num = np.gradient(M_arr, s_nodes)
    flux_inner = alpha**2*sigma_old*dsigma_num - mu_phi*dM_num
    div_term = np.gradient(s_nodes*flux_inner, s_nodes)
    rhs = 2.0/(np.sqrt(s_nodes**2 + epsr**2)*sigma_old)*div_term
    sigma_new = np.maximum(sigma_old + dt*rhs, 1e-12)
    sigma_tilde.vector()[:] = np.log(sigma_new)
    sigma_tilde.vector()[-1] = np.log(beta/alpha)

==================== Axisymmetric surface plots====================

theta = np.linspace(0, 2np.pi, 200)
for gamma, shapes in results.items():
fig = plt.figure(figsize=(10, 4))
for (r_vals, z_vals), t in zip(shapes, times):
R, TH = np.meshgrid(r_vals, theta)
X = Rnp.cos(TH)
Yp = R*np.sin(TH)
Z = np.tile(z_vals, (len(theta), 1))
ax = fig.add_subplot(1, len(times), list(times).index(t)+1, projection=‘3d’)
ax.plot_surface(X, Yp, Z, color=‘salmon’, alpha=0.9, edgecolor=‘none’)
ax.set_title(f"γ̄={gamma}, t̄={t}“)
ax.set_axis_off()
plt.suptitle(f"Membrane Budding Evolution (γ̄={gamma})”)
plt.show()

dokken · October 29, 2025, 7:49pm

Please format your code with 3x` encapsulation, i.e.

```python
#add code here
```

AOlz · October 29, 2025, 8:07pm

Thank you for your kind response.

# -------------------------------------------------------------
# Budding Mechanics
# By A.Oluzeph

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

# ==================== Nondimensional parameters ====================
alpha = 1.5
beta = 0.75
mu_phi = -np.pi
gamma_vals = [0.0, 0.1]         # line tension (γ̄)
f_v = 1e-3                      # traction force
kbar = 1.0                      # bending rigidity
rho = 10.0                      # domain radius
dt = 1e-3                       # nondimensional time increment
times = [0.0, 0.91, 1.45]       # times for profiles
Ns = 800
epsr = 5e-4
delta = 1e-4                    # asymptotic origin shift

# ==================== Mesh and FunctionSpace ====================
mesh1d = IntervalMesh(Ns-1, delta, rho)

# 2nd-order Lagrange element and 7-field mixed element
P2 = FiniteElement("CG", "interval", 2)
mixed_elem = MixedElement([P2] * 7)
W = FunctionSpace(mesh1d, mixed_elem)
print("✅ Mixed space created — DOFs:", W.dim())

# Unknowns [r, z, ψ, M, L, a, λ] and test functions
Y = Function(W)
v = TestFunction(W)
r, z, psi, M, L, a, lam = split(Y)
vr, vz, vpsi, vM, vL, va, vlam = split(v)
x = SpatialCoordinate(mesh1d)[0]

# ==================== Protein distribution (σ̄ = exp(σ̃)) ====================
Vsig = FunctionSpace(mesh1d, "CG", 1)
sigma_tilde = Function(Vsig)
sigma_tilde.vector()[:] = 0.0
sigma_bar = exp(sigma_tilde)
C = mu_phi * sigma_bar
r_safe = sqrt(r**2 + epsr**2)

# ==================== Governing equations====================
dr = r.dx(0)
dz = z.dx(0)
dpsi = 2*(M/kbar) - sin(psi)/r_safe + 2*C
dM = M.dx(0)
L_def = L - r*dM

Fv = (alpha*sigma_bar - beta)**2 + kbar*(M/kbar - C)**2 + lam \
     - (2*M/kbar - sin(psi)/r_safe + 2*C)*M
gamma_sym = Constant(0.0)

dL = r*((Fv - gamma_sym*cos(psi)/r_safe) *
        (2*M/kbar - sin(psi)/r_safe + 2*C)
        + (Fv - gamma_sym*cos(psi)/r_safe)*sin(psi)/r_safe
        + 2*M*((M/kbar) - sin(psi)/r_safe + C)*sin(psi)/r_safe)

da = 2*np.pi*r
dsigma = sigma_bar.dx(0)
dlam = 2*(mu_phi*M - alpha*(alpha*sigma_bar - beta))*dsigma

# Weak form (residual)
F = (
    (dr - cos(psi))*vr
  + (dz - sin(psi))*vz
  + (psi.dx(0) - dpsi)*vpsi
  + L_def*vM
  + (L.dx(0) - dL)*vL
  + (a.dx(0) - da)*va
  + (lam.dx(0) - dlam)*vlam
) * dx

# ==================== Boundary Condition====================
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], delta)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], rho)

left = Left()
right = Right()

# Apply BCs
bc_r0   = DirichletBC(W.sub(0), Constant(0.0), left)               # r(0) = 0
bc_psi0 = DirichletBC(W.sub(2), Constant(0.0), left)               # ψ(0) = 0
bc_L0   = DirichletBC(W.sub(4), Constant(0.0), left)               # L(0) = 0
bc_a0   = DirichletBC(W.sub(5), Constant(0.0), left)               # a(0) = 0

bc_zN   = DirichletBC(W.sub(1), Constant(0.0), right)              # z(ρ) = 0
bc_psiN = DirichletBC(W.sub(2), Constant(0.0), right)              # ψ(ρ) = 0
bc_LN   = DirichletBC(W.sub(4), Constant(0.0), right)              # L(ρ) = 0
bc_aN   = DirichletBC(W.sub(5), Constant(np.pi*rho**2), right)     # a(ρ) = πρ²

bcs = [bc_r0, bc_psi0, bc_L0, bc_a0, bc_zN, bc_psiN, bc_LN, bc_aN]

# ==================== Jacobian and Nonlinear solver ====================
# Automatic symbolic Jacobian of the residual
J = derivative(F, Y)

problem = NonlinearVariationalProblem(F, Y, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["relative_tolerance"] = 1e-6
prm["newton_solver"]["absolute_tolerance"] = 1e-8
prm["newton_solver"]["maximum_iterations"] = 30

# ==================== Time stepping and results ====================
results = {}
for gamma in gamma_vals:
    gamma_sym.assign(gamma)
    print(f"\n=== Solving for γ̄ = {gamma} ===")
    results[gamma] = []
    for t in times:
        print(f"  t̄ = {t}")
        solver.solve()
        r_vals = Y.sub(0).vector().get_local()
        z_vals = Y.sub(1).vector().get_local()
        results[gamma].append((r_vals.copy(), z_vals.copy()))

        # ---- Protein diffusion update (forward Euler) ----
        s_nodes = mesh1d.coordinates().flatten()
        sigma_old = np.exp(sigma_tilde.vector().get_local())
        dsigma_num = np.gradient(sigma_old, s_nodes)
        M_arr = Y.sub(3).vector().get_local()
        dM_num = np.gradient(M_arr, s_nodes)
        flux_inner = alpha**2*sigma_old*dsigma_num - mu_phi*dM_num
        div_term = np.gradient(s_nodes*flux_inner, s_nodes)
        rhs = 2.0/(np.sqrt(s_nodes**2 + epsr**2)*sigma_old)*div_term
        sigma_new = np.maximum(sigma_old + dt*rhs, 1e-12)
        sigma_tilde.vector()[:] = np.log(sigma_new)
        sigma_tilde.vector()[-1] = np.log(beta/alpha)

# ==================== Axisymmetric surface plots====================
theta = np.linspace(0, 2*np.pi, 200)
for gamma, shapes in results.items():
    fig = plt.figure(figsize=(10, 4))
    for (r_vals, z_vals), t in zip(shapes, times):
        R, TH = np.meshgrid(r_vals, theta)
        X = R*np.cos(TH)
        Yp = R*np.sin(TH)
        Z = np.tile(z_vals, (len(theta), 1))
        ax = fig.add_subplot(1, len(times), list(times).index(t)+1, projection='3d')
        ax.plot_surface(X, Yp, Z, color='salmon', alpha=0.9, edgecolor='none')
        ax.set_title(f"γ̄={gamma}, t̄={t}")
        ax.set_axis_off()
    plt.suptitle(f"Membrane Budding Evolution (γ̄={gamma})")
    plt.show()

dokken · October 29, 2025, 8:16pm

This problem is so far away from a minimal example that it is really hard to give any guidance. As you can observe, your Newton solver diverges

✅ Mixed space created — DOFs: 11193

=== Solving for γ̄ = 0.0 ===
  t̄ = 0.0
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.812e+02 (tol = 1.000e-08) r (rel) = 1.000e+00 (tol = 1.000e-06)
  Newton iteration 1: r (abs) = -nan (tol = 1.000e-08) r (rel) = -nan (tol = 1.000e-06)

This can be because your initial guess is far away from the solution, or a load of other issues.
See Default absolute tolerance and relative tolerance - #4 by nate for tips and tricks.

I would especially pay attention to:

Does your Newton solver fail on the first step producing NaNs?
You likely have an initial guess leading to singularities in the Jacobian/residual. E.g: for a solution variable u
Zero initial guess with coefficients of the type 1/u or √u
Piecewise constant initial guess when invoking ∇u−1 or √εII(u)

In general I would advice you to build up your model step by step.

Topic		Replies	Views
Error: Newton solver did not converge Errors	1	803	June 12, 2022
Newton Solver did not converge mathematics	1	1321	August 13, 2020
Solver cannot converge General	3	148	June 11, 2024
RuntimeError: Newton solver did not converge because maximum number of iterations reached Errors	1	2333	April 20, 2022
What is the reason for nonconvergence of newton solver?	3	709	September 7, 2020