Application of GAMG preconditioner for singular block

Linear Algebra
1

Hi all,

we are currently trying to create a block-preconditioner where the blocks to be preconditioned are pure Neumann Laplacians. With the code below we try to incorporate the null space, but the convergence of MINRES is weird as can be seen from the code snippet below.
We would appreciate any help in setting the null space correctly or finding out, which options to pass to the PETSc solver. Also, we are not sure how to orthogonalize the RHS correctly as nsp.orthogonalize(b1) does not seem to work in dolfinx.

Thanks, Kai

    null_vec = b1.copy() # A1_block.createVecLeft()
    null_vec.array[:] = np.sqrt(1.0/null_vec.getSize())
    null_vec.normalize()
    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    A1_block.setNullSpace(nsp)
    #nsp.orthogonalize(b1)

print('pol, type:', ksp1.getType(), 'residual:', ksp1.getResidualNorm(), 'iterations: ', ksp1.getIterationNumber(), '\n\tconvergence reason', ksp1.getConvergedReason(), 'manually computed residual:', res1.norm())
pol, type: minres residual: 7603541464.314032 iterations:  2 
	convergence reason 2 manually computed residual: 0.03956033832717597
2

As you have not supplied a minimal code example to reproduce the issue, it is quite hard for anyone to give you any guidance.

Could you make a minimal code example that backs up this statement?

3

Thanks for the quick reply! I post the minimal example below. The line

nsp.orthogonalize(b1)

raises

AttributeError: 'petsc4py.PETSc.NullSpace' object has no attribute 'orthogonalize'

Minimal example:

import numpy as np
from dolfinx.fem import Function, FunctionSpace, assemble_matrix_block, assemble_vector_block, form
from dolfinx.mesh import CellType, create_rectangle
from ufl import FiniteElement, TrialFunction, TestFunction, variable, diff, dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc

# Parameters
epsilon = 1.0e-03
dt = 1.0e-05
t = 0
T = 1e-04

# Mesh and function spaces
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([10, 2.5])], [200, 100], CellType.triangle)
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
MEc, MEmu = FunctionSpace(mesh, P1), FunctionSpace(mesh, P1)

# Functions
dc1, dmu1 = TrialFunction(MEc), TrialFunction(MEmu)
q1, v1 = TestFunction(MEc), TestFunction(MEmu)
c1, c01 = Function(MEc), Function(MEc)
mu1, mu01 = Function(MEmu), Function(MEmu)

# Initial conditions
mu1.x.array[:] = 0.0
c1.interpolate(lambda x: 0.35 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
c1.x.scatter_forward()

def setup_potential(c01):
    # Potential
    c01 = variable(c01)
    f = 3.5 * c01**2 * (1 - c01)**2
    dfdc01 = diff(f, c01)
    return dfdc01

def assemble_IMEX(dt, dc1, q1, dmu1, v1, c01, dfdc01):
    # Weak form
    Lone_block = form([[dt*inner(grad(dmu1),grad(q1))*dx, inner(dc1, q1)*dx], 
        [inner(dmu1, v1)*dx, - epsilon*inner(grad(dc1),grad(v1))*dx]])
    fone_block = form([inner(c01, q1)*dx, inner(dfdc01, v1)*dx])
    
    # Assemble linear system
    A1_block = assemble_matrix_block(Lone_block)
    A1_block.assemble()
    b1 = assemble_vector_block(fone_block, Lone_block)
    
    # Preconditioner
    a_p_11 = form(epsilon*inner(grad(dc1),grad(v1))*dx)
    a_p = [[Lone_block[0][0], None],
       [None, a_p_11]]
    P1 = assemble_matrix_block(a_p)
    P1.assemble()
    
    ##### Code snippet 1
    # Nullspace
    null_vec = b1.copy() # A1_block.createVecLeft()
    null_vec.array[:] = np.sqrt(1.0/null_vec.getSize())
    null_vec.normalize()
    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    A1_block.setNullSpace(nsp)
    #nsp.orthogonalize(b1)
    #####
    
    # Krylov solver
    ksp1 = PETSc.KSP()
    ksp1.create(PETSc.COMM_WORLD)
    #ksp1.setType("minres")
    ksp1.setTolerances(rtol=1e-06)
    ksp1.setOperators(A1_block, P1)
    ksp1.getPC().setType("gamg")
    
    return ksp1, b1, A1_block


dfdc01 = setup_potential(c01)

n = len(c1.vector.array)

# Loop over time
while t < T:
    c01.x.array[:] = c1.x.array
    mu01.x.array[:] = mu1.x.array
    
    ksp1, b1, A1 = assemble_IMEX(dt, dc1, q1, dmu1, v1, c01, dfdc01)
    
    # Solve linear system
    x1 = A1.createVecRight()
    ksp1.solve(b1, x1)
    
    res1 = A1*x1 - b1
    
    ##### Code snippet 2
    print('pol, type:', ksp1.getType(), 'residual:', ksp1.getResidualNorm(), 'iterations: ', ksp1.getIterationNumber(), '\n\tconvergence reason', ksp1.getConvergedReason(), 'manually computed residual:', res1.norm())
    #####
    
    t += dt
    
    mu1.x.array[:] = x1.array_r[:n]
    c1.x.array[:] = x1.array_r[n:2*n]
5

See for instance: dolfinx/test_krylov_solver.py at main · FEniCS/dolfinx · GitHub

1 Like