Hello I was trying to debug my code and found the following issue…
I am solving a nonlinear system of equations in a 2d mesh with second order partial derivative (so order for basis function no less than two is necessary) and I used custom Newton solver to solve the algebraic system. But the solver continuously returned me inf which I had no idea. So I simplified the expression as much as possible.
It seems when there is something like (ux.dx(0))**2 and ux.dx(0).dx(0), the iteration returns me inf. I attached the code and the expressions.
Note these are not the expressions I am working on. They are just to reproduce the inf. When playing with these expressions, I only wanted to see for which I could have a norm that was not infinite.
import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
a = 3.
lcar = 2.
gmsh.initialize()
gmsh.model.add("mesh")
p0 = gmsh.model.geo.addPoint(a, 0, 0, lcar)
p1 = gmsh.model.geo.addPoint(a/2, np.sqrt(3.)*a/2, 0, lcar)
p2 = gmsh.model.geo.addPoint(-a/2, np.sqrt(3)*a/2, 0, lcar)
p3 = gmsh.model.geo.addPoint(-a, 0, 0, lcar)
p4 = gmsh.model.geo.addPoint(-a/2, -np.sqrt(3)*a/2, 0, lcar)
p5 = gmsh.model.geo.addPoint(a/2, -np.sqrt(3)*a/2, 0, lcar)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)
line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
domain = gmsh.model.geo.addPlaneSurface([line_loop])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0], 1)
gmsh.model.addPhysicalGroup(1, [l1], 2)
gmsh.model.addPhysicalGroup(1, [l2], 3)
gmsh.model.addPhysicalGroup(1, [l3], 4)
gmsh.model.addPhysicalGroup(1, [l4], 5)
gmsh.model.addPhysicalGroup(1, [l5], 6)
gmsh.model.addPhysicalGroup(2, [domain], 0)
gmsh.model.mesh.generate()
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2
)
gmsh.finalize()
boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 2]
boundaries_3 = boundaries.indices[boundaries.values == 3]
boundaries_4 = boundaries.indices[boundaries.values == 4]
boundaries_5 = boundaries.indices[boundaries.values == 5]
boundaries_6 = boundaries.indices[boundaries.values == 6]
Y = 1
nu_p = 0.2
tau = 0
V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))
vxvy = ufl.TestFunction(V)
(vx, vy) = vxvy[0], vxvy[1]
duxuy = ufl.TrialFunction(V)
uxuy = dolfinx.fem.Function(V)
(ux, uy) = uxuy[0], uxuy[1]
F = (ux)**1 * vx * ufl.dx + (uy)**1 * vy * ufl.dx
F = (ux)**2 * vx * ufl.dx + (uy)**1 * vy * ufl.dx
F = (ux.dx(0))**1 * vx * ufl.dx + (uy.dx(1))**2 * vy * ufl.dx # inf
F = (ux.dx(0))**1 * vx * ufl.dx + (uy.dx(1))**1 * vy * ufl.dx
F = (ux.dx(0))**2 * vx * ufl.dx + (uy.dx(1))**1 * vy * ufl.dx # inf
F = (ux.dx(1))**1 * vx * ufl.dx + (uy.dx(0))**1 * vy * ufl.dx
F = (ux.dx(1))**2 * vx * ufl.dx + (uy.dx(0))**1 * vy * ufl.dx # inf
F = (ux.dx(1).dx(0))**1 * vx * ufl.dx + (uy.dx(0))**1 * vy * ufl.dx # inf
F = (ux.dx(0).dx(0))**1 * vx * ufl.dx + (uy.dx(0))**1 * vy * ufl.dx # inf
F = (ux.dx(0))**1 * vx * ufl.dx + (uy)**1 * vy * ufl.dx
F = (ux.dx(0).dx(0))**1 * vx * ufl.dx + (uy)**1 * vy * ufl.dx # inf
F = (ux.dx(0))**2 * vx * ufl.dx + (uy)**1 * vy * ufl.dx # inf
uxuy.x.array[:] = 2.
residual = dolfinx.fem.form(F)
J = ufl.derivative(F, uxuy)
jacobian = dolfinx.fem.form(J)
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = petsc4py.PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
dx = dolfinx.fem.Function(V)
ii = 0
coords = V.tabulate_dof_coordinates()[:, 0]
sort_order = np.argsort(coords)
max_iterations = 25
solutions = np.zeros((max_iterations + 1, len(coords)))
solutions[0] = uxuy.x.array[sort_order]
ii = 0
while ii < max_iterations:
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, jacobian)
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD_VALUES, mode=petsc4py.PETSc.ScatterMode.REVERSE)
L.scale(-1)
L.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT_VALUES, mode=petsc4py.PETSc.ScatterMode.FORWARD)
solver.solve(L, dx.vector)
dx.x.scatter_forward()
uxuy.x.array[:] += dx.x.array
ii += 1
correction_norm = dx.vector.norm(0)
print(f"Iteration {ii}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
solutions[ii, :] = uxuy.x.array[sort_order]