### Why am I getting r (rel) = -nan with the Newton solver

221

views

1

I am trying to solve a PDE but I am getting an error,

```
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
```

This is created with the following code.

```
from fenics import *
mesh = IntervalMesh(20, 0, 1)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
rho = Function(V)
rho_n1 = Function(V)
rho_n2 = Function(V)
dt = Constant(0.1)
dr = dx
F = ((rho - 2 * rho_n1 + rho_n2) / dt**2) * v * dr \
+ ((rho - rho_n1) / dt) * nabla_grad(v)[0] * dr
solve(F == 0, rho)
```

Any tips on how to fix this issue would be appreciated.

Community: FEniCS Project

You don't have any trial function

written
6 months ago by
hirshikesh

### 1 Answer

6

"You don't have any trial function"

That shouldn't be a problem. The derivative is calculated automatically. So what's happening in CameronDevine's code is the same as would happen if we wrote

Maybe think about your boundary conditions. At the moment you implicitly prescribe zero Neumann boundary conditions. If I add a Dirichlet BC it converges. Although I can't say anything about how physically meaningful this is:

That shouldn't be a problem. The derivative is calculated automatically. So what's happening in CameronDevine's code is the same as would happen if we wrote

```
drho = TrialFunction(V)
dF = derivative(F,rho,drho)
solve(F == 0, rho, J=dF)
```

The real problem I think is, that rho=0 everywhere is the solution to the problem which is why the initial error is already zero. If the error is zero, then the relative error, i.e. the error divided by the absolute error is nan.Maybe think about your boundary conditions. At the moment you implicitly prescribe zero Neumann boundary conditions. If I add a Dirichlet BC it converges. Although I can't say anything about how physically meaningful this is:

```
from fenics import *
mesh = IntervalMesh(20, 0, 1)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
rho = Function(V)
rho_n1 = Function(V)
rho_n2 = Function(V)
dt = Constant(0.1)
dr = dx
F = ((rho - 2 * rho_n1 + rho_n2) / dt**2) * v * dr \
+ ((rho - rho_n1) / dt) * nabla_grad(v)[0] * dr
def right(x, on_boundary):
return near(x[0],1) and on_boundary
dbc = DirichletBC(V, Constant(1.0), right)
drho = TrialFunction(V)
dF = derivative(F,rho,drho)
solve(F == 0, rho, bcs = dbc, J=dF)
```

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