Rayleigh-Benard Convection

This problem involves a variable-temperature incompressible fluid. Variations in the fluid temperature are assumed to affect the momentum balance through a buoyant term (the Boussinesq approximation), leading to a Navier-Stokes equation with a nonlinear coupling to a convection-diffusion equation for temperature.

We will set up the problem using Taylor-Hood elements for the Navier-Stokes part, and piecewise linear elements for the temperature.

from firedrake import *

N = 128

M = UnitSquareMesh(N, N)

V = VectorFunctionSpace(M, "CG", 2)
W = FunctionSpace(M, "CG", 1)
Q = FunctionSpace(M, "CG", 1)
Z = V * W * Q

upT = Function(Z)
u, p, T = split(upT)
v, q, S = TestFunctions(Z)

Two key physical parameters are the Rayleigh number (Ra), which measures the ratio of energy from buoyant forces to viscous dissipation and heat conduction and the Prandtl number (Pr), which measures the ratio of viscosity to heat conduction.

Ra = Constant(200.0)
Pr = Constant(6.8)

Along with gravity, which points down.

g = Constant((0, -1))

F = (
    inner(grad(u), grad(v))*dx
    + inner(dot(grad(u), u), v)*dx
    - inner(p, div(v))*dx
    - (Ra/Pr)*inner(T*g, v)*dx
    + inner(div(u), q)*dx
    + inner(dot(grad(T), u), S)*dx
    + 1/Pr * inner(grad(T), grad(S))*dx
)

There are two common versions of this problem. In one case, heat is applied from bottom to top so that the temperature gradient is enforced parallel to the gravitation. In this case, the temperature difference is applied horizontally, perpendicular to gravity. It tends to make prettier pictures for low Rayleigh numbers, but also tends to take more Newton iterations since the coupling terms in the Jacobian are a bit stronger. Switching to the first case would be a simple change of bits of the boundary associated with the second and third boundary conditions below:

bcs = [
    DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3, 4)),
    DirichletBC(Z.sub(2), Constant(1.0), (1,)),
    DirichletBC(Z.sub(2), Constant(0.0), (2,))
]

Like Navier-Stokes, the pressure is only defined up to a constant.:

nullspace = MixedVectorSpaceBasis(
    Z, [Z.sub(0), VectorSpaceBasis(constant=True), Z.sub(2)])

First off, we’ll solve the full system using a direct solver. As previously, we use MUMPS, so wrap the solve in try/except to avoid errors if it is not available.

from firedrake.petsc import PETSc

try:
   solve(F == 0, upT, bcs=bcs, nullspace=nullspace,
         solver_parameters={"mat_type": "aij",
                            "snes_monitor": None,
                            "ksp_type": "gmres",
                            "pc_type": "lu",
                            "pc_factor_mat_solver_type": "mumps"})
except PETSc.Error as e:
    if e.ierr == 92:
        warning("MUMPS not installed, skipping direct solve")
    else:
        raise e

For our next trick, we will use a fieldsplit preconditioner. This time, rather than using a Schur complement, we will use a multiplicative type (effectively block Gauss-Seidel). As ever, this has more options, so we’ll use a parameters dictionary. We use matrix-free actions for the coupled operator, and solve the linearised system with GMRES preconditioned with a multiplicative fieldsplit.

parameters = {"mat_type": "matfree",
              "snes_monitor": None,
              "ksp_type": "gmres",
              "pc_type": "fieldsplit",
              "pc_fieldsplit_type": "multiplicative",

We want to split the Navier-Stokes part off from the temperature variable.

"pc_fieldsplit_0_fields": "0,1",
"pc_fieldsplit_1_fields": "2",

We’ll invert the Navier-Stokes block with MUMPS:

"fieldsplit_0_ksp_type": "preonly",
"fieldsplit_0_pc_type": "python",
"fieldsplit_0_pc_python_type": "firedrake.AssembledPC",
"fieldsplit_0_assembled_pc_type": "lu",
"fieldsplit_0_assembled_pc_factor_mat_solver_type": "mumps",

the temperature block will also be inverted directly, but with plain LU.:

"fieldsplit_1_ksp_type": "preonly",
"fieldsplit_1_pc_type": "python",
"fieldsplit_1_pc_python_type": "firedrake.AssembledPC",
"fieldsplit_1_assembled_pc_type": "lu"}

Now for the solve.

upT.assign(0)
try:
    solve(F == 0, upT, bcs=bcs, nullspace=nullspace,
          solver_parameters=parameters)
except PETSc.Error as e:
    if e.ierr == 92:
        warning("MUMPS not installed, skipping assembled fieldsplit solve")
    else:
        raise e

Finally, we’ll demonstrate recursive fieldsplitting. We’ll use the same multiplicative fieldsplit preconditioner for the velocity-pressure and temperature blocks, but we’ll precondition the Navier-Stokes part with PCDPC using a lower Schur complement factorisation, and approximately invert the temperature block using algebraic multigrid. There are lots of parameters here, so let’s run through them. Since there are many options here, in particular for the nested subsolves, we specify options using nested, rather than flat, dictionaries. The solver parameters dictionary can either be a flat dictionary of key-value pairs, where both the keys and the values are strings, or it can be nested. In the latter case, the value should be a dictionary, of options and the key is prepended to all keys in the dictionary before passing to the solver.

parameters = {"mat_type": "matfree",
              "snes_monitor": None,

We’ll use inexact GMRES solves to invert the Navier-Stokes block, so the preconditioner as a whole is not stationary, hence we need flexible GMRES.

"ksp_type": "fgmres",
"ksp_gmres_modifiedgramschmidt": True,
"pc_type": "fieldsplit",
"pc_fieldsplit_type": "multiplicative",

Again we split off Navier-Stokes from the temperature block

"pc_fieldsplit_0_fields": "0,1",
"pc_fieldsplit_1_fields": "2",

which we solve inexactly using preconditioned GMRES.

"fieldsplit_0": {
    "ksp_type": "gmres",
    "ksp_gmres_modifiedgramschmidt": True,
    "ksp_rtol": 1e-2,
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "schur",
    "pc_fieldsplit_schur_fact_type": "lower",

Invert the velocity block with a single V-cycle of algebraic multigrid:

"fieldsplit_0": {
    "ksp_type": "preonly",
    "pc_type": "python",
    "pc_python_type": "firedrake.AssembledPC",
    "assembled_pc_type": "hypre"
},

and approximate the Schur complement inverse with PCD.

"fieldsplit_1": {
     "ksp_type": "preonly",
     "pc_type": "python",
     "pc_python_type": "firedrake.PCDPC",

We need to configure the pressure mass and Poisson solves, along with how to apply the convection-diffusion operator. For the latter, we will use an assembled operator this time round.

         "pcd_Mp_ksp_type": "preonly",
         "pcd_Mp_pc_type": "ilu",
         "pcd_Kp_ksp_type": "preonly",
         "pcd_Kp_pc_type": "hypre",
         "pcd_Fp_mat_type": "aij"
    }
},

Now for the temperature block, we use a moderately coarse tolerance for algebraic multigrid preconditioned GMRES.

     "fieldsplit_1": {
          "ksp_type": "gmres",
          "ksp_rtol": "1e-4",
          "pc_type": "python",
          "pc_python_type": "firedrake.AssembledPC",
          "assembled_pc_type": "hypre"
     }
}

And we’re done with all the options. All that’s left is to solve the problem. Recall that the PCD preconditioner needs to know where the velocity space lives in the velocity-pressure block, which we provide through the application context argument. It also needs to know the Reynolds number, which defaults to 1.0, which happens to work for our problem setup. We haven’t added the Rayleigh or Prandtl numbers to the dictionary since our known preconditioners don’t actually require them, although doing so would be quite easy.:

appctx = {"velocity_space": 0}
upT.assign(0)

solve(F == 0, upT, bcs=bcs, nullspace=nullspace,
      solver_parameters=parameters, appctx=appctx)

Finally, we’ll output the results for visualisation.

u, p, T = upT.split()
u.rename("Velocity")
p.rename("Pressure")
T.rename("Temperature")

File("benard.pvd").write(u, p, T)

A runnable python script implementing this demo file is available here.