Dear community,
I am currently working on Cahn-Hilliard simulations, using the DOLFINx Python demo as the base code. My primary modification involves changing the equation for f.
During the simulation, I encountered an issue with NaN values in the convergence rate at certain steps:
[dolfinx] [info] Newton iteration 1: r (abs) = nan (tol = 1e-10), r (rel) = nan (tol = 1e-06)
Initially, I suspected the error might be caused by negative arguments in the logarithmic function within f. To address this, I attempted the solution described in Square Root and Natural Logarithm not Working. However, even after implementing this fix, the code still failed.
The error consistently occurs at t=210. For reference, I used the latest version of the library available on Conda (v0.9.0). Below, I have attached the code for further investigation.Thank in advance for any help.
Code example:
from petsc4py import PETSc
import dolfinx
from mpi4py import MPI
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log
from dolfinx.fem import Function, functionspace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.mesh import CellType, create_rectangle
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
import ufl.algebra
log.set_output_file("log.txt")
# system parameters
temp = 483
chi = 0.008756
N1 = 8500
N2 = 4600
phi1_init = 0.2
Rg = 1.13 * 10e-8
grad_coef = (chi/3) ** 0.5
M = 1
# Master parameters:
dt = 10
final_time = 3000
theta = 0.5
mesh_resolution = 80
domain_size = 20
EPS = 1e-15
dtype = np.float64 # force typing for mesh and finite elements
msh = create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (domain_size, domain_size)], [mesh_resolution, mesh_resolution], CellType.triangle, dtype=np.real(dtype(0)).dtype)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=np.real(dtype(0)).dtype)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME)
u0 = Function(ME)
phi1, mu = ufl.split(u)
phi1_0, mu_0 = ufl.split(u0)
u.x.array[:] = 0.0
rng = np.random.default_rng(13)
u.sub(0).interpolate(lambda x: phi1_init + 0.02 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()
phi1 = ufl.variable(phi1)
f = phi1 / N1 * ufl.ln(phi1 + EPS) + (1 - phi1) / N2 * ufl.ln(1 - phi1 + EPS) + chi * phi1 * (1 - phi1)
dfdc = ufl.diff(f, phi1)
mu_mid = (1.0 - theta) * mu_0 + theta * mu
# Weak statement of the equations
F0 = inner(phi1, q) * dx - inner(phi1_0, q) * dx + dt * M * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - grad_coef ** 2 * inner(grad(phi1), grad(v)) * dx
F = F0 + F1
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-6
solver.max_it = 50
solver.error_on_nonconvergence = True
solver.relaxation_parameter = 1
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()
# Step in time
t = 0
V0, dofs = ME.sub(0).collapse()
phi1 = u.sub(0)
u0.x.array[:] = u.x.array
log.set_log_level(log.LogLevel.INFO)
while t < final_time:
t += dt
r = solver.solve(u)
print(f"t' {t}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
