Hi all,
For a Navier-Stokes problem, I want to separately precondition the u block (Jacobi) and the p block (algebraic multigrid on the Schur complement), while using mixed spaces.
I believe I have everything set up correctly, and am (finally 
) no longer getting the PETSc warning āThere are unused database optionsā. But when I activate solver verbosity (setting ksp_view), then none of my own settings are mentioned in the output. Instead it mentions LU or MUMPS preconditioners (depending on mpi run with 1 or multiple cores).
Here is a MWE of a lid driven cavity:
from dolfinx import mesh, fem
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import ufl
from ufl import div, nabla_grad, dot, inner, outer, dx, sym
from basix.ufl import element, mixed_element
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
# User parameters
Re = 10 # Reynolds number
N = 5 # Elements in x & y
Pu = 2 # Polynomial order u
Pp = 1 # Polynomial order p
# Function space
L = 1 # Domain size
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(0,0),(L,L)], [N,N], mesh.CellType.quadrilateral)
Ve = element("Lagrange", domain.basix_cell(), Pu, shape=(domain.geometry.dim,))
Qe = element("Lagrange", domain.basix_cell(), Pp)
W_el = mixed_element([Ve, Qe])
W = fem.functionspace(domain, W_el)
V, WV_map = W.sub(0).collapse() # DOF map from W to V
Q, WQ_map = W.sub(1).collapse() # DOF map from W to Q
# Trial and test functions
w_trial = ufl.TrialFunction(W)
v, q = ufl.TestFunctions(W)
u_p = fem.Function(W)
u, p = ufl.split(u_p)
# Navier-Stokes weak form
nu = fem.Constant(domain, PETSc.ScalarType(1/Re)) # Viscosity
B = - inner( outer(u,u), nabla_grad(v) )*dx \
+ dot( q, div(u) )*dx \
- dot( div(v), p )*dx \
+ inner( 2*nu*sym(nabla_grad(u)), sym(nabla_grad(v)) )*dx
# Boundary conditions
def node_origin(x):
return np.logical_and(np.isclose(x[0],0),np.isclose(x[1],0))
dof_origin_p = fem.locate_dofs_geometrical((W.sub(1), Q), node_origin)
bc_zero_p = fem.dirichletbc(fem.Function(Q), dof_origin_p, W.sub(1))
def walls(x):
return np.isclose(x[0],0) | np.isclose(x[0], L) | np.isclose(x[1],0)
facets_wall = mesh.locate_entities_boundary(domain, domain.geometry.dim-1, walls)
dofs_walls = fem.locate_dofs_topological((W.sub(0), V), domain.geometry.dim-1, facets_wall)
bc_zero_u = fem.dirichletbc(fem.Function(V), dofs_walls, W.sub(0))
def lid(x):
return np.isclose(x[1], L)
facets_lid = mesh.locate_entities_boundary(domain, domain.geometry.dim-1, lid)
dofs_lid = fem.locate_dofs_topological((W.sub(0), V), domain.geometry.dim-1, facets_lid)
u_unit = fem.Function(V)
u_unit.interpolate(lambda x: (np.ones(x.shape[1]),np.zeros(x.shape[1])) )
bc_unit_u = fem.dirichletbc(u_unit, dofs_lid, W.sub(0))
bcs = [bc_zero_p,bc_zero_u,bc_unit_u]
##########
# Below is where the solver and its settings are
##########
J = ufl.derivative(B,u_p,w_trial)
ns_problem = NonlinearProblem(B,u_p,bcs=bcs,J=J)
newton_solver = NewtonSolver(domain.comm,ns_problem)
opts = PETSc.Options()
ksp = newton_solver.krylov_solver
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
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}ksp_view"] = None
pc = ksp.getPC() # Preconditioning object
pc.setType("fieldsplit") # Different preconditioning for different blocks
IS_u = PETSc.IS().createGeneral(WV_map, domain.comm)
pc.setFieldSplitIS(('u',IS_u))
IS_p = PETSc.IS().createGeneral(WQ_map, domain.comm)
pc.setFieldSplitIS(('p',IS_p))
ksp_u, ksp_p = pc.getFieldSplitSubKSP()
pc.setFieldSplitType(PETSc.PC.CompositeType.SCHUR)
pc.setFieldSplitSchurFactType(PETSc.PC.FieldSplitSchurFactType.UPPER)
pc.setFieldSplitSchurPreType(PETSc.PC.FieldSplitSchurPreType.SELFP)
option_prefix_u = ksp_u.getOptionsPrefix()
opts[f"{option_prefix_u}ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix_u}pc_type"] = "jacobi"
option_prefix_p = ksp_p.getOptionsPrefix()
opts[f"{option_prefix_p}ksp_type"] = "preonly" # Precondition only once in Newton iterations
opts[f"{option_prefix_p}pc_type"] = "hypre" # Hypre (algebraic multigrid) preconditioning for p-block
opts[f"{option_prefix_p}pc_hypre_type"] = "boomeramg" # Choice of implementation
opts[f"{option_prefix_p}pc_hypre_boomeramg_coarsen_type"] = "pmis" # Grid coarsening strategy
opts[f"{option_prefix_p}pc_hypre_boomeramg_interp_type"] = "FF1" # Projection on to coarser grid
opts[f"{option_prefix_p}pc_hypre_boomeramg_strong_threshold"] = "0.5" # Coarsening limit
ksp_u.setFromOptions()
ksp_p.setFromOptions()
ksp.setFromOptions()
##########
# Above is where the solver and its settings are
##########
it,_ = newton_solver.solve(u_p)
print(f"Solved in {it} iterations.")
And this is the output for a single-core run:
KSP Object: (nls_solve_) 1 MPI process
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zero
tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
right preconditioning
using UNPRECONDITIONED norm type for convergence test
PC Object: (nls_solve_) 1 MPI process
type: lu
out-of-place factorization
tolerance for zero pivot 2.22045e-14
matrix ordering: nd
factor fill ratio given 5., needed 1.47955
Factored matrix follows:
Mat Object: (nls_solve_) 1 MPI process
type: seqaij
rows=278, cols=278
package used to perform factorization: petsc
total: nonzeros=14328, allocated nonzeros=14328
using I-node routines: found 103 nodes, limit used is 5
linear system matrix = precond matrix:
Mat Object: 1 MPI process
type: seqaij
rows=278, cols=278
total: nonzeros=9684, allocated nonzeros=9684
total number of mallocs used during MatSetValues calls=0
using I-node routines: found 142 nodes, limit used is 5
Any help would be much appreciated!
Best,
Stein
