Hello,
After exploring the idea of approximating the schur complement with an augmented lagrangein method, I came across a grad-div implementation in freefem (GitHub - prj-/moulin2019al: Augmented Lagrangian Preconditioner for Hydrodynamic Stability Analysis). I have the impression that this can be converted to fenicsx. In fact, using splitfield multiplicative we obtain the lower block triangular part of the saddle point matrix. Then we add Sp -1 = -(γ + Re-1)Mp-1
in the pressure block (null) to obtain our preconditioner (using PCFieldSplitGetSubKSP and KSPSetOperators
).
However, I get the error
"0 SNES Function norm 1.142225359199e+01
Linear ns_ solve did not converge due to DIVERGED_PC_FAILED iterations 0
PC failed due to SUBPC_ERROR ".
I’ve tried replacing the velocity and pressure block solvers with direct LU solvers, but nothing changes, even if I remove the Sp matrix. I think it’s coming from the multiplicative pc? I’d like to point out that with the schur +Selfp method or by doing a direct LU without gmres, the code works just fine.
Do you have any idea where this is coming from?
Here is my simple, usable code:
from mpi4py import MPI
import dolfinx
import time
from ufl import split, TestFunctions, TestFunction, derivative, adjoint, action, Identity, grad, transpose, inner, dx, replace, div, Measure, TrialFunction, sqrt, dot
import ufl
import basix.ufl
from petsc4py import PETSc
from petsc4py.PETSc import Options
from dolfinx.fem.petsc import NonlinearProblem, apply_lifting, assemble_vector, set_bc, assemble_matrix, create_matrix, create_vector
from dolfinx.fem import Function, functionspace, Constant, dirichletbc, locate_dofs_topological, locate_dofs_geometrical, form, assemble_scalar
from dolfinx.mesh import create_box, create_rectangle, locate_entities_boundary, CellType
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import mesh
import numpy as np
from dolfinx.io.gmshio import model_to_mesh
import scipy
import gmsh
comm = MPI.COMM_WORLD
gmsh.initialize()
tags = {}
tags["inlet"] = 1
tags["outlet"] = 2
tags["wall"] = 3
model = gmsh.model
model_rank = 0
comm = MPI.COMM_WORLD
if comm.rank == model_rank:
factory = model.occ
s1 = factory.addRectangle(-1, 0, 0, 1, 1)
s2 = factory.addRectangle(0, -1, 0, 7.5, 2)
s3 = 3
factory.fuse([(2, s1)], [(2, s2)], tag=s3)
factory.synchronize()
ov = model.getBoundary([(2, s3)])
l1, l2, l3, l4, l5, l6 = (val[1] for val in ov) # counterclockwise (l6 = inlet)
# Tag interior
model.addPhysicalGroup(2, [s3], tag=0)
# Tag boundaries
model.addPhysicalGroup(1, [l6], tag=tags["inlet"])
model.setPhysicalName(1, tags["inlet"], "inlet")
model.addPhysicalGroup(1, [l4], tag=tags["outlet"])
model.setPhysicalName(1, tags["outlet"], "outlet")
model.addPhysicalGroup(1, [l1, l2, l3, l5], tag=tags["wall"])
model.setPhysicalName(1, tags["wall"], "wall")
# Set uniform mesh size at all points
mesh_size = 0.025
model.mesh.setSize(model.getEntities(0), mesh_size)
# Generate 2D mesh
model.mesh.generate(dim=2)
# Convert Gmsh mesh to DOLFINx
mesh, _, facet_tags = model_to_mesh(model, MPI.COMM_WORLD, rank=model_rank, gdim=2)
facet_tags.name = "facet_tags"
gmsh.finalize()
solver = "non lu"
precond = "mal"
x = ufl.SpatialCoordinate(mesh)
fdim = mesh.topology.dim-1
metadata = {"quadrature_degree": 4}
dx = Measure("dx", domain=mesh, metadata=metadata)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + mesh.topology.index_map(mesh.topology.dim).num_ghosts
cell_indices = np.arange(num_cells, dtype=np.int32)
h_values = mesh.h(dim=mesh.topology.dim, entities=cell_indices); h_avg = np.mean(h_values)
v_max_inlet = Constant(mesh, PETSc.ScalarType(1))
cfl = 3
dt = cfl * h_avg / v_max_inlet
k = Constant(mesh, PETSc.ScalarType(dt))
Re = 100
rho_ = Constant(mesh, PETSc.ScalarType(1))
mu_ = Constant(mesh, PETSc.ScalarType(rho_*v_max_inlet*2/Re))
nu_ = Constant(mesh, PETSc.ScalarType(mu_/rho_))
V_element = basix.ufl.element("CG", mesh.basix_cell(), 2, shape = (mesh.topology.dim,))
P_element = basix.ufl.element("CG", mesh.basix_cell(), 1)
U_element = basix.ufl.mixed_element([V_element,P_element])
U = functionspace(mesh, U_element); U0 = U.sub(0);P0 = U.sub(1)
w_v, w_p = TestFunctions(U)
uh = Function(U)
v, p = split(uh)
uh_prev = Function(U)
v_prev, p_prev = split(uh_prev)
V_v = U.sub(0).collapse()[0]
def v_inflow_eval(x):
values = np.zeros((2, x.shape[1]))
values[0] = 4.0 * x[1] * (1.0 - x[1])
values[1] = np.zeros(x.shape[1])
return values
def v0_expr(x):
values = np.zeros((mesh.topology.dim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = v_max_inlet # v_x = v_max_inlet, v_y = 0, v_z = 0
return values
v_wall = Function(V_v)
v_wall.x.array[:] = 0.0
v_inflow = Function(V_v)
v_inflow.interpolate(v_inflow_eval)
# Localize boundary dofs for velocity
fdim = mesh.topology.dim - 1
wall_dofs_v = locate_dofs_topological((U0,V_v), fdim, facet_tags.find(tags["wall"]))
inlet_dofs_v = locate_dofs_topological((U0,V_v), fdim, facet_tags.find(tags["inlet"]))
bcs = [dirichletbc(v_inflow, inlet_dofs_v, U0), dirichletbc(v_wall, wall_dofs_v, U0)]
#uh.sub(0).interpolate(v0_expr)
#uh.sub(1).interpolate(lambda x: np.zeros(x.shape[1]))
F_PV = (-div(v) * w_p * dx + nu_ * inner(grad(v),grad(w_v)) * dx + inner(grad(v)*v, w_v) * dx - p * div(w_v) * dx)
if precond == "mal":
gamma_ = 1.0
F_PV += - gamma_*div(v)*div(w_v) *dx
class NonlinearPDE_SNESProblem:
def __init__(self, F, u, bc):
V = u.function_space
du = TrialFunction(V)
self.L = form(F)
self.a = form(derivative(F, u, du))
self.bc = bc
self._F, self._J = None, None
self.u = u
def F(self, snes, x, F):
"""Assemble residual vector."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
with F.localForm() as f_local:
f_local.set(0.0)
assemble_vector(F, self.L)
apply_lifting(F, [self.a], bcs=[self.bc], x0=[x], alpha=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, self.bc, x, -1.0)
def J(self, snes, x, J, P):
"""Assemble Jacobian matrix."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
J.zeroEntries()
assemble_matrix(J, self.a, bcs=self.bc)
J.assemble()
#Snes solver
problem = NonlinearPDE_SNESProblem(F_PV, uh, bcs)
J = create_matrix(problem.a)
b = create_vector(problem.L)
# Create Newton solver and solve
snes = PETSc.SNES().create()
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
problem_prefix = "ns_"
snes.setOptionsPrefix(problem_prefix)
opts = PETSc.Options()
opts.prefixPush(problem_prefix)
opts["snes_type"] = "newtonls"
opts['snes_linesearch_type'] = 'bt'
opts['snes_linesearch_order']= 2
opts['snes_monitor'] = None
if solver == "Direct":
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
opts.prefixPop() # ns_
snes.setFromOptions()
else :
U_v, map_u = U.sub(0).collapse() # sous-espace vitesse + indices
U_p, map_p = U.sub(1).collapse()
U_indexmap = U.dofmap.index_map
localdofs_u = np.setdiff1d( U_indexmap.local_to_global(np.array(map_u)), U_indexmap.ghosts ) # getting the global velocity dofs that live on this core
localdofs_p = np.setdiff1d( U_indexmap.local_to_global(np.array(map_p)), U_indexmap.ghosts ) # getting the global pressure dofs that live on this core
snes.getKSP().setType("fgmres")
snes.getKSP().getPC().setType("fieldsplit")
opts["ksp_max_it"] = 2000
opts["ksp_converged_reason"] = None
opts["ksp_rtol"] = 1e-3
#opts["ksp_monitor"] = None
opts["ksp_gmres_restart"] = 200
opts["ksp_pc_side"] = "right"
ISu = PETSc.IS().createGeneral(np.array(localdofs_u,dtype=np.int32), comm=mesh.comm)
ISp = PETSc.IS().createGeneral(np.array(localdofs_p,dtype=np.int32), comm=mesh.comm)
snes.getKSP().getPC().setFieldSplitIS(("u", ISu))
snes.getKSP().getPC().setFieldSplitIS(("p", ISp))
if precond == "Simple":
opts["pc_fieldsplit_type"] = "schur"
opts["pc_fieldsplit_schur_fact_type"]= "full"
opts["pc_fieldsplit_schur_precondition"]= "selfp"
opts["fieldsplit_u_ksp_type"]= "preonly"
opts["fieldsplit_u_pc_type"]= "hypre"
opts["fieldsplit_u_pc_hypre_type"]= "boomeramg"
opts["fieldsplit_u_pc_hypre_boomeramg_no_CF"]= None
opts["fieldsplit_u_pc_hypre_boomeramg_coarsen_type"]= "HMIS"
opts["fieldsplit_u_pc_hypre_boomeramg_interp_type"]= "ext+i"
opts["fieldsplit_u_pc_hypre_boomeramg_P_max"]= 4
opts["fieldsplit_u_pc_hypre_boomeramg_agg_nl"]= 1
opts["fieldsplit_u_pc_hypre_boomeramg_agg_num_paths"]= 2
opts["fieldsplit_p_ksp_type"]= "preonly"
opts["fieldsplit_p_pc_type"]= "jacobi"
opts.prefixPop() # ns_
snes.setFromOptions()
if precond == "mal":
opts["pc_fieldsplit_type"] = "multiplicative"
opts["fieldsplit_u_ksp_type"] = "gmres"
opts["fieldsplit_u_ksp_pc_side"] = "right"
opts["fieldsplit_u_ksp_rtol"] = "1e-1"
opts["fieldsplit_u_ksp_gmres_restart"] = 50
opts["fieldsplit_u_pc_type"] = "asm"
opts["fieldsplit_u_pc_asm_overlap"] = "1"
opts["fieldsplit_u_sub_pc_type"] = "lu"
opts["fieldsplit_u_sub_pc_factor_mat_solver_type"] = "mumps"
#opts["fieldsplit_u_ksp_type"] = "preonly"
#opts["fieldsplit_u_pc_type"] = "lu"
#opts["fieldsplit_u_pc_factor_mat_solver_type"] = "mumps"
opts["fieldsplit_p_ksp_type"] = "cg"
opts["fieldsplit_p_ksp_max_it"] = "5"
opts["fieldsplit_p_pc_type"] = "jacobi"
#opts["fieldsplit_p_ksp_type"] = "preonly"
#opts["fieldsplit_p_pc_type"] = "lu"
#opts["fieldsplit_p_pc_factor_mat_solver_type"] = "mumps"
opts.prefixPop()
snes.setFromOptions()
#snes.setUp() # construit SNES et KSP principal
#snes.getKSP().setUp() # construit le KSP global
#snes.getKSP().getPC().setUp() # construit le PC FieldSplit
subksp = snes.getKSP().getPC().getFieldSplitSubKSP()
vel_ksp, pres_ksp = subksp
q_p, r_p = ufl.TrialFunction(U_p), ufl.TestFunction(U_p)
S_form = form (-1/(gamma_ + 1/Re) * q_p*r_p *dx)
S = assemble_matrix(S_form)
S.assemble()
pres_ksp.setOperators(S)
snes.getKSP().getPC().view()
x = uh.x.petsc_vec.copy()
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
t0 = time.time()
snes.solve(None, x)
t1 = time.time()
if comm.rank ==0:
print("Temps écoulé: ", t1-t0)
x.copy(uh.x.petsc_vec)
uh.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
uh.x.scatter_forward()
snes.getKSP().getPC().view()