I am trying to solve a nonlinear Poisson equation. I hit some unexpected behaviour for a certain set of inputs, and I am unable to make sense of it. The NewtonSolver behaviour is puzzling to me, this is the solver output. Residuals are huge and then collapse to zero. Needless to say the solution is not ok.
2023-08-08 19:39:59.834 ( 238.182s) [main ] petsc.cpp:675 INFO| PETSc Krylov solver starting to solve system.
2023-08-08 19:39:59.841 ( 238.189s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = 2.87352e+115 (tol = 1e-10) r (rel) = 1.63805e+114(tol = 1e-06)
2023-08-08 19:39:59.843 ( 238.191s) [main ] petsc.cpp:675 INFO| PETSc Krylov solver starting to solve system.
2023-08-08 19:39:59.851 ( 238.199s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 4: r (abs) = 0 (tol = 1e-10) r (rel) = 0(tol = 1e-06)
2023-08-08 19:39:59.851 ( 238.199s) [main ] NewtonSolver.cpp:255 INFO| Newton solver finished in 4 iterations and 10008 linear solver iterations.
This is a MWE
import ufl
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh, fem, io, nls, log
from dolfinx.plot import create_vtk_mesh
def q(u):
return ufl.exp(u)
L, W, S = [1.66666667, 1.66666667, 1. ]
domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[L, W]]), [20,20], cell_type=mesh.CellType.quadrilateral)
x = ufl.SpatialCoordinate(domain)
u_ufl = x[1]*ufl.sin(2*3.14*x[0])*ufl.cos(S*3.14*x[1])
u_ufl_np = f"x[1]*np.sin(2*np.pi*x[0])*np.cos({S}*np.pi*x[1])"
f = - ufl.div(q(u_ufl)*ufl.grad(u_ufl))
V = fem.FunctionSpace(domain, ("CG", 2))
u_exact = lambda x: eval(u_ufl_np)
u_D = fem.Function(V)
u_D.interpolate(u_exact)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))
uh = fem.Function(V)
v = ufl.TestFunction(V)
F = q(uh)*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - f*v*ufl.dx
problem = fem.petsc.NonlinearProblem(F, uh, bcs=[bc])
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
#solver.rtol = 1e-1
solver.report = True
##########
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
##########
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")
For inputs as L, W, S = [1.66666667, 1.1, 1. ] works fine, for example.
Any hints please?
