- firedrake - imperial.ac.uk

Firedrake on the Isambard HPC system
by Acreman, David 05 Dec '19

05 Dec '19

******************* This email from D.M.Acreman(a)exeter.ac.uk originated outside Imperial. ******************* Hi all, If anyone is interesting in running Firedrake on the Isambard HPC system there are now instructions available at https://github.com/firedrakeproject/isambard (the procedure is similar to setting up Firedrake on Archer). Isambard is one of the EPSRC tier-2 services and will be available to use during the changeover between Archer and Archer-2. Isambard has over 10,000 cores provided by ARM64 processors which have proven particularly effective at running memory-bound applications. If you would like to apply for an account please follow the instructions at https://gw4-isambard.github.io/docs/user-guide/requestaccount.html Regards, Dave Dr. David Acreman, Senior Research Fellow, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Stocker Road, Exeter EX4 4QL

1 0

Assign Initial Value to Mixed Function Spaces
by Thoms 2015 25 Nov '19

25 Nov '19

Dear all, currently, I'm trying to solve a DAE problem (the Poisson-Nernst-Planck system) with 4 variables in Firedrake. The function spaces in the problem are defined as follows: *V0 = FunctionSpace(mesh, "DG",2)* *X = V0*V0*V0*V0* *C0,C1,C2,C3 = Function(X)v0,v1,v2,v3 = TestFunctions(X) * Because the problem is non-linear, I need to prescribe the initial values before running the solver. And the initial values are defined as follows: *ic_C = Constant(1)* *C0.assign(ic_C)* *C1.assign(ic_C)* *C2.assign(ic_C)* *C3.assign(ic_C)* An issue comes up during this step, which says the 'indexed' object (the C0, C1 and so on) has no attribute 'assign': *--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-20-803df79d259a> in <module> 7 ic_Phi = Constant(Phi_init)#project(Phi_0, X) 8 ----> 9 C0.assign(ic_C0); 10 #C1.assign(ic_C1); 11 #C2.assign(ic_C2); AttributeError: 'Indexed' object has no attribute 'assign'* The questions: 1. Do I defined the function spaces correctly? 2. If the Function Spaces are correctly defined, then how to assign the value to the indexed function spaces (the mixed function spaces)? Thank you in advance! Thomas

2 2

VFS
by Cotter, Colin J 23 Nov '19

23 Nov '19

Dear Firedrakers, How do the indices work if you take VectorFunctionSpace(mesh, "BDM", 1) ? Is the first index the dim copies and the second index the vector component of the copies of the BDM spaces? all the best --cjc Professor Colin Cotter (he/him) Department of Mathematics 755, Huxley Building Imperial College London South Kensington Campus United Kingdom of Great Britain and Northern Ireland +44 2075943468

2 3

Don't leave it til the last minute to stock up!
by Adam Martin 18 Nov '19

18 Nov '19

Indigo - essentials for industry FREE CHOCOLATE AND 2020 WALL PLANNER! WITH ALL ORDERS 10% OFF ALL ORDERS UNTIL CHRISTMAS!! [ www.indigo.co ]( http://r.indigosmtp.com/mk/cl/f/Ga7VGa_ixLAdr3DZJAJFp6lC3RH8xb6SDXrpwl_fa24… ) Dear Customer Christmas is just around the corner so now is a great time to check stocks and order ready for the start of 2020! If you require any industrial supplies please do not hesitate to give me a call Until Christmas we are giving away a free 2020 Wall Planner, helping you to have a super organzed and productive 2020! ALL orders also receive a fantastic Chocolate Set to say thank you for your support duriung 2019 We have had an incredible year and look forward to serving you IN 2020 For your reference we supply a massive range of industrial supplies at low, low prices. Products include; Temporary Floor, Carpet and Surface Protection Materials including Self Adhesive Carpet Protection Film and Correx Sheets and Rolls, our dedicated page can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/WUge0VkeH0wYNbgBgBfJlS97X0DHYRS8P2NRM-r9heQ… ) Packaging Supplies including Bubble Wrap, Stretch Film, Packaging Tape, Vermiculite, Document Enclsoed Envelopes etc, our dedicated page can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/mvz52KX76c_OT6ySV_ywxziOWe0_A_5kORZbBhh8CLj… ) Single and Double Sided Adhesive Tapes including Anti Slip Tape, Gaffer Tape, Masking Tape, Double Sided Tapes, Magnetic Tape and more, our dedicated page can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/2-xIrG6dWIypKro6gj24P3Obq6XdSWuOfA49K4_TSw8… ) Hook and Loop Tape and Fasteners, 3M Dual Lock and more can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/yz7hhvGbGMe8tsz4O9nS5UcTPPFKWPZLx_lLfm1beku… ) Asbestos Supplies including Asbestos Sacks can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/4fG_Wd5KQ_tWnX6Jqv_4NyS4f92LUvgpfBsCXm7jwtD… ) Magnets and Magnetic Materials can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/6NlpB7OFTpU8s2a_1afWfw_2ZtxD96roB6nt3mjB1pM… ) Polythene Sheeting including Flame Retardant, Visqueen, Tuffreel and Megafilm can be found [ here ]( http://r.indigosmtp.com/mk/cl/f/21QodMaxgizfc51mamtucUmqKvD4-0A0P7G7TagiFzz… ) We are pleased to let you know that our brand new website is now live! Hundreds of improvements have been made many of which should make it easier for you to find and order the products you need. Check out [ www.indigo.co ]( http://r.indigosmtp.com/mk/cl/f/wDer51WTvwfbKEUxASpADlWo5H4Gu6vwfE453o6mOcd… ) now We also supply a full range of Work Wear, Safety Signs, Swann Morton Blades and Handles, Janitorial Supplies and much more Don't forget all prices are backed up with our price beater guarantee, if you find it cheaper, tell us and we will beat it* please give us a call if you have any requirements or need a sample on 01268 768 768 We hope you have a fantastic Christmas and Happy New Year! Many thanks Adam Martin adam(a)indigo.co T: 01268 768 768 Visit our NEW website at [ www.indigo.co ]( http://r.indigosmtp.com/mk/cl/f/yrUE04KFa37K61LfR8_VFbcnY-CM1OSTCjZBCtqqPto… ) GDPR - We currently have your email on our system, if you wish us to delete this please simply un-subscribe from below or contact adam.martin(a)indigoplc.com and he will get it remove manually and you will receive no further emails in the future. We do not sell your data, you will only receive emails from the Indigo Group but we will remove upon request. Thank you This email was sent to firedrake(a)imperial.ac.uk You received this email because you are registered with indigo. We hope you do not mind the odd email! [ Unsubscribe here ]( http://r.indigosmtp.com/mk/un/snD75NxzXyRC2uiwTfTy05XzL0SdE4VrIFUwD0qWOa4eJ… ) © 2018 Indigo Industrial Supplies Limited

1 0

Preconditioner for the Stokes-Brinkmann problem
by Salazar De Troya, Miguel 24 Oct '19

24 Oct '19

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)

1 0

Point evaluation
by Rodrigo Piccinini 22 Oct '19

22 Oct '19

Hello, I'm trying to evaluate a Function at a point, but I get the error bellow. Could anybody help me on that? The log files are attached. My code can be found here: https://github.com/rbpiccinini/codes/blob/master/eigendrake/spe10_transient… CompilationError: Command "['mpicc', '-std=c99', '-fPIC', '-Wall', > '-march=native', '-O3', '-ffast-math', > '-I/home/piccinini/apps/firedrake/src/firedrake/firedrake', > '-I/home/piccinini/apps/firedrake/include', > '-I/home/piccinini/apps/firedrake/lib/python3.7/site-packages/petsc/include', > '-o', > '/home/piccinini/apps/firedrake/.cache/pyop2/0a/2ff88afc9cdab5dd593941b370e84c_p6855.so.tmp', > '/home/piccinini/apps/firedrake/.cache/pyop2/0a/2ff88afc9cdab5dd593941b370e84c_p6855.c', > '-shared', '-L/home/piccinini/apps/firedrake/lib', '-lspatialindex_c', > '-Wl,-rpath,/home/piccinini/apps/firedrake/lib']" return error status 1. > Unable to compile code > Thanks -- Rodrigo

2 4

Installation on Ubuntu: 18.04.3 LTS
by Bhavesh Shrimali 18 Oct '19

18 Oct '19

Hi all, I tried installing firedrake on my laptop with `Ubuntu: 18.04` but have been unsuccessful so far. These are the warnings I see in the `log` file =============================================================================== ***** WARNING : F90 (set to mpif90) found in environment variables - ignoring use ./configure F90=$F90 if you really want to use that value ****** =============================================================================== =============================================================================== ***** WARNING : Using default optimization C flags -g -O You might consider manually setting optimal optimization flags for your system with COPTFLAGS="optimization flags" see config/examples/arch-*-opt.py for examples =============================================================================== =============================================================================== ***** WARNING : Using default C++ optimization flags -g -O You might consider manually setting optimal optimization flags for your system with CXXOPTFLAGS="optimization flags" see config/examples/arch-*-opt.py for examples =============================================================================== =============================================================================== ***** WARNING : Using default FORTRAN optimization flags -g -O You might consider manually setting optimal optimization flags for your system with FOPTFLAGS="optimization flags" see config/examples/arch-*-opt.py for examples =============================================================================== and it fails to build PETSc beyond this point. ******************************************************************************* UNABLE to CONFIGURE with GIVEN OPTIONS (see configure.log for details): ------------------------------------------------------------------------------- Your libraries are from OpenMPI but it appears your mpiexec is from MPICH ******************************************************************************* I am attaching the entire `log` file for reference. Any leads into what could be possibly going wrong would be helpful. Thanks Bhavesh -- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2(a)illinois.edu <bshrima2(a)illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*

2 1

Talks from Firedrake 19
by Lawrence Mitchell 11 Oct '19

11 Oct '19

Dear all, the Firedrake '19 page (https://firedrakeproject.org/firedrake_19.html) has been updated to add talk slides from some of the talks from attendees. Unfortunately I was not able to make the in-room recording work, and so there is no audio. For those of you who presented and would like your slides to also appear, please send them to me (or provide a link to a website). Thanks, Lawrence

1 0

using different petsc version
by Florian Bruckner 09 Oct '19

09 Oct '19

Dear Firedrake-Team, this may be a question that is very easy to answer, but I didn't get it working by myself. I would like to use a different petsc version which contains some small changes. I am using firedrake inside a Docker container, where i run the firedrake-install script. I tried using PETSC_DIR with --honour-petsc-dir, but this fails with some missing include files. Is PETSC_DIR only the source directory of my own petsc version, or do I have to compile it before? Is it somehow possible to just specify an individual url / git-repository to use a different petsc version? thanks for any advice Florian

2 1

Re: [firedrake] [petsc-users] Stokes-Brinkmann equation preconditioner
by Salazar De Troya, Miguel 07 Oct '19

07 Oct '19

From: Matthew Knepley <knepley(a)gmail.com> Date: Friday, October 4, 2019 at 3:19 AM To: Lawrence Mitchell <wence(a)gmx.li> Cc: "Salazar De Troya, Miguel" <salazardetro1(a)llnl.gov>, "petsc-users(a)mcs.anl.gov" <petsc-users(a)mcs.anl.gov> Subject: Re: [petsc-users] Stokes-Brinkmann equation preconditioner On Fri, Oct 4, 2019 at 6:04 AM Lawrence Mitchell <wence(a)gmx.li<mailto:wence@gmx.li>> wrote: > On 4 Oct 2019, at 10:46, Matthew Knepley via petsc-users <petsc-users(a)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(a)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)

1 0