SNES solver giving up without converging

dolfinx
1

Hi all,
I’m trying to use the following non-linear solver to solve a variational plasticity problem with a mixed element state:

from dolfinx import la
from dolfinx.fem import form
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_matrix, create_vector, set_bc)
from ufl import TestFunction, TrialFunction, derivative
from petsc4py import PETSc

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):
        
        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
        apply_lifting(F, [self.a], bcs=[self.bc], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(F, self.bc, x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a, bcs=self.bc)
        J.assemble()

The solver settings are the following:

snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonrsls')
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(atol=1e-10,rtol=1.0e-9, max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-9)
snes.getKSP().getPC().setType("lu")

My method involves minimising a functional J_p, which encompasses free energy, lagrange multipliers and dissipation. The state consists of u, dp, N, lagrange1, lagrange2; [vector, scalar, tensor, scalar, scalar] respectively.

J_p =  elastic_strain_energy(u,full_plastic_strain)*dx 
J_p += dt*dissipation_potential(dp/dt)*dx 
J_p += dt*lagrange1*ufl.tr(N)*dx
J_p += dt*lagrange2*(ufl.inner(N,N)-2/3)*dx

Upon solving for a single timestep, all I get is:

time = 0.0000
  0 SNES Function norm 8.106060082376e+08

It seems to give up after a single iteration with an enormous residual, so I tried reducing the problem to linear elasticity by driving all the non-linearly dependent variables to zero.

J_p =  elastic_strain_energy(u,full_plastic_strain)*dx
J_p += ufl.inner(N,N)*dx
J_p += dp*dp*dx
J_p += lagrange1*lagrange1*dx
J_p += lagrange2*lagrange2*dx

Upon solving this, the solver does converge to the elastic solution, but still starts with the exact same huge residual.

time = 0.0000
  0 SNES Function norm 8.106060082376e+08 
  1 SNES Function norm 5.235158616409e-07

The problem may well be in my set up of the functional, J_p, but I’m more interested in why this solver doesn’t carry on iterating up to its max_it = 20 in the first case, it just seems to give up and spit out that large norm. Further, I would be interested to know why they both have exact same large norm initially. N.B. I have used this solver for linear viscosity with a very similar approach and it works fine, converging, using up all its iterations, sensible looking norms at each step.

A MWE in this context is quite difficult to make small, as the model is fairly involved, but I can provide one if needed.

Kind regards,

Magnus