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)