Newton solver residual returns NaN after 2 iterations

Legacy DOLFIN
1

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

2

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.

3

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?

4

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

5

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