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

221
views
1
6 months ago by

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

6
6 months ago by
"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
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)​