### Newton Solver not converging with nonlinear system

187

views

2

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 \
+ (1 + Γb) * grad(rho1)[0] * grad(v1)[0] * dr \
+ grad(rho2)[0] * grad(v2)[0] * dr \
+ Γb * grad(rho3)[0] * grad(v3)[0] * dr
for i in range(1, num_steps + 1):
solve(F == 0, rho)
rho_n.assign(rho)
```

Community: FEniCS Project

### 1 Answer

2

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:

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

Please login to add an answer/comment or follow this question.

Edit: My mistake, it's time dependent.