The Stokes equations

\begin{align*} -\Delta \mathbf{u} + \nabla p &= \mathbf{f} &&\text{in } \Omega,\\ \nabla \cdot \mathbf{u} &= 0 &&\text{in } \Omega,\\ \mathbf{u} &= \mathbf{0}&&\text{on } \partial\Omega. \end{align*}

In this tutorial you will learn how to:

• Create manufactured solutions with UFL
• Use block-preconditioners

Defining a manufactured solution

We will use a known analytical solution to the Stokes equations in this tutorial. We define the exact velocity and pressure as the following:

def u_ex(x):
    sinx = sin(pi * x[0])
    siny = sin(pi * x[1])
    cosx = cos(pi * x[0])
    cosy = cos(pi * x[1])
    c_factor = 2 * pi * sinx * siny
    return c_factor * as_vector((cosy * sinx, - cosx * siny))

def p_ex(x):
    return sin(2 * pi * x[0]) * sin(2 * pi * x[1])

Here, the input to each function is the coordinates (x,y) of the problem. These will be defined by using x = ufl.SpatialCoordinate(domain).

We use the strong formulation of the PDE to compute the source function $\mathbf{f}$ using UFL operators

def source(x):
    u, p = u_ex(x), p_ex(x)
    return - div(grad(u)) + grad(p)

Defining the variational form

We will solve the PDE by creating a set of variational forms, one for each component of the problem. This leads to a blocked discrete system:

$$\begin{align} A w &= b,\\ \begin{pmatrix} A_{\mathbf{u},\mathbf{u}} & A_{\mathbf{u},p} \\ A_{p,\mathbf{u}} & 0 \end{pmatrix} \begin{pmatrix} u\\ p \end{pmatrix} &= \begin{pmatrix}\mathbf{f}\\ 0 \end{pmatrix} \end{align}$$

def create_bilinear_form(V, Q):
    u, p = TrialFunction(V), TrialFunction(Q)
    v, q = TestFunction(V), TestFunction(Q)
    a_uu = inner(grad(u), grad(v)) * dx
    a_up = inner(p, div(v)) * dx
    a_pu = inner(div(u), q) * dx
    return fem.form([[a_uu, a_up], [a_pu, None]])

def create_linear_form(V, Q):
    v, q = TestFunction(V), TestFunction(Q)
    domain = V.mesh
    x = SpatialCoordinate(domain)
    f = source(x)
    return fem.form([inner(f, v) * dx,
                     inner(fem.Constant(domain, 0.), q) * dx])

Create a block preconditioner

We create a nested matrix P to use as the preconditioner. The top-left block of P is the top-left block of A. The bottom-right diagonal entry is a mass matrix.

def create_preconditioner(Q, a, bcs):
    p, q = TrialFunction(Q), TestFunction(Q)
    a_p11 = fem.form(inner(p, q) * dx)
    a_p = fem.form([[a[0][0], None],
                    [None, a_p11]])
    P = fem.petsc.assemble_matrix_nest(a_p, bcs)
    P.assemble()
    return P

Assemble the nested system

def assemble_system(lhs_form, rhs_form, bcs):
    A = fem.petsc.assemble_matrix_nest(lhs_form, bcs=bcs)
    A.assemble()

    b = fem.petsc.assemble_vector_nest(rhs_form)
    fem.petsc.apply_lifting_nest(b, lhs_form, bcs=bcs)
    for b_sub in b.getNestSubVecs():
        b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD,
                          mode=PETSc.ScatterMode.REVERSE)
    spaces = fem.extract_function_spaces(rhs_form)
    bcs0 = fem.bcs_by_block(spaces, bcs)
    fem.petsc.set_bc_nest(b, bcs0)
    return A, b

PETSc Krylov Subspace solver

In legacy DOLFIN, convenience functions were provided to interact with linear algebra packages such as PETSc. In DOLFINx, we instead supply users with the appropriate data types, so that the user can access all of the features of the linear package rather than being constrained to functions in our wrapper. One can also leverage the detailed documentation of PETSc. For blocked systems, see: https://petsc.org/release/docs/manual/ksp/?highlight=matnest#solving-block-matrices

def create_block_solver(A, b, P, comm):
    ksp = PETSc.KSP().create(comm)
    ksp.setOperators(A, P)
    ksp.setType("minres")
    ksp.setTolerances(rtol=1e-9)
    ksp.getPC().setType("fieldsplit")
    ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

    nested_IS = P.getNestISs()
    ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]),
                                ("p", nested_IS[0][1]))

    # Set the preconditioners for each block
    ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
    ksp_u.setType("preonly")
    ksp_u.getPC().setType("gamg")
    ksp_p.setType("preonly")
    ksp_p.getPC().setType("jacobi")

    # Monitor the convergence of the KSP
    ksp.setFromOptions()
    return ksp

Solving the Stokes problem with a block-preconditioner

def solve_stokes(u_element, p_element, domain):
    V = fem.FunctionSpace(domain, u_element)
    Q = fem.FunctionSpace(domain, p_element)

    lhs_form = create_bilinear_form(V, Q)
    rhs_form = create_linear_form(V, Q)

    bcs = create_velocity_bc(V)
    nsp = create_nullspace(rhs_form)
    A, b = assemble_system(lhs_form, rhs_form, bcs)
    assert nsp.test(A)
    A.setNullSpace(nsp)

    P = create_preconditioner(Q, lhs_form, bcs)
    ksp = create_block_solver(A, b, P, domain.comm)

    u, p = fem.Function(V), fem.Function(Q)
    w = PETSc.Vec().createNest([u.vector, p.vector])
    ksp.solve(b, w)
    assert ksp.getConvergedReason() > 0
    u.x.scatter_forward()
    p.x.scatter_forward()
    return compute_errors(u, p)

We now use the Stokes solver we have defined to experiment with a range of element pairs that can be used. First, we define a function that takes a pair of elements as input and plots a graph showing the error as $h$ is decreased.

Crouzeix-Raviart

Alternatively, the same order convergence can be achieved using fewer degrees of freedom if a Crouzeix-Raviart {cite}crouzeix-raviart element is used for the velocity space.

element_u = VectorElement("CR", "triangle", 1)
element_p = FiniteElement("DG", "triangle", 0)
error_plot(element_u, element_p, 1)

Piecewise linear pressure space

When using a piecewise linear pressure space, we could again try using a velocity space one degree higher, but we would again observe that there is no convergence. In order to achieve convergence, we can augment the quadratic space with a cubic bubble function on the triangle {cite}crouzeix-falk.

enriched_element = (
    FiniteElement("Lagrange", "triangle", 2)
    + FiniteElement("Bubble", "triangle", 3))
element_u = VectorElement(enriched_element)
element_p = FiniteElement("DG", "triangle", 1)
error_plot(element_u, element_p, 2)

Taylor-Hood Element

Continuous Lagrange spaces for pressure and velocity.

element_u = VectorElement("CG", "triangle", 3)
element_p = FiniteElement("CG", "triangle", 2)
error_plot(element_u, element_p, 3)

When using a piecewise quadratic space, we can use a cubic velocity space augmented with quartic bubbles.

We have to define this velocity element as a custom element (it cannot be created as an enriched element, as the basis functions of degree 3 Lagrange and degree 4 bubbles are not linearly independent). More examples of how custom elements can be created can be found in the Basix documentation.

Defining the polynomial space

When creating a custom element, we must input the coefficients that define a basis of the set of polynomials that our element spans. In this example, we will represent the 9 functions in our space in terms of the 10 orthonormal polynomials of degree $\leqslant3$ on a quadrilateral, so we create a 9 by 10 matrix.

A polynomial $f$ on the triangle can be written as $$f(x,y)=\sum_i\left(\int_0^1\int_0^{1-y}f(x, y)q_i(x, y) \,\mathrm{d}x\,\mathrm{d}y\right)q_i(x,y),$$ where $q_0$, $q_1$, ... are the orthonormal polynomials. The entries of our coefficient matrix are these integrals.

wcoeffs = np.zeros((12, 15))
pts, wts = basix.make_quadrature(basix.CellType.triangle, 8)
poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.triangle, 4, pts)
x = pts[:, 0]
y = pts[:, 1]
for j, f in enumerate([
    1, x, y, x ** 2 * y, x * y ** 2, (1 - x - y) ** 2 * y, (1 - x - y) * y ** 2,
    x ** 2 * (1 - x - y), x * (1 - x - y) ** 2,
    x * y * (1 - x - y), x ** 2 * y * (1 - x - y), x * y ** 2 * (1 - x - y)
]):
    for i in range(15):
        wcoeffs[j, i] = sum(f * poly[i, :] * wts)

Interpolation

Next, we compute the points and matrices that define how functions can be interpolated into this space. For this element, all the DOFs are point evaluations, so we create lists of these points and (reshaped) identity matrices.

x = [[], [], [], []]
x[0].append(np.array([[0.0, 0.0]]))
x[0].append(np.array([[1.0, 0.0]]))
x[0].append(np.array([[0.0, 1.0]]))
x[1].append(np.array([[2 / 3, 1 / 3], [1 / 3, 2 / 3]]))
x[1].append(np.array([[0.0, 1 / 3], [0.0, 2 / 3]]))
x[1].append(np.array([[1 / 3, 0.0], [2 / 3, 0.0]]))
x[2].append(np.array([[1 / 4, 1 / 4], [1 / 2, 1 / 4], [1 / 4, 1 / 2]]))

M = [[], [], [], []]
for _ in range(3):
    M[0].append(np.array([[[[1.]]]]))
for _ in range(3):
    M[1].append(np.array([[[[1.], [0.]]], [[[0.], [1.]]]]))
M[2].append(np.array([[[[1.], [0.], [0.]]], [[[0.], [1.], [0.]]], [[[0.], [0.], [1.]]]]))

Creating the element

We now create the element by passing in the information we created above, as well as the cell type, value shape, number of derivatives used by the DOFs, map type, whether the element is discontinuous, the highest degree Lagrange space that is a subspace of the element, and the polynomial degree of the element.