Firedrake Documentation

Firedrake is an automated system designed for solving partial differential equations (PDEs) using the finite element method (FEM). It leverages sophisticated code generation to allow mathematicians, scientists, and engineers to create high-performance simulations with high productivity.

Key capabilities include:

• PDE Specification: Use the Unified Form Language (UFL) from the FEniCS Project to express any PDE.
• Solvers: Seamless coupling with PETSc for programmable solvers.
• Meshing: Support for unstructured meshes (triangular, quadrilateral, and tetrahedral) as well as layered meshes (triangular wedges or hexahedra).
• Finite Elements: A vast range of finite element spaces.
• Optimization: Automatic optimization including vectorization and sum factorization for high-order elements.
• Preconditioning: Support for geometric multigrid and customizable operator preconditioners.
• Advanced Methods: Support for static condensation, hybridisation, and HDG methods.

For multi-field problems like the Stokes equations, geometric multigrid can be composed with PETSc's fieldsplit preconditioning. This allows you to use different solvers for different fields (e.g., using MG to invert the velocity block in a Schur complement system).

Using Auxiliary Operators in Multigrid

When using fieldsplit with a Schur complement, you may need to provide a bilinear form for an auxiliary operator (like a pressure mass matrix) that does not appear in the original system. You can do this by subclassing AuxiliaryOperatorPC and implementing the form method.

Coupled Geometric Multigrid

# Example: Subclassing AuxiliaryOperatorPC for a Schur complement
class Mass(AuxiliaryOperatorPC):
    def form(self, pc, test, trial):
        a = 1/nu * inner(test, trial)*dx
        bcs = None
        return (a, bcs)

# Using it in solver parameters
parameters = {
    "ksp_type": "gmres",
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "schur",
    "fieldsplit_1_pc_type": "python",
    "fieldsplit_1_pc_python_type": "geometric_multigrid.Mass",
    # ... other parameters ...
}

To achieve Reynolds-robustness (where error estimates and Krylov iteration counts are independent of the Reynolds number $\mathrm{Re}$ and mesh size), use a strategy combining three components:

1. H(div)-conforming discretisation: Use an element pair like Brezzi–Douglas–Marini (BDM) for velocity and Discontinuous Galerkin (DG) for pressure. This ensures the divergence-free constraint is captured exactly at the discrete level.
2. Augmented Lagrangian technique: Add a penalty term $\frac{\gamma}{2} \int_\Omega (\nabla \cdot u)^2 ,dx$ to the Lagrangian to control the pressure Schur complement. For large $\gamma$, the Schur complement is driven to a scalar multiple of the pressure mass matrix $Q$, which is easily invertible.
3. Parameter-robust geometric multigrid: Use a multigrid preconditioner based on vertex-star space decomposition to handle the augmented velocity block $A_\gamma$ efficiently.

This approach is particularly effective for the stationary incompressible Navier-Stokes equations where classical solvers like PCD or LSC degrade as $\mathrm{Re}$ increases.

UFL distinguishes between 'primal' quantities (objects in a space $V$) and 'dual' quantities (objects in the dual space $V^*$).

• Functions/Coefficients: A Function (or Coefficient) represents a known function in the primal space. When integrated (e.g., f = c * dx), it results in a scalar Form.
• Arguments/Trial/Test Functions: Argument, TrialFunction, and TestFunction represent unknown functions in the primal space.
• Cofunctions: When you assemble() a linear form involving an unknown (like TrialFunction(V) * dx), the resulting object is a Cofunction. A Cofunction is the dual equivalent of a Coefficient and represents a known object in the dual space.
• Coarguments: To represent an unknown object in the dual space, use a Coargument. You can create one by calling Argument() on a dual space or by calling Coargument() directly on a dual space.

# Primal unknown (TrialFunction) becomes a Cofunction when assembled
a = TrialFunction(V)
f_1 = a * dx
cf = assemble(f_1)  # type Cofunction

# Dual unknown (Coargument)
u = Argument(V.dual(), 2)  # type Coargument
w = Coargument(V.dual(), 3)  # type Coargument

The MassInvPC preconditioner approximates the inverse of the Schur complement using a pressure mass inverse. This is effective for constant viscosity problems.

Key behaviors:

• Viscosity Weighting: For variable but low-contrast viscosity, you can pass a dictionary containing the viscosity (e.g., "mu") into the solve call. If not provided, it defaults to 1.0.
• Matrix Assembly: By default, the mass matrix is assembled. However, you can use an unassembled mass matrix by setting "fieldsplit_1_Mp_mat_type": "matfree" within the preconditioner parameters.
• Unassembled Mass Matrix: If using a matrix-free mass matrix, standard preconditioners like ilu cannot be used. Instead, you must use a Krylov solver (e.g., cg) with no preconditioner (pc_type: none).

# Using an unassembled mass matrix with CG
parameters["fieldsplit_1_Mp_mat_type"] = "matfree"
parameters["fieldsplit_1_Mp_pc_type"] = "ksp"
parameters["fieldsplit_1_Mp_ksp_ksp_type"] = "cg"
parameters["fieldsplit_1_Mp_ksp_pc_type"] = "none"

solve(a == L, up, bcs=bcs, nullspace=nullspace, solver_parameters=parameters)

If you are working with operators in a MixedFunctionSpace that are singular, you must specify the null space for each diagonal block separately. This is done using a MixedVectorSpaceBasis.

How to use MixedVectorSpaceBasis:

1. Constructor Arguments: The first argument must be the MixedFunctionSpace you are building the basis for. The second argument is a list of VectorSpaceBasis objects (one for each diagonal block).
2. Partial Null Spaces: You do not need to provide a null space for every block. For blocks where the null space is not a concern, pass an indexed function space using the .sub(i) method of the mixed space at the appropriate position in the list.

Example: For a mixed space $W = V \times Q$ where only $V$ has a constant null space, you would pass the constant basis for $V$ and the sub-space of $W$ for $Q$.

V = ...
Q = ...
W = V*Q
v_basis = VectorSpaceBasis(constant=True)

# Define null space for the first block (V) and use the sub-space for the second (Q)
nullspace = MixedVectorSpaceBasis(W, [v_basis, W.sub(1)])