Non linear solver

variational formulation
1

Hi all,
I’d like to ask you an advice regarding the choice of the non linear solver because the two solvers I choose does not converge as I refine my mesh. The problem I’m trying to solve, with a non linear variational approach, is a hyperelastic model for torsion deformation.
I’ve already tested the “snes” solver and the “newton_solver” in the following way:

    FF, dF = self.FirstAndSecondVariation()
    problem = NonlinearVariationalProblem(FF, self.up, self.bcs, dF)
    solver = NonlinearVariationalSolver(problem)
    solver.parameters.update(params)
    # at the beginnig of my code
    params = {'nonlinear_solver': 'snes',
           'snes_solver':
            {
                'linear_solver'           : 'mumps',
                'absolute_tolerance'      : 1e-10,
                'relative_tolerance'      : 1e-10,
                'maximum_iterations'      : 20,
             }
         }

or

FF, dF = self.FirstAndSecondVariation()
problem = NonlinearVariationalProblem(FF, self.up, self.bcs, dF)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-10
prm['newton_solver']['relative_tolerance'] = 1E-10
prm['newton_solver']['maximum_iterations'] = 20
prm['newton_solver']['relaxation_parameter'] = 1.0

Thank you for your help!

1 Like
2

SNES should be pretty robust as far as hyperelasticity goes. How ill-conditioned is your problem?

Another approach could be to use a suitable Krylov solver (cg/bicgstab) with a multigrid preconditioner (petsc_amg) instead of a direct solver (mumps.superlu_dist), since for larger problems the latter may run out of memory. You could do something similar to https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/python/demo/undocumented/elasticity/demo_elasticity.py#lines-35 on your assembled Jacobian (dF)…

1 Like
3

I would implement the nonlinear solver loop by myself, in dolfin-x style something like

# create solver
ksp = PETSc.KSP().create(comm)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")

while it < maxiter:
    
    # assemble rhs vector
    r = assemble_vector(varform)
    apply_lifting(r, [jac], [dbcs], x0=[u.vector], scale=-1.0)
    r.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    set_bc(r, dbcs, x0=u.vector, scale=-1.0)
    
    # assemble system matrix
    K = assemble_matrix(jac, dbcs)
    K.assemble()

    # solve linear system
    ksp.setOperators(K)
    ksp.solve(-r, del_u.vector)
    
    # get residual and increment norm
    struct_res_norm = r.norm()
    struct_inc_norm = del_u.vector.norm()
    
    # update solution
    u.vector.axpy(1.0, del_u.vector)
    u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    
    if comm.rank == 0:
        print("%i     %.4e         %.4e" % (it, struct_res_norm, struct_inc_norm))
        sys.stdout.flush()
    
    # increase iteration index
    it += 1

    # check if converged
    if struct_res_norm <= tolres and struct_inc_norm <= tolinc:
        break

else:

    if comm.rank == 0:
        print("Newton did not converge after %i iterations!" % (it))
        sys.stdout.flush()

and then would do a “pseudo-transient continuation”-like regularization (PTC) on the system matrix, like K += alpha * I, in PETSc done with

K.shift(alpha)

where you initialize alpha by some scalar value and update it after each iteration with

alpha *= struct_res_norm/struct_res_norm_last

Like this, the closer you get to the solution, the more the method behaves like a standard Newton.
I also use an adaptive scheme where alpha is only set to non-zero if there is a NAN in the res norm (can be checked with numpy’s isnan(struct_res_norm)).
Especially for hyperelasticity, this literally always helps me. :slight_smile:

Best,
Marc

4 Likes