Troubles setting up PETSc options with dolfinx 0.9.0

dolfinx
1

Hi all,

I’m having trouble migrating a DolfinX v0.6.0 code to DolfinX v0.9.0 (nonlinear NS solver with Schur complement as preconditioning).

This MWE runs under fenics-dolfinx 0.6.0 with petsc4py/petsc 3.19.6 (the conda env has been created from conda create -n fenicsx060 python=3.10 fenics-dolfinx=0.6.0 -c conda-forge):

import numpy as np
import os
import time


#####################################################################
# FEniCS import

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import gmshio
from dolfinx import fem, mesh
from dolfinx.fem.petsc import NonlinearProblem
from ufl import derivative, FacetNormal, TestFunctions, sym, dot, nabla_div, nabla_grad, dx, inner, outer
from dolfinx.nls.petsc import NewtonSolver
import ufl
#####################################################################

def norm_L2(comm, v):
    return np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(inner(v, v) * dx)), op=MPI.SUM))

def g_expr(x):

    return np.vstack((1 - np.exp((Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) * x[0]) * np.cos(2 * np.pi * x[1]),
    (Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) / (2 * np.pi) * np.exp((Re / 2 - np.sqrt(Re**2 / 4 + 4 * np.pi**2)) * x[0]) * np.sin(2 * np.pi * x[1])))

def dir_marker(x):
    return (
        np.isclose(x[0], 0.0)
        | np.isclose(x[0], 1.0)
        | np.isclose(x[1], 0.0)
        | np.isclose(x[1], 1.0)
    )



#############################################################################################################
# VARIABLES

msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)

pts = msh.geometry.x
npts = len(pts)
gdim = msh.geometry.dim

#############################################################################################################
# CONSTANTES

Re = 50
dt = 1e-2

#############################################################################################################
# FEM

def epsilon(u):
	return sym(nabla_grad(u))

# Function spaces for the velocity and for the pressure

Ve = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
Qe = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
W  = fem.FunctionSpace(msh, Ve*Qe)

V, V_map = W.sub(0).collapse()
Q, Q_map = W.sub(1).collapse()

size_V = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
size_Q = Q.dofmap.index_map.size_local * Q.dofmap.index_map_bs

v, q = TestFunctions(W)

wnP1 = fem.Function(W)
wn = fem.Function(W)
bc = fem.Function(W)

unP1, p = ufl.split(wnP1)
un, _ = ufl.split(wn)

A =	( dot(unP1 - un, v) / dt * dx
    - inner(outer(unP1, unP1), nabla_grad(v)) * dx
    + dot(q, nabla_div(unP1)) * dx
    - dot(nabla_div(v), p) * dx
    + inner(epsilon(unP1), epsilon(v))/Re * dx
    )

#############################################################################################################
# Boundaries

u_D = fem.Function(W)
u_D.sub(0).interpolate(g_expr)
dir_points = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, dir_marker)
dir_vel_dofs = fem.locate_dofs_topological(W, msh.topology.dim - 1, dir_points)
bc_u = fem.dirichletbc(u_D, dir_vel_dofs)

bc = [bc_u]



problem = NonlinearProblem(A, wnP1, bcs=bc)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver # Handle for PETSc interaction
ksp.setType("gmres")
pc = ksp.getPC() # Preconditioning object
pc.setType("fieldsplit") # Different preconditioning for different blocks

# Split fields for applying different preconditioners
u_dofs = V_map 
p_dofs = Q_map 
ISu = PETSc.IS().createGeneral(u_dofs, msh.comm)
ISp = PETSc.IS().createGeneral(p_dofs, msh.comm)
pc.setFieldSplitIS(("u",ISu))
pc.setFieldSplitIS(("p",ISp))

# Set preconditioner options
opts = PETSc.Options()
opts.clear()
option_prefix = ksp.getOptionsPrefix()

print(f"Option prefix: {option_prefix}")

opts[f"{option_prefix}ksp_pc_side"] = "right" # Use right preconditioning: ( K M^-1 ) M x = b. This preserves the residual.
opts[f"{option_prefix}pc_fieldsplit_type"] = "schur" # Based on schur complement of K: possible since only two blocks
opts[f"{option_prefix}pc_fieldsplit_schur_fact_type"] = "upper" # Different flavors: diag, lower, upper or full
opts[f"{option_prefix}pc_fieldsplit_schur_preconditioning"] = "user" # Default preconditioner for one of the schur blocks.
opts[f"{option_prefix}fieldsplit_u_ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix}fieldsplit_u_pc_type"] = "jacobi" # Jacobi (diagonals) preconditioning for u-block
opts[f"{option_prefix}fieldsplit_p_ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix}fieldsplit_p_pc_type"] = "hypre" # Hypre (algebraic multigrid) preconditioning for p-block
opts[f"{option_prefix}fieldsplit_p_pc_hypre_type"] = "boomeramg" # Choice of implementation
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_coarsen_type"] = "pmis" # Grid coarsening strategy
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_interp_type"] = "FF1" # Projection on to coarser grid
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_strong_threshold"] = "0.5" # Coarsening limit

ksp.setFromOptions()

# print("Options PETSc définies :")
# for key, value in opts.getAll().items():
#     print(f"{key} = {value}")

start = time.time()
solver.solve(wnP1)
end = time.time()
print(f"Elapsed time: {end - start:.4f} seconds for {npts} points")

and gives the expected output:

Option prefix: nls_solve_
Elapsed time: 1.8762 seconds for 10201 points

When I run this new MWE with fenics-dolfinx 0.9.0 and petsc 3.22.3 / petsc4py 3.22.4 (the conda env has been created using conda create -n fenicsx090 fenics-dolfinx -c conda-forge ):

import numpy as np
import os
import time


#####################################################################
# FEniCS import

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import gmshio
from dolfinx import default_real_type, fem, mesh
from dolfinx.fem.petsc import NonlinearProblem
from ufl import derivative, FacetNormal, TestFunctions, sym, dot, nabla_div, nabla_grad, dx, inner, outer
from basix.ufl import element, mixed_element
from dolfinx.nls.petsc import NewtonSolver
import ufl
#####################################################################

def norm_L2(comm, v):
    return np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(inner(v, v) * dx)), op=MPI.SUM))

def g_expr(x):

    return np.vstack((1 - np.exp((Re.value / 2 - np.sqrt(Re.value**2 / 4 + 4 * np.pi**2)) * x[0]) * np.cos(2 * np.pi * x[1]),
    (Re / 2 - np.sqrt(Re.value**2 / 4 + 4 * np.pi**2)) / (2 * np.pi) * np.exp((Re.value / 2 - np.sqrt(Re.value**2 / 4 + 4 * np.pi**2)) * x[0]) * np.sin(2 * np.pi * x[1])))

def dir_marker(x):
    return (
        np.isclose(x[0], 0.0)
        | np.isclose(x[0], 1.0)
        | np.isclose(x[1], 0.0)
        | np.isclose(x[1], 1.0)
    )



#############################################################################################################
# VARIABLES

msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)

pts = msh.geometry.x
npts = len(pts)
gdim = msh.geometry.dim

#############################################################################################################
# CONSTANTES

Re = fem.Constant(msh, default_real_type(50))
dt = fem.Constant(msh, default_real_type(1e-2))

#############################################################################################################
# FEM

def epsilon(u):
	return sym(nabla_grad(u))

# Function spaces for the velocity and for the pressure
V_el = element("Lagrange", msh.basix_cell(), 2, shape=(gdim,))
Q_el = element("Lagrange", msh.basix_cell(), 1)
W_el = mixed_element([V_el, Q_el])

W = fem.functionspace(msh, W_el)
V, V_map = W.sub(0).collapse()
Q, Q_map = W.sub(1).collapse()

size_V = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
size_Q = Q.dofmap.index_map.size_local * Q.dofmap.index_map_bs

v, q = TestFunctions(W)

wnP1 = fem.Function(W)
wn = fem.Function(W)
bc = fem.Function(W)

unP1, p = ufl.split(wnP1)
un, _ = ufl.split(wn)

A =	( dot(unP1 - un, v) / dt * dx
    - inner(outer(unP1, unP1), nabla_grad(v)) * dx
    + dot(q, nabla_div(unP1)) * dx
    - dot(nabla_div(v), p) * dx
    + inner(epsilon(unP1), epsilon(v))/Re * dx
    )

#############################################################################################################
# Boundaries

u_D = fem.Function(W)
u_D.sub(0).interpolate(g_expr)
dir_points = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, dir_marker)
dir_vel_dofs = fem.locate_dofs_topological(W, msh.topology.dim - 1, dir_points)
bc_u = fem.dirichletbc(u_D, dir_vel_dofs)

bc = [bc_u]



problem = NonlinearProblem(A, wnP1, bcs=bc)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver # Handle for PETSc interaction
ksp.setType("gmres")
pc = ksp.getPC() # Preconditioning object
pc.setType("fieldsplit") # Different preconditioning for different blocks

# Split fields for applying different preconditioners
u_dofs = V_map 
p_dofs = Q_map 
ISu = PETSc.IS().createGeneral(u_dofs, msh.comm)
ISp = PETSc.IS().createGeneral(p_dofs, msh.comm)
pc.setFieldSplitIS(("u",ISu))
pc.setFieldSplitIS(("p",ISp))

# Set preconditioner options
opts = PETSc.Options()
opts.clear()
option_prefix = ksp.getOptionsPrefix()

print(f"Option prefix: {option_prefix}")

opts[f"{option_prefix}ksp_pc_side"] = "right" # Use right preconditioning: ( K M^-1 ) M x = b. This preserves the residual.
opts[f"{option_prefix}pc_fieldsplit_type"] = "schur" # Based on schur complement of K: possible since only two blocks
opts[f"{option_prefix}pc_fieldsplit_schur_fact_type"] = "upper" # Different flavors: diag, lower, upper or full
opts[f"{option_prefix}pc_fieldsplit_schur_precondition"] = "user" # Default preconditioner for one of the schur blocks.
opts[f"{option_prefix}fieldsplit_u_ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix}fieldsplit_u_pc_type"] = "jacobi" # Jacobi (diagonals) preconditioning for u-block
opts[f"{option_prefix}fieldsplit_p_ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix}fieldsplit_p_pc_type"] = "hypre" # Hypre (algebraic multigrid) preconditioning for p-block
opts[f"{option_prefix}fieldsplit_p_pc_hypre_type"] = "boomeramg" # Choice of implementation
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_coarsen_type"] = "pmis" # Grid coarsening strategy
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_interp_type"] = "FF1" # Projection on to coarser grid
opts[f"{option_prefix}fieldsplit_p_pc_hypre_boomeramg_strong_threshold"] = "0.5" # Coarsening limit

ksp.setFromOptions()

# print("Options PETSc définies :")
# for key, value in opts.getAll().items():
#     print(f"{key} = {value}")

start = time.time()
solver.solve(wnP1)
end = time.time()
print(f"Elapsed time: {end - start:.4f} seconds for {npts} points")

The output specifies that the same options are unused:

Option prefix: nls_solve_
Elapsed time: 5.0500 seconds for 10201 points
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 12 unused database options. They are:
Option left: name:-nls_solve_fieldsplit_p_ksp_type value: preonly source: code
Option left: name:-nls_solve_fieldsplit_p_pc_hypre_boomeramg_coarsen_type value: pmis source: code
Option left: name:-nls_solve_fieldsplit_p_pc_hypre_boomeramg_interp_type value: FF1 source: code
Option left: name:-nls_solve_fieldsplit_p_pc_hypre_boomeramg_strong_threshold value: 0.5 source: code
Option left: name:-nls_solve_fieldsplit_p_pc_hypre_type value: boomeramg source: code
Option left: name:-nls_solve_fieldsplit_p_pc_type value: hypre source: code
Option left: name:-nls_solve_fieldsplit_u_ksp_type value: preonly source: code
Option left: name:-nls_solve_fieldsplit_u_pc_type value: jacobi source: code
Option left: name:-nls_solve_ksp_pc_side value: right source: code
Option left: name:-nls_solve_pc_fieldsplit_schur_fact_type value: upper source: code
Option left: name:-nls_solve_pc_fieldsplit_schur_preconditioning value: user source: code
Option left: name:-nls_solve_pc_fieldsplit_type value: schur source: code

I believe this is related to the PETSc versions, but from what I’ve seen on their website, the prefix format hasn’t changed, so I’m giving it a try here.

The solution given in this topic:

works only on the first four options.

Thanks!