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<mailto:wence@gmx.li>> wrote:

On 4 Oct 2019, at 10:46, Matthew Knepley via petsc-users <petsc-users@mcs.anl.gov<mailto: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<mailto: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<http://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)