# 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 :class:`~.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 :ref:`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 `__.