Hello,
I am restating this question here for better context. Last email was a bit sloppy with all the stuff coming from the petsc mail list.
I am trying to solve the Stokes equation with the Brinkman term to simulate a solid material. My intention is to implement the preconditioner in this paper: https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426
(section 2.6) where they solve for the velocity and substitute that expression in the pressure equation. They end up solving a system of the type B K^-1 B^T, i.e. the Schur complement of the problem. For this system of equations, they argue that the preconditioner
in page 11 is perfect for a given constant Brinkman penalty term.
Because I am solving for velocity and pressure without doing any substitution, I thought I could use a PC fieldsplit type Schur (full factorization) and provide the preconditioner in the paper to solve the
Schur complement. In the paper, the action of the PC is defined by:
M^{-1} q = \phi_q + \mu q
\mu is some number, and \phi_q solves the BVP
-div \alpha^{-1} grad \phi_q = q
grad \phi_q \dot n = 0
\int \phi_q = 0
Therefore, the PC for the fieldsplit_1 needs to:
1. Discretise and solve this BVP to compute \phi_q
2. Add \mu q.
I.e. the action of PCApply is:
y <- \mu x + Laplace^{-1} x
I followed the example in https://www.firedrakeproject.org/_modules/firedrake/preconditioners/pcd.html
to try to implement this preconditioner . The code is attached below. The fieldsplit_1_ solve takes many iterations to converge whenever a Gaussian bump or
a jump function models the \alpha parameter. It behaves better for \alpha constant (an assumption in the paper for the preconditioner to be perfect). I am solving the pure Neumann problem by passing the nullspace to the KSP. Are there mistakes in the way I
am implementing the preconditioner? Or maybe the lack of convergence is just purely because of the method?
Thanks
Miguel
Miguel A. Salazar de Troya
Postdoctoral Researcher, Lawrence Livermore National Laboratory
B141
Rm: 1085-5
Ph: 1(925) 422-6411
from firedrake import *
from firedrake.petsc import PETSc
from firedrake.preconditioners.base import PCBase
class BrinkmannPC(PCBase):
r""" Preconditioner for the Stokes-Brinkmann problem
Fill details
"""
def initialize(self, pc):
from firedrake import TrialFunction, TestFunction, dx, \
assemble, inner, grad, split, Constant, parameters
from firedrake.assemble import allocate_matrix, create_assembly_callable
if pc.getType() != "python":
raise ValueError("Expecting PC type python")
prefix = pc.getOptionsPrefix() + "brink_"
# we assume P has things stuffed inside of it
_, P = pc.getOperators()
context = P.getPythonContext()
test, trial = context.a.arguments()
if test.function_space() != trial.function_space():
raise ValueError("Pressure space test and trial space differ")
Q = test.function_space()
p = TrialFunction(Q)
q = TestFunction(Q)
# Regularisation to avoid having to think about nullspaces.
alpha = context.appctx.get("alpha", Constant(1.0))
stiffness = inner(1.0 / alpha * grad(p), grad(q))*dx
opts = PETSc.Options()
# we're inverting Kp, so default them to assembled.
# Mp only needs its action, so default it to mat-free.
# These can of course be overridden.
default = parameters["default_matrix_type"]
#default = "aij"
Kp_mat_type = opts.getString(prefix+"Kp_mat_type", default)
self.Mp_mat_type = opts.getString(prefix+"Mp_mat_type", "matfree")
Kp = assemble(stiffness, form_compiler_parameters=context.fc_params,
mat_type=Kp_mat_type,
options_prefix=prefix + "Kp_")
Kp_petsc = Kp.M.handle
L = q * dx
self.rhs = assemble(L)
from firedrake import VectorSpaceBasis, Function
x_axis = Function(Q).interpolate(Constant(1.0))
nullspace = VectorSpaceBasis(vecs=[x_axis])
nullspace.orthonormalize()
nullspace_petsc = nullspace.nullspace()
Kp_petsc.setNearNullSpace(nullspace_petsc)
self.nullspace = nullspace
Kksp = PETSc.KSP().create(comm=pc.comm)
Kksp.incrementTabLevel(1, parent=pc)
Kksp.setOptionsPrefix(prefix + "Kp_")
Kksp.setOperators(Kp_petsc)
Kksp.setFromOptions()
Kksp.setUp()
self.Kksp = Kksp
mu = context.appctx.get("mu", Constant(1.0))
mass = mu*p*q*dx
self.Mp = allocate_matrix(mass, form_compiler_parameters=context.fc_params,
mat_type=self.Mp_mat_type,
options_prefix=prefix + "Mp_")
self._assemble_Mp = create_assembly_callable(mass, tensor=self.Mp,
form_compiler_parameters=context.fc_params,
mat_type=self.Mp_mat_type)
self._assemble_Mp()
Mpmat = self.Mp.petscmat
self.workspace = [Mpmat.createVecLeft() for i in (0, 1)]
def update(self, pc):
pass
def apply(self, pc, x, y):
a, b = self.workspace
with self.rhs.dat.vec_wo as rhs_dat:
x.copy(rhs_dat)
self.nullspace.orthogonalize(self.rhs)
with self.rhs.dat.vec_wo as rhs_dat:
self.Kksp.solve(rhs_dat, a)
self.Mp.petscmat.mult(x, y)
y.axpy(1.0, a)
def applyTranspose(self, pc, x, y):
a, b = self.workspace
with self.rhs.dat.vec_wo as rhs_dat:
x.copy(rhs_dat)
self.nullspace.orthogonalize(self.rhs)
with self.rhs.dat.vec_wo as rhs_dat:
self.Kksp.solveTranspose(rhs_dat, a)
self.Mp.petscmat.mult(x, y)
y.axpy(1.0, a)
def view(self, pc, viewer=None):
super(BrinkmannPC, self).view(pc, viewer)
viewer.printfASCII("KSP solver for K^-1:\n")
self.Kksp.view(viewer)
N = 100
mesh = UnitSquareMesh(N, N)
V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1)
Z = V * W
print("# DOFS {}".format(Z.dim()))
w = Function(Z)
u, p = split(w)
v, q = TestFunctions(Z)
mu = Constant(1.0)
x, y = SpatialCoordinate(mesh)
from ufl import And
alpha = conditional(And(y > 0.4, And(x < 0.2, y < 0.8)), Constant(1e2), Constant(1e-1))
#alpha = Constant(1e1)
#alpha = 1e2*exp((0.2 - (x - 0.5)*(x - 0.5) - (y - 0.5)*(y - 0.5)) / 0.05)
alpha_interp = Function(W)
alpha_interp.interpolate(alpha)
File("alpha.pvd").write(alpha_interp)
F = mu*inner(grad(u), grad(v))*dx - p*div(v)*dx - div(u)*q*dx + alpha*inner(u, v)*dx
bc_value = as_vector([0, x * (x - 1)])
bcs = [DirichletBC(Z.sub(0), bc_value, (4,)),
DirichletBC(Z.sub(0), zero(mesh.geometric_dimension()), (1, 2)),
DirichletBC(Z.sub(0), bc_value, (3,))]
def create_solver(solver_parameters, appctx=None, nullspace=nullspace):
p = {}
if solver_parameters is not None:
p.update(solver_parameters)
# Default to linear SNES
p.setdefault("snes_type", "ksponly")
p.setdefault("ksp_atol", 1e-6)
problem = NonlinearVariationalProblem(F, w, bcs=bcs)
solver = NonlinearVariationalSolver(problem,
solver_parameters=p, appctx=appctx, nullspace=nullspace)
return solver
fieldsplit_0_lu = {
"ksp_type": "preonly",
"ksp_max_it": 1,
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
fieldsplit_1_lu = fieldsplit_0_lu
fieldsplit_0_python = {
"ksp_type": "preonly",
"pc_type": "python",
"pc_python_type": "firedrake.AssembledPC",
"assembled_pc_type": "hypre",
}
fieldsplit_1_python = {
"ksp_type" : "cg",
"ksp_atol" : 1e-2,
"pc_type": "python",
"pc_python_type": "__main__.BrinkmannPC",
"ksp_converged_reason" : None,
"brink_Kp_ksp_type": "preonly",
"brink_Kp_pc_type": "lu",
"brink_Kp_pc_factor_mat_solver_type" : "mumps",
}
iterative = True
params = {
"mat_type" : "matfree" if iterative else "aij",
"ksp_monitor_true_residual": None,
"ksp_converged_reason": None,
"ksp_max_it" : 100,
"ksp_atol" : 1e-6,
"ksp_type" : "fgmres",
"pc_type" : "fieldsplit",
"pc_fieldsplit_type" : "schur",
"pc_fieldsplit_schur_fact_type": "full",
"pc_fieldsplit_schur_precondition": "selfp" if not iterative else "a11",
"fieldsplit_0": fieldsplit_0_python if iterative else fieldsplit_0_lu,
"fieldsplit_1" : fieldsplit_1_python if iterative else fieldsplit_1_lu,
}
w.assign(0)
appctx = {"alpha" : alpha, "mu" : mu}
nullspace = MixedVectorSpaceBasis(
Z, [Z.sub(0), VectorSpaceBasis(constant=True)])
solver = create_solver(params, appctx=appctx, nullspace=nullspace)
solver.solve()
u, p = w.split()
u.rename("velocity")
p.rename("pressure")
File("brinkman_problem.pvd").write(u, p)