# Return infinite norms in the custom newton solver

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()

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)

line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])

gmsh.model.geo.synchronize()
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.scale(-1)

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]
``````

These formulations don’t really make sense. Particularly with no prescription of boundary conditions. It would be helpful to see the original system that you’re attempting to discretise.

@Yankang_Liu I interpreted your message to flag the post as a request to close it. If my interpretation is wrong, please get in touch with me and I’ll reopen it.