Steady-State Subduction Zone Setup

Steady-State Subduction Zone Setup#

Authors: Kidus Teshome, Cian Wilson

Recalling our implementation strategy we are following a similar workflow to that seen repeatedly in the background examples.

  1. we will describe the subduction zone geometry and tesselate it into non-overlapping triangles to create a mesh

  2. we will declare function spaces for the temperature, wedge velocity and pressure, and slab velocity and pressure

  3. using these function spaces we will declare trial and test functions

  4. we will define Dirichlet boundary conditions at the boundaries as described in the introduction

  5. we will describe discrete weak forms for temperature and each of the coupled velocity-pressure systems that will be used to assemble the matrices (and vectors) to be solved

  6. we will set up matrices and solvers for the discrete systems of equations

  7. we will solve the matrix problems

In the previous four notebooks we have implemented basic functionality for all the steps. What remains is to implement specific cases that include the temperature weak forms and its coupled solution to the flow in the subduction zone. In the following notebooks we do that for the steady-state cases that allow us to test our method against the case 1 benchmark solution of Wilson & van Keken, 2023.

Preamble#

Let’s start by adding the path to the modules in the python folder to the system path (so we can find the our custom modules).

import sys, os
basedir = ''
if "__file__" in globals(): basedir = os.path.dirname(__file__)
sys.path.append(os.path.join(basedir, os.path.pardir, os.path.pardir, 'python'))

Let’s also load the module generated by the previous notebooks to get access to the parameters and functions defined there.

from sz_problems.sz_params import default_params, allsz_params
from sz_problems.sz_slab import create_slab
from sz_problems.sz_geometry import create_sz_geometry
from sz_problems.sz_problem import SubductionProblem

Then let’s load all the required modules at the beginning.

import geometry as geo
import utils
from mpi4py import MPI
import dolfinx as df
import dolfinx.fem.petsc
from petsc4py import PETSc
import numpy as np
import scipy as sp
import ufl
import basix.ufl as bu
import matplotlib.pyplot as pl
import copy
import pyvista as pv
import pathlib
output_folder = pathlib.Path(os.path.join(basedir, "output"))
output_folder.mkdir(exist_ok=True, parents=True)

SteadySubductionProblem class#

We continue building on the SubductionProblem class implemented in notebooks/03_sz_problems/3.2e_sz_problem.ipynb, deriving a SteadySubductionProblem class that implements the temperature equations for a steady-state, subduction zone. We will couple these equations to the isoviscous and dislocation creep rheology flow solutions in the next two notebooks.

5. Equations#

The function SubductionProblem.stokes_forms implements the Stokes equations for all our problems. Here we need to override the equivalent (and currently unimplemented) function for temperature (temperature_forms) given the steady-state temperature advection-diffusion equation

(128)#\[\begin{equation} \rho c_p \vec{v} \cdot \nabla T = \nabla \cdot \left( k \nabla T \right) + H \end{equation}\]

which we wish to convert into bilinear, \(S_T = S_T(T_t, T_a)\), and linear, \(f_T = f_T(T_t)\), forms, such that:

(129)#\[\begin{equation} S_T T = f_T \end{equation}\]

Due to the variation of the material parameters and the velocity functions across the domain, \(S_T\) is compiled from several integrals of different subregions of the domain

(130)#\[\begin{align} S_T =& \int_{\text{wedge}} \left[ T_t \rho_m \vec{v}_w\cdot\nabla T_a + \nabla T_t \cdot k_m\nabla T_a \right] dx \\ & + \int_{\text{slab}} \left[ T_t \rho_m \vec{v}_s\cdot\nabla T_a + \nabla T_t \cdot k_m\nabla T_a \right] dx \\ & + \int_{\text{crust}} \left[ \nabla T_t \cdot k_c\nabla T_a \right] dx \end{align}\]

Meanwhile \(f_T\) depends on whether the case has over-riding continental

(131)#\[\begin{equation} f_T = \int_{\text{upper crust}} T_t H_1 dx + \int_{\text{lower crust}} T_t H_2 dx \end{equation}\]

or oceanic

(132)#\[\begin{equation} f_T = 0 \end{equation}\]

crust.

We add the function temperature_forms to the SteadyIsoSubductionProblem class (newly derived from the base SubductionProblem class).

class SteadySubductionProblem(SubductionProblem):
    def temperature_forms(self):
        """
        Return the forms ST and fT for the matrix problem ST*T = fT for the steady-state 
        temperature advection-diffusion problem.

        Returns:
          * ST - lhs bilinear form for the temperature problem
          * fT - rhs linear form for the temperature problem
          * rT - residual linear form for the temperature problem
        """
        with df.common.Timer("Forms Temperature"):
            # set the crustal conductivity
            kc   = self.kc
            if self.sztype=='oceanic':
                # if we are oceanic then we use the mantle value
                kc   = self.km
            
            # advection diffusion in the slab
            STs = (self.T_t*self.rhom*self.cp*ufl.inner(self.vs_i, ufl.grad(self.T_a)) + \
                ufl.inner(ufl.grad(self.T_a), self.km*ufl.grad(self.T_t)))*self.dx(self.slab_rids)
            # advection diffusion in the wedge
            STw = (self.T_t*self.rhom*self.cp*ufl.inner(self.vw_i, ufl.grad(self.T_a)) + \
                ufl.inner(ufl.grad(self.T_a), self.km*ufl.grad(self.T_t)))*self.dx(self.wedge_rids)
            # just diffusion in the crust
            STc = ufl.inner(ufl.grad(self.T_a), kc*ufl.grad(self.T_t))*self.dx(self.crust_rids)
            # the complete bilinear form
            ST  = STs + STw + STc
            if self.sztype=='continental':
                # if the sztype is 'continental' then put radiogenic heating in the rhs form
                lc_rids = tuple([self.geom.crustal_layers['Crust']['rid']])
                uc_rids = tuple([self.geom.crustal_layers['UpperCrust']['rid']])
                fT  = self.T_t*self.H1*self.dx(uc_rids) + self.T_t*self.H2*self.dx(lc_rids)
            else:
                # if the sztype is 'oceanic' then create a zero rhs form
                zero_c = df.fem.Constant(self.mesh, df.default_scalar_type(0.0))
                fT = self.T_t*zero_c*self.dx
            # residual form
            # (created as a list of forms so we can assemble into a nest vector)
            rT = df.fem.form([ufl.action(ST, self.T_i) - fT])
        # return the forms
        return df.fem.form(ST), df.fem.form(fT), df.fem.form(rT)

6. Matrix-Vector System#

We will use the TemperatureSolver class, implemented in notebooks/3.2e_sz_problem.ipynb, that wraps a PETSc KSP linear solver to handle the assembly and solution of the temperature system in the next two notebooks.

Finish up#

Convert this notebook to a python module (saving first and ignoring markdown cells and those tagged as “main” or “ipy”).

from ipylab import JupyterFrontEnd
app = JupyterFrontEnd()
app.commands.execute('docmanager:save')
!jupyter nbconvert --TagRemovePreprocessor.enabled=True --TagRemovePreprocessor.remove_cell_tags="['main', 'ipy']" --TemplateExporter.exclude_markdown=True --TemplateExporter.exclude_input_prompt=True --TemplateExporter.exclude_output_prompt=True --NbConvertApp.export_format=script --ClearOutputPreprocessor.enabled=True --FilesWriter.build_directory=../../python/sz_problems --NbConvertApp.output_base=sz_steady_problem 3.3a_sz_steady_problem.ipynb

[NbConvertApp] Converting notebook 3.3a_sz_steady_problem.ipynb to script

[NbConvertApp] Writing 3103 bytes to ../../python/sz_problems/sz_steady_problem.py
Contents