Notebook 6: Rayleigh-Bénard Convection (time-stepping example)

Run this notebook on Binder

We’ll look at a convection problem which couples Stokes Flow with time-dependent advection/diffusion to give simple Rayleigh-Bénard convection model.

\[ -\nabla \cdot \left[ \frac{\eta}{2}\left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) - p \mathbf{I} \right] = -\rho_0 \alpha T \mathbf{g} \]
\[ \nabla \cdot \mathbf{u} = 0 \]

\(\eta\) is viscosity, \(p\) is pressure, \(\rho_0\) is a reference density, \(\alpha\) is thermal expansivity, and \(T\) is the temperature. Here we explicitly express density variations in terms of temperature variations.

Thermal evolution is given by

\[ \frac{\partial T}{\partial t} - \mathbf{u}\cdot\nabla T = \kappa \nabla^2 T \]
where the velocity, \(\mathbf{u}\) is the result of the Stokes flow calculation. \(\kappa\) is the thermal diffusivity (compare this with Notebook 4).

The starting point is our previous notebook where we solved for Stokes flow in a cylindrical annulus geometry. We then add an advection-diffusion solver to evolve temperature. The Stokes buoyancy force is proportional to the temperature anomaly, and the velocity solution is fed back into the temperature advection term. The timestepping loop is written by hand because usually you will want to do some analysis or output some checkpoints.

To read more about the applications of simple mantle convection models like this one, see (for example) Schubert et al, 2001.

#|  echo: false  # Hide in html version

# This is required to fix pyvista
# (visualisation) crashes in interactive notebooks (including on binder)

import nest_asyncio

nest_asyncio.apply()

#| output: false # Suppress warnings in html version

import numpy as np
import sympy
import underworld3 as uw

uw.__file__

res = 5
r_o = 1.0
r_i = 0.55

norm_r = (r_o - r_i)
r_o /= norm_r
r_i /= norm_r


rayleigh_number = 1.0e4

meshball = uw.meshing.Annulus(
    radiusOuter=r_o,
    radiusInner=r_i,
    cellSize=1 / res,
    qdegree=3,
)

# Coordinate directions etc
x, y = meshball.CoordinateSystem.X
r, th = meshball.CoordinateSystem.xR
unit_rvec = meshball.CoordinateSystem.unit_e_0

# Orientation of surface normals
Gamma_N = unit_rvec

# Mesh variables for the unknowns

v_soln = uw.discretisation.MeshVariable("V0", meshball, 2, degree=2, varsymbol=r"{v_0}")
p_soln = uw.discretisation.MeshVariable("p", meshball, 1, degree=1, continuous=True)
t_soln = uw.discretisation.MeshVariable("T", meshball, 1, degree=3, continuous=True)

Create linked solvers

We create the Stokes solver as we did in the previous notebook. The buoyancy force is proportional to the temperature anomaly (t_soln). Solvers can either be provided with unknowns as pre-defined meshVariables, or they will define their own. When solvers are coupled, explicitly defining unknowns makes everything clearer.

The advection-diffusion solver evolved t_soln using the Stokes velocity v_soln in the fluid-transport term.

Curved, free-slip boundaries

In the annulus, a free slip boundary corresponds to zero radial velocity. However, in this mesh, \(v_r\) is not one of the unknowns (\(\mathbf{v} = (v_x, v_y)\)). We apply a non linear boundary condition that penalises \(v_r\) on the boundary as discussed previously in Example 5.

stokes = uw.systems.Stokes(
    meshball,
    velocityField=v_soln,
    pressureField=p_soln,
)

stokes.bodyforce = rayleigh_number * t_soln.sym * unit_rvec

stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
stokes.tolerance = 1.0e-3

stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"

stokes.add_natural_bc(1000000 * Gamma_N.dot(v_soln.sym) * Gamma_N, "Upper")

if r_i != 0.0:
    stokes.add_natural_bc(1000000 * Gamma_N.dot(v_soln.sym) * Gamma_N, "Lower")

# Create solver for the energy equation (Advection-Diffusion of temperature)

adv_diff = uw.systems.AdvDiffusion(
    meshball,
    u_Field=t_soln,
    V_fn=v_soln,
    order=2,
    verbose=False,
)

adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel
adv_diff.constitutive_model.Parameters.diffusivity = 1

## Boundary conditions for this solver

adv_diff.add_dirichlet_bc(+1.0, "Lower")
adv_diff.add_dirichlet_bc(-0.0, "Upper")

# adv_diff.petsc_options.setValue("snes_monitor", None)
# adv_diff.petsc_options.setValue("ksp_monitor", None)

uw.systems.AdvDiffusion.view()

uw.systems.SemiLagragian_DDt.view()

Underworld expressions

Note that the parameters in the consititutive models are defined as uw.expressions because

  1. We generally would like to display them symbolically in equations rather than see their floating point value

  2. They can represent complex expressions that “blow up” the symbolic forms when we inspect the equations

  3. We would like to be able to substitute new values in the equations and use the lazy evaluation of expressions to avoid having to redefine expressions in multiple locations where these symbols appear

  4. For future reference, in optimisation problems, we differentiate expressions with respect to their parameters and want the results also to be evaluated in a lazy fashion.

display(type(stokes.constitutive_model.Parameters.shear_viscosity_0))
display(type(adv_diff.constitutive_model.Parameters.diffusivity))

stokes.constitutive_model.Parameters.shear_viscosity_0.sym

Initial condition

We need to set an initial condition for the temperature field as the coupled system is an initial value problem. Choose whatever works but remember that the boundary conditions will over-rule values you set on the lower and upper boundaries.

# Initial temperature

init_t = 0.1 * sympy.sin(3 * th) * sympy.cos(np.pi * (r - r_i) / (r_o - r_i)) + (
    r_o - r
) / (r_o - r_i)

t_soln.array[...] = uw.function.evaluate(init_t, t_soln.coords)

Initial velocity solve

The first solve allows us to determine the magnitude of the velocity field and is useful to keep separated to check convergence rates etc.

For non-linear problems, we usually need an initial guess using a reasonably close linear problem.

zero_init_guess is used to reset any information in the vector of unknowns (i.e. do not use any initial information if zero_init_guess==True).

stokes.solve(zero_init_guess=True)

# Keep the initialisation separate
# so we can run the loop again in a notebook

max_steps = 25
timestep = 0
elapsed_time = 0.0

adv_diff.view(class_documentation=True)

# init_t = 0.1 * sympy.sin(3 * th) * sympy.cos(np.pi * (r - r_i) / (r_o - r_i)) + (
#     r_o - r
# ) / (r_o - r_i)


# t_soln.array[...] = uw.function.evaluate(
#         init_t + 0.0001 * t_soln.sym[0], t_soln.coords
#     )

adv_diff.estimate_dt().squeeze()

adv_diff.solve(timestep=0.01, zero_init_guess=True)

# Null space ?

for step in range(0, max_steps):

    stokes.solve(zero_init_guess=True)
    delta_t = 3 * adv_diff.estimate_dt()
    adv_diff.solve(timestep=delta_t, zero_init_guess=False, verbose=False)

    timestep += 1
    elapsed_time += delta_t

    if timestep % 5 == 0:
        print(f"Timestep: {timestep}, time {elapsed_time}")

# visualise it


if uw.mpi.size == 1:
    import pyvista as pv
    import underworld3.visualisation as vis

    pvmesh = vis.mesh_to_pv_mesh(meshball)
    pvmesh.point_data["P"] = vis.scalar_fn_to_pv_points(pvmesh, p_soln.sym)
    pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym)
    pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln.sym)

    pvmesh_t = vis.meshVariable_to_pv_mesh_object(t_soln)
    pvmesh_t.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln.sym)

    pvmesh_v = vis.meshVariable_to_pv_mesh_object(v_soln)
    pvmesh_v.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh_v, v_soln.sym)


    skip = 10
    points = np.zeros((meshball._centroids[::skip].shape[0], 3))
    points[:, 0] = meshball._centroids[::skip, 0]
    points[:, 1] = meshball._centroids[::skip, 1]
    point_cloud = pv.PolyData(points)

    pvstream = pvmesh.streamlines_from_source(
        point_cloud,
        vectors="V",
        integration_direction="both",
        integrator_type=45,
        surface_streamlines=True,
        initial_step_length=0.01,
        # max_time=1.0,
        max_steps=100,
    )

    pl = pv.Plotter(window_size=(750, 750))

    pl.add_mesh(
        pvmesh_t,
        cmap="RdBu_r",
        edge_color="Grey",
        edge_opacity=0.33,
        scalars="T",
        show_edges=True,
        use_transparency=False,
        opacity=1.0,
        show_scalar_bar=False,
    )

    pl.add_mesh(
        pvstream,
        opacity=0.5,
        show_scalar_bar=False,
        cmap="Greens",
        render_lines_as_tubes=False,
        line_width=1,
        
    )

    pl.add_arrows(pvmesh_v.points, pvmesh_v.point_data["V"], mag=1e-3 )


    pl.export_html("html5/annulus_convection_plot.html")
    # pl.show(cpos="xy", jupyter_backend="trame")

#| fig-cap: "Interactive Image: Convection model output"
from IPython.display import IFrame

IFrame(src="html5/annulus_convection_plot.html", width=500, height=400)

Exercise - Null space

Based on our previous notebook, can you see how to calculate and (if necessary) remove rigid-body the rotation null-space from the solution ?

The use of a coarse-level singular-value decomposition for the velocity solver should help, in this case, but sometimes you can see that there is a rigid body rotation (look at the streamlines). It’s wise to check and quantify the presence of the null space.

    stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"

Exercise - Heat flux

Could you calculate the radial heat flux field ? Its surface average value plotted against time tells you if you have reached a steady state.

Hint:

\[ Q_\textrm{surf} = \nabla T \cdot \hat{r} + T (\mathbf{v} \cdot \hat{r} ) \]
    Q_surf = -meshball.vector.gradient(t_soln.sym).dot(unit_rvec) +\
                    t_soln.sym[0] * v_soln.sym.dot(unit_rvec)

References

Schubert, G., Turcotte, D. L., & Olson, P. (2001). Mantle Convection in the Earth and Planets (1st ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511612879