Set bounds in a NonLinearProblem

dolfinx
6

Sorry to continue a thread already marked as solved, but I am not sure I have figured out the right way of doing this and need some more help.

Here is a MWE adapted from the unit test for nonlinear pde snes.

# Copyright (C) 2018 Garth N. Wells
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
"""Unit tests for Newton solver assembly"""

import numpy as np

import ufl
from dolfinx import cpp as _cpp
from dolfinx import la
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
                         locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_matrix, create_vector, set_bc)
from dolfinx.mesh import create_unit_square
from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner

from mpi4py import MPI
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()

def test_nonlinear_pde_snes(bound):
    """Test Newton solver for a simple nonlinear PDE"""
    # Create mesh and function space
    mesh = create_unit_square(MPI.COMM_WORLD, 12, 15)
    V = FunctionSpace(mesh, ("Lagrange", 1))
    u = Function(V)
    u.x.array[:] = 15
    u.name = 'solution'
    from dolfinx import io
    with io.XDMFFile(MPI.COMM_WORLD,'u.xdmf','w') as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(u,0)
    v = TestFunction(V)
    F = inner(5.0, v) * dx - ufl.sqrt(u * u) * inner(
        grad(u), grad(v)) * dx - inner(u, v) * dx
    
    u_bc = Function(V)
    u_bc.x.array[:] = 1.0
    bc = dirichletbc(u_bc, locate_dofs_geometrical(V, lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                              np.isclose(x[0], 1.0))))
    
    # Create nonlinear problem
    problem = NonlinearPDE_SNESProblem(F, u, bc)
    
    b = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
    J = create_matrix(problem.a)
    
    # Create Newton solver and solve
    snes = PETSc.SNES().create()
    opts = PETSc.Options()
    opts['snes_monitor'] = None
    snes.setFromOptions()
    snes.setType('vinewtonssls')
    snes.setFunction(problem.F, b)
    snes.setJacobian(problem.J, J)
    snes.setTolerances(atol=1e-10,rtol=1.0e-9, max_it=10)
    snes.getKSP().setType("preonly")
    snes.getKSP().setTolerances(rtol=1.0e-9)
    snes.getKSP().getPC().setType("lu")
    if bound:
        u_min, u_max = Function(V), Function(V)
        u_min.x.array[:] = 1.0
        u_max.x.array[:] = 1.5
        snes.setVariableBounds(u_min.vector, u_max.vector)
    snes.solve(None, u.vector)
    print('Minimum in u array:',np.min(u.x.array),'\nMaximum in u array:', np.max(u.x.array))

    from dolfinx import io
    with io.XDMFFile(MPI.COMM_WORLD,'u.xdmf','a') as xdmf:
        xdmf.write_function(u,1)

print('No bounds')
test_nonlinear_pde_snes(bound=False)
print('\nBounds set between 1.0 and 1.5')
test_nonlinear_pde_snes(bound=True)

When I run this code, I get the following results,

No bounds
0 SNES Function norm 9.083961017374e+02
1 SNES Function norm 2.603033612549e-01
2 SNES Function norm 6.870718004914e-02
3 SNES Function norm 1.678228243832e-03
4 SNES Function norm 8.126048345290e-07
Minimum in u array: 1.0
Maximum in u array: 1.3851364823572738

Bounds set between 1.0 and 1.5
0 SNES Function norm 2.202572332184e+00
1 SNES Function norm 4.139002904298e-01
Minimum in u array: -0.012085010497842469
Maximum in u array: 0.03760750671222456

Considering that the final solution without bounds lies between 1.0 and 1.5, I expected the solver to converge to that solution without any issues even with those bounds set. But this doesn’t seem to be the case. Is there something wrong with the way I am setting the bounds on the problem?

Also, the two solvers that PETSc provides for problems with bounds, vinewtonssls and vinewtonrsls give different results.

Could someone please share their expertise in this matter?

Thank you
Krishna Kenja