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)

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?

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.