One step forward, two steps back.
I managed to get both the linear and nonlinear solver to apply the preconditioner by starting the option setting code with: ksp.setOptionsPrefix(""). No idea why.
Everything seems to work, until I run with mpi and more than 1 core. For a tiny problem it says:
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 55, Out of memory
So I would still much appreciate any help…
Here is the updated code snippet:
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)
ksp = newton_solver.krylov_solver
opts = PETSc.Options()
ksp.setOptionsPrefix("") # This fixes the issue with the nonlinear 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
IS_u = PETSc.IS().createGeneral(WV_map, domain.comm).sort()
IS_p = PETSc.IS().createGeneral(WQ_map, domain.comm).sort()
pc = ksp.getPC() # Preconditioning object
pc.setType("fieldsplit") # Different preconditioning for different blocks
pc.setFieldSplitIS(('u',IS_u),('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.")