Notebook 6: Time-dependence

Convection solution in an annulus with fixed inner, outer temperatures and free slip boundaries For details, see the notebook code.

We’ll look at a convection problem which couples Stokes Flow with time-dependent advection/diffusion.

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.

import underworld3 as uw
import numpy as np
import sympy

res = 10
r_o = 1.0
r_i = 0.55

rayleigh_number = 3.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(10000 * Gamma_N.dot(v_soln.sym) *  Gamma_N, "Upper")

if r_i != 0.0:
    stokes.add_natural_bc(10000 * 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")

Underworld expressions

Note that

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

underworld3.function.expressions.UWexpression

underworld3.function.expressions.UWexpression

\(\displaystyle 1\)

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)

with meshball.access(t_soln):
    t_soln.data[:,0] = 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 = 50
timestep = 0
elapsed_time = 0.0

adv_diff.solve(timestep=0.1)

# Null space ?

for step in range(0, max_steps):

    stokes.solve(zero_init_guess=False)
    delta_t = adv_diff.estimate_dt() 
    adv_diff.solve(timestep=delta_t)

    timestep += 1
    elapsed_time += delta_t

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

Timestep: 5, time 0.007655378472205266
Timestep: 10, time 0.01387213547299994
Timestep: 15, time 0.017197945909787137
Timestep: 20, time 0.019855789890665786
Timestep: 25, time 0.02250450167934171
Timestep: 30, time 0.02546333347384266
Timestep: 35, time 0.028720014807997008
Timestep: 40, time 0.032113891988035236
Timestep: 45, time 0.035501720967277196
Timestep: 50, time 0.038852418750145654

# 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)

    
    skip = 1
    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=500, 
    )
   

    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.3, show_scalar_bar=False, cmap="Greens", render_lines_as_tubes=False)

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

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

Interactive Image: Convection model output

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 it’s wise to check anyway.

    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)