From: Matthew Knepley <knepley@gmail.com>
Date: Friday, October 4, 2019 at 3:19 AM
To: Lawrence Mitchell <wence@gmx.li>
Cc: "Salazar De Troya, Miguel" <salazardetro1@llnl.gov>, "petsc-users@mcs.anl.gov" <petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] Stokes-Brinkmann equation preconditioner
On Fri, Oct 4, 2019 at 6:04 AM Lawrence Mitchell <wence@gmx.li> wrote:
> On 4 Oct 2019, at 10:46, Matthew Knepley via petsc-users <petsc-users@mcs.anl.gov> wrote:
>
> On Thu, Oct 3, 2019 at 6:34 PM Salazar De Troya, Miguel via petsc-users <petsc-users@mcs.anl.gov> wrote:
> 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)
>
>
> Link does not work for me.
Try https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426 (mail clients are awful)
>
> 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)
>
>
> Yes, this will form the exact factorization and a matrix-free form of the Schur complement.
>
> and provide the preconditioner in the paper to solve the Schur complement.
>
>
> Yes, you can provide a user-defined PC for the Schur complement.
>
> My question is, should I provide this preconditioner through PCFieldSplitSetSchurPre or through fieldsplit_1_pc_type (probably through the Firedrake interface as inhttps://www.firedrakeproject.org/demos/stokes.py.html) ?
>
>
> The name PCFieldSplitSetSchurPre seems to be very misleading. You do not use it to provide a _preconditioner_. You use it to determine
> the _preconditioning matrix_ from which the actual preconditioner is built. The preconditioner itself is defined using -fieldpslit_1_pc_type.
> Since I do not know what the preconditioner looks like, I cannot say what preconditioner matrix you would want. Since Firedrake can construct
> any operator for you, you might not care about the matrix we pass to you.
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
So your 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
(Aside, surely a primal-dual error has been made in the analysis here)
You could, I think, do this by providing the discretisation of this laplacian and using an additive PCComposite (although I don't know what PC you would use to just scale the input, it's easy to write one though).
Richardson R(x, y) is
y += mu (b - Ax)
so if b = 0 and A = I, you get that. It seems easier just to make a PCSHELL with VecAXPY I guess.
Thanks,
Matt
Cheers,
Lawrence
Why is there a need to use a PCComposite instead of applying “y <- \mu x + Laplace^{-1} x” directly in the same preconditioner? I have tried to do the latter without success, maybe that is why. I have used Firedrake and followed the example
in
https://www.firedrakeproject.org/_modules/firedrake/preconditioners/pcd.html . I hope it is ok to include it here given that it uses a lot of the petsc4py routines. The code is attached below. The implementation of the preconditioner for the Schur complement
must be wrong somewhere because 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. I hope it is okay to ask here if there is something wrong in how I am applying “y <- \mu x + Laplace^{-1} x”
Thanks
Miguel
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)
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"]
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
ctsp = Function(Q).interpolate(Constant(1.0))
nullspace = VectorSpaceBasis(vecs=[ctsp])
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(1e5), Constant(1e5))
alpha = Constant(1e3)
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):
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)
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}
solver = create_solver(params, appctx=appctx)
solver.solve()
u, p = w.split()
u.rename("velocity")
p.rename("pressure")
File("out.pvd").write(u, p)