Interface problem from tutorial with SNES solver

[bbbyyfff]

Dear community,

I am trying to solve the interface problem in tutorial 3.3 with an SNES solver. I have already tried the linear ksp solver whereas the problem does not converge when I tried to use an SNES solver. The MWE for this problem is attached here:

import gmsh
import numpy as np
import ufl
import dolfinx
from dolfinx import fem, io, mesh

from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Vec, KSP, InsertMode, ScatterMode
from petsc4py import PETSc

import meshio
import dolfinx.io.gmshio as gmshio

from matplotlib import pyplot as plt

import multiphenicsx.fem as mfem

from typing import List, Tuple, Optional

# Geometrical parameters
L = 1.
sz = 1. / 100.

gmsh.initialize()
gmsh.model.add("interface")

p1 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, sz)
p2 = gmsh.model.geo.addPoint(L/2, 0.0, 0.0, sz)
p3 = gmsh.model.geo.addPoint(L, 0.0, 0.0, sz)
p4 = gmsh.model.geo.addPoint(L, L, 0.0, sz)
p5 = gmsh.model.geo.addPoint(L/2, L, 0.0, sz)
p6 = gmsh.model.geo.addPoint(0.0, L, 0.0, sz)

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, p6)
l6 = gmsh.model.geo.addLine(p6, p1)
l7 = gmsh.model.geo.addLine(p2, p5)

cloop_left = gmsh.model.geo.addCurveLoop([l1, l7, l5, l6])
cloop_right = gmsh.model.geo.addCurveLoop([l2, l3, l4, -l7])

s_left = gmsh.model.geo.addPlaneSurface([cloop_left])
s_right = gmsh.model.geo.addPlaneSurface([cloop_right])
gmsh.model.geo.synchronize()
# gmsh.model.mesh.setRecombine(2, s_left)
# gmsh.model.mesh.setRecombine(2, s_right)


gmsh.model.addPhysicalGroup(1, [l1, l2, l3, l4, l5, l6], 1)
gmsh.model.addPhysicalGroup(1, [l7], 2)
gmsh.model.addPhysicalGroup(2, [s_left], 1)
gmsh.model.addPhysicalGroup(2, [s_right], 2)

gmsh.model.mesh.generate(2)

partitioner = mesh.create_cell_partitioner(mesh.GhostMode.shared_facet)
msh, subdomains, boundaries = gmshio.model_to_mesh(gmsh.model, comm=MPI.COMM_WORLD, rank=0, gdim=2, partitioner=partitioner)
gmsh.finalize()

cells_Omega1 = subdomains.indices[subdomains.values == 1]  # Subdomain Omega 1
cells_Omega2 = subdomains.indices[subdomains.values == 2]  # Subdomain Omega 2

facets_partial_Omega = boundaries.indices[boundaries.values == 1] # Exterior boundary
facets_Gamma = boundaries.indices[boundaries.values == 2] # Interface

dx = ufl.Measure("dx")(subdomain_data = subdomains)
dS = ufl.Measure("dS")(subdomain_data = boundaries)

V = fem.FunctionSpace(msh, ("CG", 2))
V1 = V.clone()
V2 = V.clone()
M = V.clone()

dofs_V1_Omega1 = fem.locate_dofs_topological(V1, subdomains.dim, cells_Omega1)
dofs_V2_Omega2 = fem.locate_dofs_topological(V2, subdomains.dim, cells_Omega2)
dofs_M_Gamma = fem.locate_dofs_topological(M, boundaries.dim, facets_Gamma)

restriction_V1_Omega1 = mfem.DofMapRestriction(V1.dofmap, dofs_V1_Omega1)
restriction_V2_Omega2 = mfem.DofMapRestriction(V2.dofmap, dofs_V2_Omega2)
restriction_M_Gamma = mfem.DofMapRestriction(M.dofmap, dofs_M_Gamma)
restriction = [restriction_V1_Omega1, restriction_V2_Omega2, restriction_M_Gamma]

(u1, u2, l) = (fem.Function(V1), fem.Function(V2), fem.Function(M))
(du1, du2, dl) = (ufl.TrialFunction(V1), ufl.TrialFunction(V2), ufl.TrialFunction(M))
(v1, v2, m) = (ufl.TestFunction(V1), ufl.TestFunction(V2), ufl.TestFunction(M))


# Define residual energy
zero = dolfinx.fem.Constant(msh, ScalarType(0))

sub_expr_u1 = fem.Expression(zero, V1.element.interpolation_points())
u1.interpolate(sub_expr_u1)

sub_expr_u2 = fem.Expression(zero, V2.element.interpolation_points())
u2.interpolate(sub_expr_u2)

R = [ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx(1) + ufl.inner(l("-"), v1("-")) * dS(2) - ufl.inner(1, v1) * dx(1),
     ufl.inner(ufl.grad(u2), ufl.grad(v2)) * dx(2) - ufl.inner(l("+"), v2("+")) * dS(2) - ufl.inner(1, v2) * dx(2),
     ufl.inner(u1("-"), m("-")) * dS(2) - ufl.inner(u2("+"), m("+")) * dS(2) + ufl.inner(zero, m("-")) * dS(2)]

dR = [[ufl.derivative(R[0], u1), ufl.derivative(R[0], u2), ufl.derivative(R[0], l)],
      [ufl.derivative(R[1], u1), ufl.derivative(R[1], u2), ufl.derivative(R[1], l)],
      [ufl.derivative(R[2], u1), ufl.derivative(R[2], u2), ufl.derivative(R[2], l)]]

# Define boundary conditions
dofs_V1_partial_Omega = dolfinx.fem.locate_dofs_topological(V1, boundaries.dim, facets_partial_Omega)
dofs_V2_partial_Omega = dolfinx.fem.locate_dofs_topological(V2, boundaries.dim, facets_partial_Omega)
bc1 = dolfinx.fem.dirichletbc(zero, dofs_V1_partial_Omega, V1)
bc2 = dolfinx.fem.dirichletbc(zero, dofs_V2_partial_Omega, V2)
bcs = [bc1, bc2]

# Define SNES class
class SNES_Interface_Problem:
    '''
    Nonlinear problem class compatible with PETSC.SNES solver.
    '''

    def __init__(self, R: List[ufl.Form], dR: List[List[ufl.Form]], 
                  solutions: Tuple[fem.Function, fem.Function, fem.Function],
                  bcs: List[fem.DirichletBCMetaClass],
                  P: Optional[List[List[ufl.Form]]]=None):
        '''

        Parameters
        ----------
        R : List[ufl.Form]
            DESCRIPTION Residual.
        dR : List[List[ufl.Form]]
            DESCRIPTION Jacobian.
        solutions : Tuple[fem.Function, fem.Function, fem.Function, fem.Function]
            DESCRIPTION Solution functions.
        bcs : List[dolfinx.DirichletBC]
            DESCRIPTION Dirichlet boundary conditions.
        Returns
        -------
        None.

        '''

        self._R = fem.form(R)
        self._dR = fem.form(dR)
        self._obj_vec = mfem.petsc.create_vector_block(self._R, restriction)
        self._solutions = solutions
        self._bcs = bcs
        self._P = P

    def create_snes_solution(self):
        x = mfem.petsc.create_vector_block(self._R, restriction=restriction)
        with mfem.petsc.BlockVecSubVectorWrapper(x, [V1.dofmap, V2.dofmap, M.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.vector.localForm() as sub_solution_local:
                    x_wrapper_local[:] = sub_solution_local
        return x

    def update_solutions(self, x: PETSc.Vec):
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        with mfem.petsc.BlockVecSubVectorWrapper(x, [V1.dofmap, V2.dofmap, M.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.vector.localForm() as sub_solution_local:
                    x_wrapper_local[:] = sub_solution_local

    def obj(self, snes: PETSc.SNES, x: PETSc.Vec):
        self.R(snes, x, self._obj_vec)
        return self._obj_vec.norm()

    def R(self, snes: PETSc.SNES, x: PETSc.Vec, R_vec: PETSc.Vec):
        with R_vec.localForm() as R_vec_local:
            R_vec_local.set(0.0)
        mfem.petsc.assemble_vector_block(
            R_vec, self._R, self._dR, self._bcs, x0=x, scale=-1.0,
            restriction = restriction, restriction_x0=restriction)
    def dR(self, snes: PETSc.SNES, x: PETSc.Vec, dR_mat: PETSc.Mat, P_mat: PETSc.Mat):
        
        dR_mat.zeroEntries()
        mfem.petsc.assemble_matrix_block(
            dR_mat, self._dR, self._bcs, diagonal=1.0,
            restriction = (restriction, restriction))
        dR_mat.assemble()
        if self._P is not None:
            P_mat.zeroEntries()
            mfem.petsc.assemble_matrix_block(
                P_mat, self._P, self._bcs, diagonal=1.0,
                restriction=(restriction, restriction))
            P_mat.assemble()
            
problem = SNES_Interface_Problem(R, dR, (u1, u2, l), bcs)
R_vec = mfem.petsc.create_vector_block(problem._R, restriction=restriction)
dR_mat = mfem.petsc.create_matrix_block(problem._dR, restriction=(restriction, restriction))

# Solve
snes = PETSc.SNES().create(msh.comm)
snes.setType('newtonls')
snes.setTolerances(atol=1e-8, rtol=1e-8, max_it=50)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.getKSP().setTolerances(atol=1e-8, rtol=1.0e-8)

PETSc.Options()['snes_linesearch_type'] = 'basic'
PETSc.Options()['snes_monitor'] = None
PETSc.Options()['snes_linesearch_monitor'] = None
PETSc.Options()['ksp_converged_reason'] = None
PETSc.Options()['ksp_monitor_true_residual'] = None
PETSc.Options()['snes_converged_reason'] = None
snes.setFromOptions()

snes.setObjective(problem.obj)
snes.setFunction(problem.R, R_vec)
snes.setJacobian(problem.dR, J=dR_mat, P=None)

solution = problem.create_snes_solution()
snes.solve(None, solution)

problem.update_solutions(solution)
solution.destroy()
snes.destroy()

What I get is something like this:

0 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  1 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  2 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  3 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  4 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  5 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  6 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  7 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  8 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
  9 SNES Function norm 2.589553928459e-02 
    0 KSP none resid norm 2.589553928459e-02 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
    1 KSP none resid norm 3.152279663127e-15 true resid norm 3.152279663127e-15 ||r(i)||/||b|| 1.217306049696e-13
  Linear solve converged due to CONVERGED_ITS iterations 1
 10 SNES Function norm 2.589553928459e-02 
Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 10

Thanks for any help!

[francesco-ballarin]

I'll have a look at this next week.

[bbbyyfff]

Thank you very much! Sorry for this lengthy code, but the only difference from the tutorial is just switching the solver.

[francesco-ballarin]

There are two mistakes in your code:

1. you probably copied and pasted from create_snes_solution to update_solutions, but watch out that the final line in update_solutions should read sub_solution_local[:] = x_wrapper_local instead of the other way around
2. you are not calling self.update_solutions(x) at the very beginning of the method R, therefore the solution computed in the Newton step is never getting copied into u1, u2, l.

With these changes, I get convergence in one nonlinear iteration.

[bbbyyfff]

Thank you so much! That solved my problem!