<html xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.msonormal0, li.msonormal0, div.msonormal0
{mso-style-name:msonormal;
mso-margin-top-alt:auto;
margin-right:0in;
mso-margin-bottom-alt:auto;
margin-left:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style>
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><b><span style="font-size:12.0pt;color:black">From: </span></b><span style="font-size:12.0pt;color:black">Matthew Knepley <knepley@gmail.com><br>
<b>Date: </b>Friday, October 4, 2019 at 3:19 AM<br>
<b>To: </b>Lawrence Mitchell <wence@gmx.li><br>
<b>Cc: </b>"Salazar De Troya, Miguel" <salazardetro1@llnl.gov>, "petsc-users@mcs.anl.gov" <petsc-users@mcs.anl.gov><br>
<b>Subject: </b>Re: [petsc-users] Stokes-Brinkmann equation preconditioner<o:p></o:p></span></p>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="margin-left:.5in">On Fri, Oct 4, 2019 at 6:04 AM Lawrence Mitchell <<a href="mailto:wence@gmx.li">wence@gmx.li</a>> wrote:<o:p></o:p></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<p class="MsoNormal" style="margin-left:.5in"><br>
<br>
> On 4 Oct 2019, at 10:46, Matthew Knepley via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> <br>
> On Thu, Oct 3, 2019 at 6:34 PM Salazar De Troya, Miguel via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> 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:<a href="https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426" target="_blank">https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426</a>
(section 2.6) <br>
> <br>
> <br>
> Link does not work for me.<br>
<br>
Try <a href="https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426" target="_blank">
https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426</a> (mail clients are awful)<br>
<br>
> <br>
> 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.<br>
> <br>
> <br>
> <br>
> Because I am solving for velocity and pressure without doing any substitution, I thought I could use a PC fieldsplit type Schur (full factorization)<br>
> <br>
> <br>
> Yes, this will form the exact factorization and a matrix-free form of the Schur complement.<br>
> <br>
> and provide the preconditioner in the paper to solve the Schur complement.<br>
> <br>
> <br>
> Yes, you can provide a user-defined PC for the Schur complement.<br>
> <br>
> My question is, should I provide this preconditioner through PCFieldSplitSetSchurPre or through fieldsplit_1_pc_type (probably through the Firedrake interface as inhttps://<a href="http://www.firedrakeproject.org/demos/stokes.py.html" target="_blank">www.firedrakeproject.org/demos/stokes.py.html</a>)
?<br>
> <br>
> <br>
> The name PCFieldSplitSetSchurPre seems to be very misleading. You do not use it to provide a _preconditioner_. You use it to determine<br>
> the _preconditioning matrix_ from which the actual preconditioner is built. The preconditioner itself is defined using -fieldpslit_1_pc_type.<br>
> Since I do not know what the preconditioner looks like, I cannot say what preconditioner matrix you would want. Since Firedrake can construct<br>
> any operator for you, you might not care about the matrix we pass to you.<br>
<br>
The action of the PC is defined by:<br>
<br>
M^{-1} q = \phi_q + \mu q<br>
<br>
\mu is some number, and \phi_q solves the BVP<br>
<br>
-div \alpha^{-1} grad \phi_q = q<br>
<br>
grad \phi_q \dot n = 0<br>
\int \phi_q = 0<br>
<br>
So your PC for the fieldsplit_1 needs to:<br>
<br>
1. Discretise and solve this BVP to compute \phi_q<br>
2. Add \mu q.<br>
<br>
I.e. the action of PCApply is:<br>
<br>
y <- \mu x + Laplace^{-1} x<br>
<br>
(Aside, surely a primal-dual error has been made in the analysis here)<br>
<br>
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).<o:p></o:p></p>
</blockquote>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in">Richardson R(x, y) is<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"> y += mu (b - Ax)<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in">so if b = 0 and A = I, you get that. It seems easier just to make a PCSHELL with VecAXPY I guess.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"> Thanks,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"> Matt<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:.5in"> <o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<p class="MsoNormal" style="mso-margin-top-alt:0in;margin-right:0in;margin-bottom:12.0pt;margin-left:.5in">
Cheers,<br>
<br>
Lawrence<br>
<br>
<o:p></o:p></p>
</blockquote>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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 <a href="https://www.firedrakeproject.org/_modules/firedrake/preconditioners/pcd.html">
https://www.firedrakeproject.org/_modules/firedrake/preconditioners/pcd.html</a> . 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”<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks<o:p></o:p></p>
<p class="MsoNormal">Miguel<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">from firedrake import *<o:p></o:p></p>
<p class="MsoNormal">from firedrake.petsc import PETSc<o:p></o:p></p>
<p class="MsoNormal">from firedrake.preconditioners.base import PCBase<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">class BrinkmannPC(PCBase):<o:p></o:p></p>
<p class="MsoNormal"> r""" Preconditioner for the Stokes-Brinkmann problem<o:p></o:p></p>
<p class="MsoNormal"> Fill details<o:p></o:p></p>
<p class="MsoNormal"> """<o:p></o:p></p>
<p class="MsoNormal"> def initialize(self, pc):<o:p></o:p></p>
<p class="MsoNormal"> from firedrake import TrialFunction, TestFunction, dx, \<o:p></o:p></p>
<p class="MsoNormal"> assemble, inner, grad, split, Constant, parameters<o:p></o:p></p>
<p class="MsoNormal"> from firedrake.assemble import allocate_matrix, create_assembly_callable<o:p></o:p></p>
<p class="MsoNormal"> if pc.getType() != "python":<o:p></o:p></p>
<p class="MsoNormal"> raise ValueError("Expecting PC type python")<o:p></o:p></p>
<p class="MsoNormal"> prefix = pc.getOptionsPrefix() + "brink_"<o:p></o:p></p>
<p class="MsoNormal"> # we assume P has things stuffed inside of it<o:p></o:p></p>
<p class="MsoNormal"> _, P = pc.getOperators()<o:p></o:p></p>
<p class="MsoNormal"> context = P.getPythonContext()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> test, trial = context.a.arguments()<o:p></o:p></p>
<p class="MsoNormal"> if test.function_space() != trial.function_space():<o:p></o:p></p>
<p class="MsoNormal"> raise ValueError("Pressure space test and trial space differ")<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> Q = test.function_space()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> p = TrialFunction(Q)<o:p></o:p></p>
<p class="MsoNormal"> q = TestFunction(Q)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> alpha = context.appctx.get("alpha", Constant(1.0))<o:p></o:p></p>
<p class="MsoNormal"> stiffness = inner(1.0 / alpha * grad(p), grad(q))*dx<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> opts = PETSc.Options()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> # we're inverting Kp, so default them to assembled.<o:p></o:p></p>
<p class="MsoNormal"> # Mp only needs its action, so default it to mat-free.<o:p></o:p></p>
<p class="MsoNormal"> # These can of course be overridden.<o:p></o:p></p>
<p class="MsoNormal"> default = parameters["default_matrix_type"]<o:p></o:p></p>
<p class="MsoNormal"> Kp_mat_type = opts.getString(prefix+"Kp_mat_type", default)<o:p></o:p></p>
<p class="MsoNormal"> self.Mp_mat_type = opts.getString(prefix+"Mp_mat_type", "matfree")<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> Kp = assemble(stiffness, form_compiler_parameters=context.fc_params,<o:p></o:p></p>
<p class="MsoNormal"> mat_type=Kp_mat_type,<o:p></o:p></p>
<p class="MsoNormal"> options_prefix=prefix + "Kp_")<o:p></o:p></p>
<p class="MsoNormal"> Kp_petsc = Kp.M.handle<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> L = q * dx<o:p></o:p></p>
<p class="MsoNormal"> self.rhs = assemble(L)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> from firedrake import VectorSpaceBasis, Function<o:p></o:p></p>
<p class="MsoNormal"> ctsp = Function(Q).interpolate(Constant(1.0))<o:p></o:p></p>
<p class="MsoNormal"> <span lang="ES">nullspace = VectorSpaceBasis(vecs=[ctsp])<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="ES"> </span>nullspace.orthonormalize()<o:p></o:p></p>
<p class="MsoNormal"> nullspace_petsc = nullspace.nullspace()<o:p></o:p></p>
<p class="MsoNormal"> Kp_petsc.setNearNullSpace(nullspace_petsc)<o:p></o:p></p>
<p class="MsoNormal"> self.nullspace = nullspace<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> Kksp = PETSc.KSP().create(comm=pc.comm)<o:p></o:p></p>
<p class="MsoNormal"> Kksp.incrementTabLevel(1, parent=pc)<o:p></o:p></p>
<p class="MsoNormal"> Kksp.setOptionsPrefix(prefix + "Kp_")<o:p></o:p></p>
<p class="MsoNormal"> Kksp.setOperators(Kp_petsc)<o:p></o:p></p>
<p class="MsoNormal"> Kksp.setFromOptions()<o:p></o:p></p>
<p class="MsoNormal"> Kksp.setUp()<o:p></o:p></p>
<p class="MsoNormal"> self.Kksp = Kksp<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> mu = context.appctx.get("mu", Constant(1.0))<o:p></o:p></p>
<p class="MsoNormal"> mass = mu*p*q*dx<o:p></o:p></p>
<p class="MsoNormal"> self.Mp = allocate_matrix(mass, form_compiler_parameters=context.fc_params,<o:p></o:p></p>
<p class="MsoNormal"> mat_type=self.Mp_mat_type,<o:p></o:p></p>
<p class="MsoNormal"> options_prefix=prefix + "Mp_")<o:p></o:p></p>
<p class="MsoNormal"> self._assemble_Mp = create_assembly_callable(mass, tensor=self.Mp,<o:p></o:p></p>
<p class="MsoNormal"> form_compiler_parameters=context.fc_params,<o:p></o:p></p>
<p class="MsoNormal"> mat_type=self.Mp_mat_type)<o:p></o:p></p>
<p class="MsoNormal"> self._assemble_Mp()<o:p></o:p></p>
<p class="MsoNormal"> Mpmat = self.Mp.petscmat<o:p></o:p></p>
<p class="MsoNormal"> self.workspace = [Mpmat.createVecLeft() for i in (0, 1)]<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def update(self, pc):<o:p></o:p></p>
<p class="MsoNormal"> pass<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def apply(self, pc, x, y):<o:p></o:p></p>
<p class="MsoNormal"> a, b = self.workspace<o:p></o:p></p>
<p class="MsoNormal"> with self.rhs.dat.vec_wo as rhs_dat:<o:p></o:p></p>
<p class="MsoNormal"> x.copy(rhs_dat)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> self.nullspace.orthogonalize(self.rhs)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> with self.rhs.dat.vec_wo as rhs_dat:<o:p></o:p></p>
<p class="MsoNormal"> self.Kksp.solve(rhs_dat, a)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> self.Mp.petscmat.mult(x, y)<o:p></o:p></p>
<p class="MsoNormal"> y.axpy(1.0, a)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def applyTranspose(self, pc, x, y):<o:p></o:p></p>
<p class="MsoNormal"> a, b = self.workspace<o:p></o:p></p>
<p class="MsoNormal"> with self.rhs.dat.vec_wo as rhs_dat:<o:p></o:p></p>
<p class="MsoNormal"> x.copy(rhs_dat)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> self.nullspace.orthogonalize(self.rhs)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> with self.rhs.dat.vec_wo as rhs_dat:<o:p></o:p></p>
<p class="MsoNormal"> self.Kksp.solveTranspose(rhs_dat, a)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> self.Mp.petscmat.mult(x, y)<o:p></o:p></p>
<p class="MsoNormal"> y.axpy(1.0, a)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def view(self, pc, viewer=None):<o:p></o:p></p>
<p class="MsoNormal"> super(BrinkmannPC, self).view(pc, viewer)<o:p></o:p></p>
<p class="MsoNormal"> viewer.printfASCII("KSP solver for K^-1:\n")<o:p></o:p></p>
<p class="MsoNormal"> self.Kksp.view(viewer)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">N = 100<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">mesh = UnitSquareMesh(N, N)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">V = VectorFunctionSpace(mesh, "CG", 2)<o:p></o:p></p>
<p class="MsoNormal">W = FunctionSpace(mesh, "CG", 1)<o:p></o:p></p>
<p class="MsoNormal">Z = V * W<o:p></o:p></p>
<p class="MsoNormal">print("# DOFS {}".format(Z.dim()))<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">w = Function(Z)<o:p></o:p></p>
<p class="MsoNormal">u, p = split(w)<o:p></o:p></p>
<p class="MsoNormal">v, q = TestFunctions(Z)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">mu = Constant(1.0)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">x, y = SpatialCoordinate(mesh)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">from ufl import And<o:p></o:p></p>
<p class="MsoNormal"><span lang="ES">alpha = conditional(And(y > 0.4, And(x < 0.2, y < 0.8)), Constant(1e5), Constant(1e5))<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="ES">alpha = Constant(1e3)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="ES">alpha = 1e2*exp((0.2 - (x - 0.5)*(x - 0.5) - (y - 0.5)*(y - 0.5)) / 0.05)<o:p></o:p></span></p>
<p class="MsoNormal">alpha_interp = Function(W)<o:p></o:p></p>
<p class="MsoNormal">alpha_interp.interpolate(alpha)<o:p></o:p></p>
<p class="MsoNormal">File("alpha.pvd").write(alpha_interp)<o:p></o:p></p>
<p class="MsoNormal"><span lang="ES">F = mu*inner(grad(u), grad(v))*dx - p*div(v)*dx - div(u)*q*dx + alpha*inner(u, v)*dx<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="ES"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="ES"><o:p> </o:p></span></p>
<p class="MsoNormal">bc_value = as_vector([0, x * (x - 1)])<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">bcs = [DirichletBC(Z.sub(0), bc_value, (4,)),<o:p></o:p></p>
<p class="MsoNormal"> DirichletBC(Z.sub(0), zero(mesh.geometric_dimension()), (1, 2)),<o:p></o:p></p>
<p class="MsoNormal"> DirichletBC(Z.sub(0), bc_value, (3,))]<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">def create_solver(solver_parameters, appctx=None):<o:p></o:p></p>
<p class="MsoNormal"> p = {}<o:p></o:p></p>
<p class="MsoNormal"> if solver_parameters is not None:<o:p></o:p></p>
<p class="MsoNormal"> p.update(solver_parameters)<o:p></o:p></p>
<p class="MsoNormal"> # Default to linear SNES<o:p></o:p></p>
<p class="MsoNormal"> p.setdefault("snes_type", "ksponly")<o:p></o:p></p>
<p class="MsoNormal"> p.setdefault("ksp_atol", 1e-6)<o:p></o:p></p>
<p class="MsoNormal"> problem = NonlinearVariationalProblem(F, w, bcs=bcs)<o:p></o:p></p>
<p class="MsoNormal"> solver = NonlinearVariationalSolver(problem,<o:p></o:p></p>
<p class="MsoNormal"> solver_parameters=p, appctx=appctx)<o:p></o:p></p>
<p class="MsoNormal"> return solver<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">fieldsplit_0_lu = {<o:p></o:p></p>
<p class="MsoNormal"> "ksp_type": "preonly",<o:p></o:p></p>
<p class="MsoNormal"> "ksp_max_it": 1,<o:p></o:p></p>
<p class="MsoNormal"> "pc_type": "lu",<o:p></o:p></p>
<p class="MsoNormal"> "pc_factor_mat_solver_type": "mumps",<o:p></o:p></p>
<p class="MsoNormal">}<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">fieldsplit_1_lu = fieldsplit_0_lu<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">fieldsplit_0_python = {<o:p></o:p></p>
<p class="MsoNormal"> "ksp_type": "preonly",<o:p></o:p></p>
<p class="MsoNormal"> "pc_type": "python",<o:p></o:p></p>
<p class="MsoNormal"> "pc_python_type": "firedrake.AssembledPC",<o:p></o:p></p>
<p class="MsoNormal"> "assembled_pc_type": "hypre",<o:p></o:p></p>
<p class="MsoNormal">}<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">fieldsplit_1_python = {<o:p></o:p></p>
<p class="MsoNormal"> "ksp_type" : "cg",<o:p></o:p></p>
<p class="MsoNormal"> "ksp_atol" : 1e-2,<o:p></o:p></p>
<p class="MsoNormal"> "pc_type": "python",<o:p></o:p></p>
<p class="MsoNormal"> "pc_python_type": "__main__.BrinkmannPC",<o:p></o:p></p>
<p class="MsoNormal"> "ksp_converged_reason" : None,<o:p></o:p></p>
<p class="MsoNormal"> "brink_Kp_ksp_type": "preonly",<o:p></o:p></p>
<p class="MsoNormal"> "brink_Kp_pc_type": "lu",<o:p></o:p></p>
<p class="MsoNormal"> "brink_Kp_pc_factor_mat_solver_type" : "mumps",<o:p></o:p></p>
<p class="MsoNormal">}<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">iterative = True<o:p></o:p></p>
<p class="MsoNormal">params = {<o:p></o:p></p>
<p class="MsoNormal"> "mat_type" : "matfree" if iterative else "aij",<o:p></o:p></p>
<p class="MsoNormal"> "ksp_monitor_true_residual": None,<o:p></o:p></p>
<p class="MsoNormal"> "ksp_converged_reason": None,<o:p></o:p></p>
<p class="MsoNormal"> "ksp_max_it" : 100,<o:p></o:p></p>
<p class="MsoNormal"> "ksp_atol" : 1e-6,<o:p></o:p></p>
<p class="MsoNormal"> "ksp_type" : "fgmres",<o:p></o:p></p>
<p class="MsoNormal"> "pc_type" : "fieldsplit",<o:p></o:p></p>
<p class="MsoNormal"> "pc_fieldsplit_type" : "schur",<o:p></o:p></p>
<p class="MsoNormal"> "pc_fieldsplit_schur_fact_type": "full",<o:p></o:p></p>
<p class="MsoNormal"> "pc_fieldsplit_schur_precondition": "selfp" if not iterative else "a11",<o:p></o:p></p>
<p class="MsoNormal"> "fieldsplit_0": fieldsplit_0_python if iterative else fieldsplit_0_lu,<o:p></o:p></p>
<p class="MsoNormal"> "fieldsplit_1" : fieldsplit_1_python if iterative else fieldsplit_1_lu,<o:p></o:p></p>
<p class="MsoNormal">}<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">w.assign(0)<o:p></o:p></p>
<p class="MsoNormal">appctx = {"alpha" : alpha, "mu" : mu}<o:p></o:p></p>
<p class="MsoNormal">solver = create_solver(params, appctx=appctx)<o:p></o:p></p>
<p class="MsoNormal"><span lang="ES">solver.solve()<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="ES">u, p = w.split()<o:p></o:p></span></p>
<p class="MsoNormal">u.rename("velocity")<o:p></o:p></p>
<p class="MsoNormal">p.rename("pressure")<o:p></o:p></p>
<p class="MsoNormal">File("out.pvd").write(u, p)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:.5in"><o:p> </o:p></p>
</div>
</div>
</body>
</html>