Skip to main content

session4

Session4: Building a convection code

Lecture slides for this session can are here.

4.1. Go through your advection-diffusion code, and make it non-dimensional:

  • Your box dimensions become 1-by-1 now (so the non-dimensional box height becomes h′=1h′=1)
  • The non-dimensional temperature will now vary between 0 and 1 (so T′m=1Tm′=1).
  • κκ vanishes from the advection-diffusion equation (so its non-dimensional value becomes κ′=1κ′=1)
  • This will affect the non-dimensional time step and velocity as well
  • Because hh, TmTm, and κκ are all reduced to 11 in this non-dimensional code, they can be removed entirely, which will simplify/shorten the code.
  • Test your code: when plotting the (re-dimensionalised) temperature through (re-dimensionalised) time, the results should be identical to those from your dimensional code.

Next, we will convert your code into a thermal convection code. Your advection-diffusion code will be the main routine that will call a Stokes solver subfunction that calculates the velocity and pressure field. I suggest to use the following recipe:

4.2. Copy the following Stokes solver code into a separate python script file, and examine and run it to see what it does.

# Stokes2D.py
#
# Subitop modelling course, Edinburgh 2017
# Jeroen van Hunen, March 2017
# purpose: calculates 2D Stokes flow
# method:
# Version: 2D, isoviscous, equidistant grid

import numpy as np
import pylab as plt
import scipy
import scipy.sparse            # Needed for some older Spyder versions
import scipy.sparse.linalg     # Needed for some older Spyder versions

# Subfunctions: 
def idp(ix,iz,nz):
    # purpose: find index value for the p component in 2-D Stokes matrix
    fout = 3*((ix-1)*nz + iz) - 3
    return fout
    
def idvx(ix,iz,nz):
    # purpose: find index value for the vx component in 2-D Stokes matrix
    fout = 3*((ix-1)*nz + iz) - 2
    return fout
    
def idvz(ix,iz,nz):
    # purpose: find index value for the vz component in 2-D Stokes matrix
    fout = 3*((ix-1)*nz + iz) - 1
    return fout
    
def getpvxvz(sol,xx):
    # purpose: extract p, vx, and vz from sol vector from matrix inversion
    # sol contains [p1vx1vz1p2vx2vz2p3vx3vz3....]
    [nz,nx]=np.shape(xx)
    pp=sol[::3]          # Extract every 3rd position as p from sol
    vx=sol[1::3]         # ... and vx
    vz=sol[2::3]         # ... and vz
    p=np.reshape(pp,(nz,nx), order='F')     # shape solvp into nx-by-nz mesh
    vx=np.reshape(vx,(nz,nx), order='F')    # idem for vx
    vz=np.reshape(vz,(nz,nx), order='F')    # idem for vz
    p=p[1:,1:]           # remove first row and column: ghost points
    meanp=np.mean(p)     # subtract mean to make average p=0
    p=p-meanp               # 

    vx=vx[0:-1,0:]       # remove ghost points from vx
    vz=vz[0:,0:-1]       # remove ghost points from vz
    return [p, vx, vz]
    
def preppvxvzplot(pp,vx,vz,xx,islip):
    # Purpose: interpolate p, vx, and vz to the base points
    #          for plotting
    # Method:  p, vx, and vz are each defined at their own location 
    #          on the staggered grid, which makes plotting difficult
    #          vx on staggered grid is vertically between base points
    #            expand vx array with top and bottom row, note that this
    #            done differently for free-slip and no-slip.
    #            Then interpolate vertically to midpoints, which are the 
    #            base points 
    #          vz on staggered grid is horizontally between base points
    #            done as for vx, but horizontally, not vertically
    #          p on staggered grid is diagonally between base points
    #            done as vx and vz, but both horizontally and vertically
    #            This implies both hor and vert interpolation
    #          p, vx, and vz now all have nx by nz points
    # Arguments:
    #     pp = raw pressure field
    #     vx is raw x-velocity field
    #     vz is raw z-velocity field
    #     islip is slip type on bnds: 1=free-slip, -1=no-slip
    
    # vertically interpolate vx to base points: 
    vxplot=np.zeros(np.shape(xx))
    vxplot[1:-1,:]=0.5*(vx[0:-1,:]+vx[1:,:])
    # and extrapolate vx from 1st/last row to top/bottom bnd.
    vxplot[0,:]=vx[0,:]
    vxplot[-1,:]=vx[-1,:]
    
    # horizontally interpolate vz to base points: 
    vzplot=np.zeros(np.shape(xx))
    vzplot[:,1:-1]=0.5*(vz[:,0:-1]+vz[:,1:])
    # and extrapolate vz from 1st/last column to lef/right bnd.
    vzplot[:,0]=vz[:,0]
    vzplot[:,-1]=vz[:,-1]
    
    # interpolate p to base points for plotting purposes:
    pplot=np.zeros(np.shape(xx))
    # bilinear interpolation of p-points to all internal points:
    pplot[1:-1,1:-1]=0.25*(pp[0:-1,0:-1]+pp[1:,0:-1]+pp[0:-1,1:]+pp[1:,1:])
    pplot[0,1:-1]=0.5*(pp[0,0:-1]+pp[0,1:])
    # Boundary points only have two nearest p-points:
    pplot[-1,1:-1]=0.5*(pp[-1,0:-1]+pp[-1,1:])
    pplot[1:-1,0]=0.5*(pp[0:-1,0]+pp[1:,0])
    pplot[1:-1,-1]=0.5*(pp[0:-1,-1]+pp[1:,-1])
    # Corner points only have one associated internal point:
    pplot[0,0]=pp[0,0]
    pplot[-1,0]=pp[-1,0]
    pplot[0,-1]=pp[0,-1]
    pplot[-1,-1]=pp[-1,-1]
    
    return [pplot,vxplot,vzplot]

# Main routine: 
nx   = 21 
nz   = 21
islip = 1  # 1=free-slip -1=no-slip

nxz=3*nx*nz  # total nr of unknowns (nx * nz * (vx+vz+p))
x        = np.linspace(0,1,nx)
z        = np.linspace(0,1,nz)
dx   = x[1]-x[0]
dz   = z[1]-z[0]
[xx,zz]=np.meshgrid(x,z)

A    = scipy.sparse.lil_matrix((nxz,nxz))
#A.setdiag(1)                    # Does not work on older Spyder versions
for i in range(0,nxz):
    A[i,i]=1.
rhs  = np.zeros(nxz)             # create rhs (buoyancy force) vector
drho = np.zeros(np.shape(xx))    # create density distribution matrix
drho[:,0:int(nx/2)]=-0.5         # Symmetric buoyancy for both odd and even 
drho[:,int(nx/2)+1:]=0.5         #   number grids: if odd-> middle column rho=0
Ra = 1e4                         # Rayleigh number

# Fill in info in matrix for Stokes_z for internal points & left/right bnd. points:
# Note: 1) other points (top/bottom bnd.pnts and ghstpnts have vz=0, which is default
#       2) Index counters: ix=1 to nx-1, and iz=1 to nz (unusual for Python)
for iz in range (2,nz):          # iz=1 & iz=nz are default (i.e. vz=0) bc's: no calc needed
    for ix in range (1,nx):      # ix=nx is ghostpoint ix=1 & nx-1 are boundary, 
                                 #     but vz still needs calculating
        # calculate indices of all relevant grid points for vz and p:
        # for vz
        vc = idvz(ix,iz,nz)      # calculate matrix index for central vz point:
        if (ix>1):
            vl = idvz(ix-1,iz,nz)# idem, for left vx point
        if (ix<nx-1):
            vr = idvz(ix+1,iz,nz)# idem, for right vz point
        vt = idvz(ix,iz+1,nz)    # idem, for top vz point
        vb = idvz(ix,iz-1,nz)    # idem, for bottom vz point
        # for p:
        pt = idp(ix+1,iz+1,nz)   # idem, for left p point
        pb = idp(ix+1,iz,nz)     # idem, for right p point
        
        # fill in matrix components:
        irow = idvz(ix,iz,nz)
        A[irow,vc] = -2/dx**2-2/dz**2 # valid for internal points only
        if (ix>1):
            A[irow,vl] = 1/dx**2
        else:
            # free-slip add correction to central point
            A[irow,vc] = A[irow,vc] + islip*1/dx**2

        if (ix<nx-1):
            A[irow,vr] = 1/dx**2
        else:
            # free-slip add correction to central point
            A[irow,vc] = A[irow,vc] + islip*1/dx**2
        A[irow,vt] = 1/dz**2
        A[irow,vb] = 1/dz**2
        A[irow,pb] = 1/dz
        A[irow,pt] = -1/dz
        
        # rhs: Ra*drho'
        avdrho  = 0.5*(drho[iz-1,ix-1]+drho[iz-1,ix])
        rhs[irow] = avdrho*Ra

# Fill in info in matrix for Stokes_x for internal points & top/bottom bnd. points:
# Note: other points (left/right bnd.pnts and ghstpnts have vx=0, which is default
for ix in range (2,nx):          # ix=1 & nx are default (i.e. vx=0) bc's: no calc, needed
    for iz in range (1,nz):      # iz=nz are ghostpoints, iz=1&nz-1 are boundaries, 
                                 #     but vx still needs calculating there
        # calculate indices of all relevant grid points for vx and p:
        # for vx
        vc = idvx(ix,iz,nz)      # calculate matrix index for central vx point:
        vl = idvx(ix-1,iz,nz)    # idem, for left vx point
        vr = idvx(ix+1,iz,nz)    # idem, for right point
        if (iz<nz-1):
            vt = idvx(ix,iz+1,nz)# idem, for top vx point
        if (iz>1):
            vb = idvx(ix,iz-1,nz)# idem, for bottom vx point
        # for p:
        pl = idp(ix,iz+1,nz) # idem, for left p point
        pr = idp(ix+1,iz+1,nz)   # idem, for right p point
        
        # fill in matrix components:
        irow = idvx(ix,iz,nz)
        A[irow,vc] = -2/dx**2-2/dz**2 # valid for internal points only
        A[irow,vl] = 1/dx**2        
        A[irow,vr] = 1/dx**2         
        if (iz<nz-1):            # top bnd.point
            A[irow,vt] = 1/dz**2
        else:
            # free-slip add correction to central point
            A[irow,vc] = A[irow,vc] + islip*1/dz**2 
        if(iz>1):                # bottom bnd.point
            A[irow,vb] = 1/dz**2
        else:
            # free-slip add correction to central point
            A[irow,vc] = A[irow,vc] + islip*1/dz**2 # free-slip
        A[irow,pl] = 1/dx
        A[irow,pr] = -1/dx
        # all rhs components here are 0: is default

# Fill in info in matrix for continuity eqn for all pressure points:
for ix in range (2,nx+1):       # pressure point ix=1 is a ghostpoint
    for iz in range (2,nz+1):   # pressure point iz=1 is a ghostpoint
        irow=idp(ix,iz,nz)
        vxl=idvx(ix-1,iz-1,nz)
        vxr=idvx(ix,iz-1,nz)
        vzb=idvz(ix-1,iz-1,nz)
        vzt=idvz(ix-1,iz,nz)
        A[irow,vxl]=-1/dx
        A[irow,vxr]=1/dx
        A[irow,vzb]=-1/dz
        A[irow,vzt]=1/dz
        A[irow,irow]=0

# fix p=0 at one point: lowerleft corner:
irow=idp(2,2,nz)               
A[irow,irow]=1
rhs[irow]=0

# Solve system:
sol=scipy.sparse.linalg.spsolve(A,rhs)

# extrac p, vx, and vz solutions from solution vector:
[pp,vx,vz] = getpvxvz(sol,xx)
vxmax=vx.max()
vzmax=vz.max()
print(vxmax)
print(vzmax)

# preparing p, vx, and vz for plotting:
[pplot,vxplot,vzplot]=preppvxvzplot(pp,vx,vz,xx,islip)

# plot solution:
plt.clf()
plt.figure(1)
plt.imshow(pplot, 
           extent=(0,1,0,1),
           #clim=(0,Tm),
           interpolation='bilinear', 
           cmap='jet')
plt.quiver(xx, zz, vxplot, vzplot, units='width')
plt.title('Stokes flow')
plt.pause(0.00005)

4.3. Convert the main routine Stokes2D into a subfunction Stokes2Dfunc that can be called from your advection-diffusion solver:

  • Instead of defining Ra, and the mesh setup, introduce this information through function arguments.
  • The nondimensional density drho will be replaced by the nondimensional temperature, and will also be provided through a function argument.
  • The routine calculates pressure and both velocity components on staggered grids and therefore not located at the same coordinates. But the sub-function preppvxvzplot interpolates these quantities to pplot,
    vxplot, and vzplot(which are defined at the xx and zz grid locations, and therefore suitable for plotting). Transfer pplot, vxplot, and vzplot back to the main routine.
  • Test your Stokes solver as a sub-function by building a short main program that a) sets up the grid and Ra,
    calls the Stokes solver, and plots the solution. For the plotting, move the plotting routine in the Stokes solver and put it in the main routine.
  • A model solution is available here.

4.4. Replace the fixed solid rotation velocity field with the velocity field that is calculated in the Stokes routine:

  • Remove the definition of the solid rotation velocity part
  • Add the Stokes solver at the beginning of the time loop
  • Since every timestep has a different velocity field, the critical advection timestep (i.e. the Courant timestep) will be different too, and needs to be recalculated every timestep.
  • Plot the velocity arrows on top of a temperature colour plot every timestep, using the plotting setup from the Stokes2D.py routine as an example.

4.5. Calculate the resulting (dimensional) topography from the convection results: topo=−σzzρgtopo=−σzzρg:

  • First calculate the non-dimensional normal surface stress −σzz−σzz using the lecture notes
  • Dimensionalize the solution, using the following parameters: κ=10−6 m2/s, η=1027/Raη=1027/Ra Pa s, and h=1000h=1000 km.
  • Finally calculate the dimensional topography using a density contrast between the solid Earth and the surface of 4000 kg/m3, and g=10 m/s2.
  • The topography variation is the most useful information in here, since the absolute value of this topography has little meaning. Therefore subtract the average from the topography value.

A model solution is available here.

4.6. Calculate an average (‘root-mean-square’) non-dimensional velocity field vrms=∫∫v2x+v2zdzdx‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√vrms=∫∫vx2+vz2dzdx. If you want to use it, a routine that simply sums up nodal point values (not very accurate and not valid for non-uniform grids) for this is provided here. It is advised to calculate the vrmsvrms on the original (staggered grid) velocity solution, and not on the solution tht is interpolated to the xx grid.

4.7. Test you code against the published benchmark of Blankenbach et al. (1989). Reported best values are

Ravrmstopomin(m)topomax (m)
10⁴42.86-2903.2254.
10⁵193.21-2004.1461.
10⁶833.99-1284.932.