# Newton solver residual returns NaN after 2 iterations

Hi, I’m new to fenics and FEM, I tried to adapt the Cahn-Hilliard Equation demo downloaded here https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html for my own application. I literally only changed the free energy functional as

``````alpha_fchem  = 1/(8.314*300)                   # coeff to convert g to phase field dimensionless g
G0 = -336668.3750 # G0 is the pure substance gibbs free energy
Omega0 = 128205.6641
Omega1 = 61896.4375
Omega2 = -290954.8125
Omega3 = -164259.3750
f = alpha_fchem*(c*G0 + (1-c)*0.0 + 8.314*300.0*(c*ln(c)+(1-c)*ln(1-c))+
c*(1-c)*(Omega0
+ Omega1*(1-2*c)
+ Omega2*(1-2*c)**2
+ Omega3*(1-2*c)**3 ))
``````

instead of the original one

``````f    = 100*c**2*(1-c)**2
``````

and I got nan residual after 2 iterations:
Newton iteration 0: r (abs) = 1.726e+04 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 7.432e+03 (tol = 1.000e-10) r (rel) = 4.305e-01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)

I also tried to set a non-zero initial guess

``````u0.vector().set_local(rand(u0.vector().size()))
u0.vector().apply("")
``````

as mentioned here https://fenicsproject.org/qa/8703/nonlinear-dependent-iteration-shows-remaining-iteration/, but I still got nan after 2 iterations. I can confirm that this free energy functional I plugged in is a valid functional, it’s fairly common in my research area. Anyone can give me a hint on how I can resolve this error? Thanks a lot

The logarithms are likely causing problems. You could add a guard to their arguments make sure their functional values don’t blow up, or maybe try a smaller relaxation parameter. Consider also: Default absolute tolerance and relative tolerance - #4 by nate.

Hi Nate, thanks for your quick response! I don’t think the log function is causing problems, cause I tried the following functional form (basically getting rid of those Omega*(1-2*c)**I terms)

``````f = alpha_fchem*(c*G0 + (1-c)*0.0 + 8.314*300.0*(c*ln(c)+(1-c)*ln(1-c)))
``````

and the Newton solver doesn’t have any problem.

However, if I use

``````f = alpha_fchem*(c*G0 + (1-c)*0.0 + 8.314*300.0*(c*ln(c)+(1-c)*ln(1-c))
+ c*(1-c)*(Omega0)
``````

then I’ll have the nan back again:
Newton iteration 0: r (abs) = 1.434e+04 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.822e-03 (tol = 1.000e-10) r (rel) = 2.666e-07 (tol = 1.000e-06)
Newton solver finished in 2 iterations and 2 linear solver iterations.
Newton iteration 0: r (abs) = 3.850e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.769e-02 (tol = 1.000e-10) r (rel) = 7.193e-04 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.294e-08 (tol = 1.000e-10) r (rel) = 3.360e-10 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Newton iteration 0: r (abs) = 1.394e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 4.629e-01 (tol = 1.000e-10) r (rel) = 3.320e-03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.290e-05 (tol = 1.000e-10) r (rel) = 9.256e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Newton iteration 0: r (abs) = 5.264e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 3: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)
Newton iteration 4: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-06)

and these Omega*(1-2*c)**I terms have physical meanings, so they cannot be just ignored when carrying out simulations. Is there anything else I can do to resolve the error?

also how to add a guard? Can you point me to a tutorial or other relevant pages, thanks

In the context of `ln(c)` something like

``````tol = dolfinx.fem.Constant(mesh, 1e-6)
c_guard = ufl.conditional(ufl.lt(c, tol), tol, c)
ln(c_guard)
``````

or

``````tol = dolfinx.fem.Constant(mesh, 1e-6)
c_guard = ufl.max_value(c, tol)
ln(c_guard)
``````

Essentially you make a “guard” to protect the logarithm call from unexpected values in the current solution vector approximation in your (e.g. Newton) iterative solve. `c` could have negative values from oscillations in your numerical approximation, for example.

1 Like