Dear PETSc team,
We are solving the Stokes equations using PETSc (via Firedrake) on a highly
anisotropic 3D domain (L_x=1, L_y=0.01, L_z=0.1).
In this setup, standard Schur complement preconditioners using a mass inverse
for pressure struggle to converge. We could solve the problem with the LSC
preconditioner (solver parameters are shown in the script).
We have the following questions:
* Why standard preconditioners struggle in such domains?
* Why is the preconditioned residual norm for the Schur complement system
much higher than the true residual norm?
* Would you recommend alternative or more robust preconditioners for such
geometries?
Thank you for your help.
Best regards,
Hardik
HARDIK KOTHARI
hardik.koth...@corintis.com<mailto:hardik.koth...@corintis.com>
Corintis SA
EPFL Innovation Park Building C
1015 Lausanne
[https://urldefense.us/v3/__https://storcor.s3.eu-central-1.amazonaws.com/logos/Logo-black.png__;!!G_uCfscf7eWS!coZIkzX5gHWdyGx54LSoh7iP6KZ7AV-0Oa7FTFrMVAzIYrBCsVJ0vqY0hLvcIqcrdQtxcmcdpr25Hhsper7CRJwNrHrjfDsq$
]
Here at Corintis we care for your privacy. That is why we have taken
appropriate measures to ensure that the data you have provided to us is always
secure
import firedrake as fd
from firedrake import inner, grad, div, dx
from firedrake.petsc import PETSc
simulation_case = "elongated"
# simulation_case = "cube"
if simulation_case == "elongated":
Lx = 1
Ly = 0.01
Lz = 0.1
nx = 100
ny = int(nx * Ly / Lx * 5)
nz = int(nx * Lz / Lx * 5)
else:
Lx = 1
Ly = 1
Lz = 1
nx = 20
ny = nx
nz = nx
mesh = fd.BoxMesh(nx, ny, nz, Lx, Ly, Lz)
V = fd.VectorFunctionSpace(mesh, "CG", 2)
W = fd.FunctionSpace(mesh, "CG", 1)
Z = V * W
PETSc.Sys.Print(f"DOFs: {Z.dim()}")
u, p = fd.TrialFunctions(Z)
v, q = fd.TestFunctions(Z)
def a(u, v, p, q):
return (inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx
L = inner(fd.Constant((0, 0, 0)), v) * dx
max_vel = 16
x, y, z = fd.SpatialCoordinate(mesh)
u_inlet = fd.as_vector(
[
max_vel * (1 - y / Ly) * y / Ly * z / Lz * (1 - z / Lz),
0,
0,
]
)
bcs = [
fd.DirichletBC(Z.sub(0), u_inlet, (1,)),
fd.DirichletBC(Z.sub(0), fd.Constant((0, 0, 0)), (3, 4, 5, 6)),
]
up = fd.Function(Z)
a_form = a(u, v, p, q)
direct_solver_settings = {
"ksp_type": "preonly",
"ksp_monitor_true_residual": None,
"ksp_monitor": None,
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_14": 200,
"mat_mumps_icntl_24": 1,
}
class MassInvPC(fd.AssembledPC):
_prefix = "Mp_"
def form(self, pc, test, trial):
_, bcs = super(MassInvPC, self).form(pc)
appctx = self.get_appctx(pc)
mu = appctx.get("mu", 1.0)
a = inner((1 / mu) * trial, test) * dx
return a, bcs
def set_nullspaces(self, pc):
# the mass matrix does not have a nullspace
pass
assembled_hypre_settings = {
"ksp_type": "cg",
# "ksp_converged_reason": None,
"ksp_max_it": 100,
"ksp_atol": 1e-15,
"ksp_rtol": 1e-15,
"pc_type": "python",
"pc_python_type": "firedrake.AssembledPC",
"assembled_pc_type": "hypre",
"assembled_pc_hypre_type": "boomeramg",
"assembled_pc_hypre_boomeramg_coarsen_type": "PMIS",
"assembled_pc_hypre_boomeramg_interp_type": "ext+i",
"assembled_pc_hypre_boomeramg_strong_threshold": 0.1,
"assembled_pc_hypre_boomeramg_agg_nl": 1,
"assembled_pc_hypre_boomeramg_relax_type_all": "l1-Gauss-Seidel",
"assembled_pc_hypre_boomeramg_smooth_type": "ILU",
"assembled_pc_hypre_boomeramg_truncfactor": 0.3,
"assembled_pc_hypre_boomeramg_max_levels": 25,
"assembled_pc_hypre_boomeramg_max_iter": 1,
}
mass_inv_settings = {
"ksp_type": "gmres",
"ksp_converged_reason": None,
"ksp_monitor_true_residual": None,
"ksp_max_it": 100,
"ksp_atol": 1e-12,
"ksp_rtol": 1e-12,
"pc_type": "python",
"pc_python_type": "__main__.MassInvPC",
"Mp_pc_type": "lu",
}
split1_lsc_iterative = {
"ksp_type": "gmres",
"ksp_monitor_true_residual": None,
"ksp_rtol": 1e-4,
"ksp_gmres_modifiedgramschmidt": None,
"ksp_converged_reason": None,
"pc_type": "lsc",
"lsc_pc_type": "cholesky",
"pc_lsc_scale_diag": True,
}
solver_parameters = {
"ksp_type": "fgmres",
"ksp_atol": 1e-7,
"ksp_rtol": 1e-6,
"ksp_monitor_true_residual": None,
"pc_type": "fieldsplit",
"pc_fieldsplit_type": "schur",
"pc_fieldsplit_schur_fact_type": "full",
"pc_fieldsplit_schur_precondition": "selfp",
"fieldsplit_0": assembled_hypre_settings,
"fieldsplit_1": split1_lsc_iterative,
}
up.assign(0)
fd.solve(a_form == L, up, bcs=bcs, solver_parameters=solver_parameters)
u, p = up.subfunctions
u.rename("Velocity")
p.rename("Pressure")
fd.VTKFile("slimStokes.pvd").write(u, p)
