Difficult to obtain a convergence result using SNES

General
1

I am trying to reproduce an example from Mathematica:
u_1'(x)-u_2(x)=0
u_2'(x)+u_1(x)-5(1-u_1^2(x))u_2(x)=0
with u_1(0)=2 and u_2(0)=0.

The solution should be like this (confirmed from Mathamtica):
u1
u2.

This is a nonlinear problem. So when solving the algebraic system, proper initial guess should be assigned. The function is shaking violently, so I choose a random array as the initial guess. To obtain this solution, I had to run the code again and again until I saw a convergence in the iteration. The reason might still lie on the shape of the functions. Also, if I shrank the interval, say [0, 4], it became easier to converge and harder when prolong the interval say [0, 15]. This is understandable from the same view: the shorter interval corresponds to smoother functions while the longer interval may introduce another peak. In the example I choose [0, 8]. I was wondering if there is a better way to solve the algebraic system. I was trying to adjust

    snes.getKSP().setType("preonly")
    snes.getKSP().getPC().setType("lu")
    snes.getKSP().getPC().setFactorSolverType("mumps")

but didn’t know how to do exactly… And I also wonder if there is another way to improve the code.

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

n = 300
(x_min, x_max) = (0, 8)

mesh = dolfinx.mesh.create_interval(mpi4py.MPI.COMM_WORLD, n, [x_min, x_max])
tdim = mesh.topology.dim
fdim = tdim - 1

facet_left = dolfinx.mesh.locate_entities_boundary(
    mesh, fdim, lambda x: np.isclose(x[0], x_min)
)

boundaries = [(1, lambda x: np.isclose(x[0], x_min)),
              (2, lambda x: np.isclose(x[0], x_max))]
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)

facet_tag = dolfinx.mesh.meshtags(mesh, fdim,
                                  facet_indices, facet_markers)
mesh.topology.create_connectivity(fdim, tdim)

ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)

U1_element = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
U2_element = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 2)

def run_monolithic():
    U_element = ufl.MixedElement(U1_element, U2_element)
    U = dolfinx.fem.functionspace(mesh, U_element)

    du1u2 = ufl.TrialFunction(U)
    u1u2 = dolfinx.fem.Function(U)
    (u1, u2) = ufl.split(u1u2)
    v1v2 = ufl.TestFunction(U)
    (v1, v2) = ufl.split(v1v2)
    u1u2.x.array[:] = 2 * np.random.rand(1,1)[0][0]

    F = u1 * v1 * ds(2) - 2 * v1 * ds(1) - u1 * v1.dx(0) * ufl.dx - u2 * v1 * ufl.dx\
        + (u2.dx(0) + u1 - 5 * (1 - u1**2) * u2) * v2 * ufl.dx
    J = ufl.derivative(F, u1u2, du1u2)

    u2_left = dolfinx.fem.Function(U.sub(1).collapse()[0])
    u2_left.x.array[:] = 0
    bdofs_U2 = dolfinx.fem.locate_dofs_topological(
        (U.sub(1), U.sub(1).collapse()[0]), fdim, facet_left)
    left_bc_u2 = dolfinx.fem.dirichletbc(u2_left, bdofs_U2, U.sub(1))
    bc = [left_bc_u2]

    class coupled2_1(object):
        def __init__(
                self, F, J, solution, bcs, P=None
        ):
            self._F = dolfinx.fem.form(F)
            self._J = dolfinx.fem.form(J)
            self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
            self._solution = solution
            self._bcs = bcs
            self._P = P

        def create_snes_solution(self):
            x = self._solution.vector.copy()
            with x.localForm() as _x, self._solution.vector.localForm() as _solution:
                _x[:] = _solution
            return x
        
        def update_solution(self, x):
            x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
            with x.localForm() as _x, self._solution.vector.localForm() as _solution:
                _solution[:] = _x
        
        def obj(
                self, snes, x
        ):
            self.F(snes, x, self._obj_vec)
            return self._obj_vec.norm()
        
        def F(
                self, snes, x, F_vec
        ):
            self.update_solution(x)
            with F_vec.localForm() as F_vec_local:
                F_vec_local.set(0.0)
            dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
            dolfinx.fem.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)
            F_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
            dolfinx.fem.set_bc(F_vec, self._bcs, x, -1.0)

        def J(
                self, snes, x, J_mat, P_mat
        ):
            J_mat.zeroEntries()
            dolfinx.fem.petsc.assemble_matrix(
                J_mat, self._J, self._bcs, diagonal=1.0
            )
            J_mat.assemble()
            if self._P is not None:
                P_mat.zeroEntries()
                dolfinx.fem.petsc.assemble_matrix(
                    P_mat, self._P, self._bcs, diagonal=1.0
                )
                P_mat.assemble()

    problem = coupled2_1(F, J, u1u2, bc)
    F_vec = dolfinx.fem.petsc.create_vector(problem._F)
    J_mat = dolfinx.fem.petsc.create_matrix(problem._J)

    snes = petsc4py.PETSc.SNES().create(mesh.comm)
    snes.setTolerances(max_it=200)
    snes.getKSP().setType("preonly")
    snes.getKSP().getPC().setType("lu")
    snes.getKSP().getPC().setFactorSolverType("mumps")
    snes.setObjective(problem.obj)
    snes.setFunction(problem.F, F_vec)
    snes.setJacobian(problem.J, J=J_mat, P=None)
    snes.setMonitor(lambda _, it, residual: print(it, residual))
    u1u2_copy = problem.create_snes_solution()
    snes.solve(None, u1u2_copy)
    problem.update_solution(u1u2_copy)
    u1u2_copy.destroy()
    snes.destroy()
    return u1u2

u1u2_m = run_monolithic()
(u1_m, u2_m) = (u1u2_m.sub(0).collapse(), u1u2_m.sub(1).collapse())

viskex.dolfinx.plot_scalar_field(u1_m, "u1_m")

viskex.dolfinx.plot_scalar_field(u2_m, "u2_m")
2

From a glance I find this… unconventional. Is there a reason why you’re doing this? In my experience nonlinear problems are usually best initialized with initial guesses that are either constant everywhere (conforming to the boundary conditions) or the solution of the previous time step (in case of dynamic problems). The convergence behavior of Newton’s method is sensitive to the initial guess, so initializing with a random array (with high-frequency oscillations) seems like a poor choice to me.

Which seems to be consistent with your next sentence:

Which sounds reasonable, if my interpretation (you keep generating random initial guesses until you generate one that leads to solution convergence) is correct.

Lastly: Most (if not all) of your code seems to be very similar to the built-in Dolfinx NewtonSolver. You could use that class to debug your problem. Its functionalities have been validated, so by comparing its output you can verify whether your nonlinear solver implementation is correct.

3

I agree with @spcvanschie that the issue is probably the initial guess.

On the code itself: @Yankang_Liu has copied the second tutorial in multiphenicsx, in which I show how to use SNES. I can guarantee that that code has been validated :wink:

Tweaking the solver of the linear system typically will not improve the convergence of the nonlinear one. You should be changing options of the nonlinear solver, for instance SNESLineSearchType — PETSc 3.20.0 documentation is a notable one, and can be set as in [petsc-users] petsc4py - Change default line search for SNES Newton Line Search

Finally, a suggestion for future reference for @Yankang_Liu. In this forum is typically better to a post a minimal example. For instance, you could have stopped your code at the line u1u2_m = run_monolithic(), since there is no need to ask users to also read the other case, as I guess the first one already shows your issue.

2 Likes
4

Thank you very much. I was just doing this exercise to solve coupled nonlinear equations. The equations I am going to solve are very long although there are only two equations in 1d. So I want to get more adept on handling this problem.

I used finite element to solve a simpler version which only contained one of the two equations (the second unknown was assigned to a given function) and found it was initial guess sensitive. I could only set the guess as x[0]**1.1 (not even x[0]) to get the correct answer… And I have to solve the equation again and again for a bunch of different parameters. I worry there may not be a convenient way to write a universal code for the process and happened to see this example on Mathematica and grab.

5

Thank you very much for your reply! I have deleted the redundant part and look into the options!

6

You could look into using Picard Iteration to deal with the nonlinearity. As far as I’m aware that’s a more stable (but slower) linearization approach than Newton’s method.

Have you tried an initial guess of all zeros?

7

Thank you very much for pointing out a different iteration! But for

Have you tried an initial guess of all zeros?

It returned cannot converge for solver.mat_it = 400. Sorry I prefer not to attach the code for you to play with since my current project is directly related to this one.