### Newton Solver not converging with nonlinear system

187
views
2
4 months ago by

I am attempting to solve a PDE, but I am getting an error that the solution has not converged in the maximum number of iterations. What can I do to resolve this issue? A minimal working example is below.

from fenics import *

# Defining an invese hyperbolic tangent because fenics doesn't have one
def atanh(x): # Only good for abs(x) < 1!
return ln((1 + x) / (1 - x)) / 2

T = 0.1
num_steps = 50

Γb = 35819.838392953796
cb = 0.7292734109596469
δb = 0.0005651728598954567
γb = 658.2106848362291

dt = T / num_steps

mesh = IntervalMesh(50, -7.512, 7.512)

V = VectorFunctionSpace(mesh, 'CG', 1, dim = 3)

v1, v2, v3 = TestFunctions(V)

rho = Function(V)
rho_n = Function(V)

rho1, rho2, rho3 = split(rho)
rho1_n, rho2_n, rho3_n = split(rho_n)

r = Expression('x[0]', degree = 1)

rho_ic = Expression(
['0', '-0.1', '0.1+0.1*((x[0]<1)&&(x[0]>-1))'],
degree = 1)

rho20 = rho_ic[1]
rho30 = rho_ic[2]

rho.assign(rho_ic)
rho_n.assign(rho_ic)

_dt = Constant(dt)
dr = dx

F = ((rho1 - rho1_n) * v1 / _dt
+ (rho2 - rho2_n) * v2 / _dt
+ (rho3 - rho3_n) * v3 / _dt) * dr \
-((cb**2 / (1 + δb)) * ((1 - rho2**2) + Γb * δb * γb**2 * (1 - rho3**2)) * v1 * atanh(rho1)
- (cb / (1 + δb) * (grad(rho2)[0] - Γb * δb * γb * grad(rho3)[0]) * v1)) * dr \
+ ((-cb * ((1 - rho2**2) / (1 - rho1**2)) * grad(rho1)[0]
+ 2 * cb * rho2 * atanh(rho1) * grad(rho2)[0]) * v2) * dr \
+ ((-cb * ((1 - rho3**2) / (1 - rho1**2)) * grad(rho1)[0]
+ 2 * cb * rho3 * atanh(rho1) * grad(rho3)[0]) * v3) * dr \

for i in range(1, num_steps + 1):
solve(F == 0, rho)
rho_n.assign(rho)
Community: FEniCS Project
Boundary conditions?

Edit: My mistake, it's time dependent.
written 4 months ago by Nate

2
4 months ago by
Looks like your system is converging to a zero steady state. The default relative and absolute tolerances appear too small for your system.

Try using the following to solve your problem, where I have relaxed the absolute tolerance of the Newton system residual:

J = derivative(F, rho)

class NLP(NonlinearProblem):
def F(self, b, x):
assemble(F, tensor=b)
def J(self, A, x):
assemble(J, tensor=A)

solver = PETScSNESSolver()
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")
PETScOptions.set("snes_atol", 1e-8)
solver.set_from_options()
for i in range(1, num_steps + 1):
solver.solve(NLP(), rho.vector())
rho_n.assign(rho)​
In the future I may want to impose a Dirichlet boundary condition. Is there any was to do this with this solver? With this additional information is there a different solver you would recommend?
written 4 months ago by Cameron Devine
You don't need to change the solver. Just assemble the problem correctly with bcs inside NLP.
written 4 months ago by Nate