Null space is not set on underlying KSP PC object when using PETScSNESSolver

[Jack Hale]

I want to use a smoothed aggregation multigrid method to precondition the solution of the iterations of my SNES solver. To get optimal convergence it is necessary to pass a null space to the preconditioner.

Garth has written an example of using the PETScPreconditioner class in conjunction with the PETScKrylovSolver class here:

https://bitbucket.org/fenics-project/dolfin/src/b67e6aa0645e4fa17eda05df21b2258566685155/demo/undocumented/elasticity/python/demo_elasticity.py?at=master

My initial (incorrect) thought on how to do something similar with PETScSNESSolver was to pull out the underlying ksp object, make a DOLFIN PETScKrylovSolver from it and then use the 'set' method of PETScPreconditioner to set the PC and the nullspace on the ksp:

problem = GeneralProblem(a, L, bc)

nl_solver = PETScSNESSolver()
snes = nl_solver.snes()
ksp = nl_solver.snes().ksp

solver = PETScKrylovSolver(ksp)
# Create near null space basis (required for smoothed aggregation
# AMG). The solution vector is passed so that it can be copied to
# generate compatible vectors for the nullspace.
null_space = build_nullspace(V, u.vector())

# Create PETSC smoothed aggregation AMG preconditioner and attach near
# null space
pc = PETScPreconditioner("petsc_amg")
pc.set_nullspace(null_space)
pc.set(solver)

nl_solver.solve(problem, u.vector())

Currently, the PETScKrylovSolver adds the null space from the PETScPreconditioner object to the matP operator context during the call to solve:

https://bitbucket.org/fenics-project/dolfin/src/b67e6aa0645e4fa17eda05df21b2258566685155/dolfin/la/PETScKrylovSolver.cpp?at=master#cl-371

i.e. the PETScPreconditioner object does NOT set the null space during the call to set.

So when you use the SNESKrylovSolver solve function which calls the underlying SNES object's solve call:

https://bitbucket.org/fenics-project/dolfin/src/b67e6aa0645e4fa17eda05df21b2258566685155/dolfin/nls/PETScSNESSolver.cpp?at=master#cl-295

it calls its own solve routine that knows nothing about the DOLFIN PETScPreconditioner object and its nullspace. Note that it does still set the type of the preconditioner and the preconditioners options, so this behaviour is confusing.

To my mind it does not make much sense that the nullspace is set in the PETScKrylovSolver solve method. It should be done earlier. Then the PETScPreconditioner really does set ALL of the attributes of the PETScKrylovSolver.

Thoughts on fixes:

We could fix this by to overloading the PETScPreconditioner.set method to take a PETScSNESSolver, extract the underlying KSP object then assign the nullspace.

Or we could modify the current PETScPreconditioner.set method so that it is responsible for setting the nullspace on the preconditioner.

Any thoughts? My gut feeling is option 2.

I really need this working so I am keen to spend the time to fix it, but I'm looking for a little guidance.

[Jack Hale]

I've managed to reproduce Garth's example using a raw petsc4py KSP object. Still working on a petsc4py SNES version.

"""This demo solves the equations of static linear elasticity for a
pulley subjected to centripetal accelerations. The solver uses
smoothed aggregation algbraric multigrid."""

# Copyright (C) 2014 Garth N. Wells
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.

from __future__ import print_function
import sys

from dolfin import *
import petsc4py
petsc4py.init(comm=mpi_comm_world())
from petsc4py import PETSc

args = [sys.argv[0]] + """
--petsc.pc_type gamg
--petsc.ksp_type gmres
--petsc.ksp_monitor
--petsc.ksp_view
--petsc.ksp_rtol 1E-6
""".split()
parameters.parse(args)


# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
    print("DOLFIN has not been configured with PETSc. Exiting.")
    exit()

# Set backend to PETSC
parameters["linear_algebra_backend"] = "PETSc"

def build_nullspace(V, x):
    """Function to build null space for 3D elasticity"""

    # Create list of vectors for null space
    nullspace_basis = [x.copy() for i in range(6)]

    # Build translational null space basis
    V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
    V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
    V.sub(2).dofmap().set(nullspace_basis[2], 1.0);

    # Build rotational null space basis
    V.sub(0).dofmap().set_x(nullspace_basis[3], -1.0, 1, V.mesh());
    V.sub(1).dofmap().set_x(nullspace_basis[3],  1.0, 0, V.mesh());
    V.sub(0).dofmap().set_x(nullspace_basis[4],  1.0, 2, V.mesh());
    V.sub(2).dofmap().set_x(nullspace_basis[4], -1.0, 0, V.mesh());
    V.sub(2).dofmap().set_x(nullspace_basis[5],  1.0, 1, V.mesh());
    V.sub(1).dofmap().set_x(nullspace_basis[5], -1.0, 2, V.mesh());

    for x in nullspace_basis:
        x.apply("insert")

    return nullspace_basis



# Load mesh and define function space
mesh = Mesh("pulley.xml.gz")

# Function to mark inner surface of pulley
def inner_surface(x, on_boundary):
    r = 3.75 - x[2]*0.17
    return (x[0]*x[0] + x[1]*x[1]) < r*r and on_boundary

# Rotation rate and mass density
omega = 300.0
rho = 10.0

# Loading due to centripetal acceleration (rho*omega^2*x_i)
f = Expression(("rho*omega*omega*x[0]", \
                "rho*omega*omega*x[1]", \
                "0.0"), omega=omega, rho=rho)

# Elasticity parameters
E = 1.0e9
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress computation
def sigma(v):
    gdim = v.geometric_dimension()
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(gdim)

# Create function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), grad(v))*dx
L = inner(f, v)*dx

# Set up boundary condition on inner surface
c = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, c, inner_surface)

# Assemble system, applying boundary conditions and preserving
# symmetry)
A, b = assemble_system(a, L, bc)

# Create solution function
u = Function(V)

# Create near null space basis (required for smoothed aggregation
# AMG). The solution vector is passed so that it can be copied to
# generate compatible vectors for the nullspace.
null_space_petsc = build_nullspace(V, u.vector())
#null_space = VectorSpaceBasis(null_space_petsc)

ksp = PETSc.KSP()
ksp.create(comm=mpi_comm_world())
ksp.setFromOptions()
ksp.setOperators(as_backend_type(A).mat(), as_backend_type(A).mat())

null_space_cast = [as_backend_type(x).vec() for x in null_space_petsc]
null_space_final = PETSc.NullSpace()
null_space_final.create(vectors=null_space_cast)
A, P = ksp.getOperators()
P.setNearNullSpace(null_space_final)

# Compute solution
ksp.solve(as_backend_type(b).vec(), as_backend_type(u.vector()).vec());

# Save solution to VTK format
File("output/elasticity.pvd", "compressed") << u

# Save colored mesh partitions in VTK format if running in parallel
if MPI.size(mesh.mpi_comm()) > 1:
    File("output/partitions.pvd") << CellFunction("size_t", mesh, \
                                           MPI.rank(mesh.mpi_comm()))

# Project and write stress field to post-processing file
W = TensorFunctionSpace(mesh, "Discontinuous Lagrange", 0)
stress = project(sigma(u), V=W)
File("output/stress.pvd") << stress

# Plot solution
#plot(u, interactive=True)

[Jack Hale]

With the new as_backend_type(A).set_nearnullspace(basis) method in 2016.1.0 this may be fixed now, we'll follow up on it soon here at Luxembourg.