Time-Dependent Subduction Zone Setup

Time-Dependent Subduction Zone Setup#

Authors: Cameron Seebeck, Cian Wilson

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

  1. we will describe the subduction zone geometry and tesselate it intro 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 space 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

We have already seen and tested the implementation of all these steps for the steady state cases. In this section we consider steps 5-7 for the time-dependent versions. First, in this notebook we derive a new base time-dependent TDSubductionProblem class from the SubductionProblem class we implemented in notebooks/03_sz_problems/3.2e_sz_problem.ipynb.

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, StokesSolverNest, TemperatureSolver

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)

TDSubductionProblem class#

We build on the SubductionProblem class implemented in notebooks/03_sz_problems/3.2e_sz_problem.ipynb, deriving a TDSubductionProblem class that implements time-dependent versions of the temperature equations.

5. Equations#

The function SubductionProblem.stokes_forms implements the Stokes equations for all our problems, including the time-dependent cases. Here we need to override the equivalent function for temperature (temperature_forms) given the time-dependent temperature advection-diffusion equation

(133)#\[\begin{equation} \rho \left( \frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right) = \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:

(134)#\[\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

(135)#\[\begin{align} S_T =& \int_{\text{wedge}} \left[ T_t\rho_m \tilde{T}^{n+1} + \Delta t^n \theta T_t \tilde{\vec{v}}_w\cdot\nabla\tilde{T}^{n+1} + \Delta t^n \theta \nabla T_t \cdot k_m\nabla\tilde{T}^{n+1} \right] dx \nonumber \\ &~ + \int_{\text{slab}} \left[ T_t\rho_m \tilde{T}^{n+1} + \Delta t^n \theta T_t \tilde{\vec{v}}_s\cdot\nabla\tilde{T}^{n+1} + \Delta t^n \theta \nabla T_t \cdot k_m\nabla\tilde{T}^{n+1} \right] dx \nonumber \\ &~ + \int_{\text{crust}} \left[ T_t\rho_c\tilde{T}^{n+1} + \Delta t^n \theta \nabla T_t \cdot k_c\nabla\tilde{T}^{n+1}\right] dx \end{align}\]

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

(136)#\[\begin{align} f_{T_o} =& \int_{\text{wedge}} \left[ T_t\rho_m \tilde{T}^n - \Delta t^n \left(1-\theta\right) T_t \tilde{\vec{v}}_w\cdot\nabla\tilde{T}^n - \Delta t^n \left(1-\theta\right) \nabla T_t \cdot k_m\nabla\tilde{T}^n \right] dx \nonumber \\ &~ + \int_{\text{slab}} \left[ T_t\rho_m \tilde{T}^n - \Delta t^n \left(1-\theta\right) T_t \tilde{\vec{v}}_s\cdot\nabla\tilde{T}^n - \Delta t^n \left(1-\theta\right) \nabla T_t \cdot k_m\nabla\tilde{T}^n \right] dx \nonumber \\ &~ + \int_{\text{crust}} \left[ T_t\rho_c \tilde{T}^n - \Delta t^n \left(1-\theta\right) \nabla T_t \cdot k_c\nabla\tilde{T}^n\right] dx \\ f_T =&~f_{T_o} \end{align}\]

or continental

(137)#\[\begin{equation} f_T = f_{T_o} + \int_{\text{upper crust}} \Delta t^n T_t H_1 dx + \int_{\text{lower crust}} \Delta t^n T_t H_2 dx \end{equation}\]

crust.

The time-dependent equations add new parameters for the time-stepping

  • \(\theta \in [0,1]\) controlling the implicitness and accuracy of our solution - common choices are \(\theta=0\) (“explicit Euler”), \(\theta=1\) (“implicit Euler”), and \(\theta=0.5\) (“Crank-Nicolson”)

  • \(\Delta t^n = t^{n+1} - t^n\) controlling the time-step between the old, \(t^n\), and new, \(t^{n+1}\) times at which the old and new solutions, \(\tilde{T}^{n}\) and \(\tilde{T}^{n+1}\) respectively, are evaluated

For convenience we add these parameters to the members of our class by overriding the members function (this isn’t strictly required but is useful for bookkeeping).

class TDSubductionProblem(SubductionProblem):
    def members(self):
        super().members()

        # timestepping options
        self.theta = None
        self.dt    = None

We then override the function temperature_forms to return the time-dependent weak forms of the temperature equation as outlined above.

class TDSubductionProblem(TDSubductionProblem):
    def temperature_forms(self):
        """
        Return the forms ST and fT for the matrix problem ST*T = fT for the time dependent
        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"):
            # integration measures that know about the cell and facet tags

            # set the crustal conductivity and density
            kc   = self.kc
            rhoc = self.rhoc
            if self.sztype=='oceanic':
                # if we are oceanic then we use the mantle values
                kc   = self.km
                rhoc = self.rhom

            # advection diffusion in the slab
            STs = self.T_t*self.rhom*self.cp*(self.T_a+self.dt*self.theta*ufl.inner(self.vs_i, ufl.grad(self.T_a)))*self.dx(self.slab_rids) + \
                self.dt*self.km*self.theta*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_a))*self.dx(self.slab_rids)
            
            # advection diffusion in the wedge
            STw = self.T_t*self.rhom*self.cp*(self.T_a+self.dt*self.theta*ufl.inner(self.vw_i, ufl.grad(self.T_a)))*self.dx(self.wedge_rids) + \
                self.dt*self.km*self.theta*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_a))*self.dx(self.wedge_rids)
            
            # just diffusion in the crust
            STc = self.T_t*rhoc*self.cp*(self.T_a)*self.dx(self.crust_rids) + \
                self.dt*kc*self.theta*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_a))*self.dx(self.crust_rids)
            
            # the complete bilinear form
            ST  = STs + STw + STc
            
            fTs = self.T_t*self.rhom*self.cp*(self.T_n-self.dt*(1-self.theta)*ufl.inner(self.vs_i, ufl.grad(self.T_n)))*self.dx(self.slab_rids) - \
                    (1-self.theta)*self.dt*self.km*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_n))*self.dx(self.slab_rids)
            
            fTw = self.T_t*self.rhom*self.cp*(self.T_n-self.dt*(1-self.theta)*ufl.inner(self.vw_i, ufl.grad(self.T_n)))*self.dx(self.wedge_rids) - \
                    (1-self.theta)*self.dt*self.km*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_n))*self.dx(self.wedge_rids)
            
            fTc = self.T_t*rhoc*self.cp*self.T_n*self.dx(self.crust_rids) - \
                    (1-self.theta)*self.dt*kc*ufl.inner(ufl.grad(self.T_t), ufl.grad(self.T_n))*self.dx(self.crust_rids)

            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']])
                
                fTc += self.T_t*self.dt*self.H1*self.dx(uc_rids) + self.T_t*self.dt*self.H2*self.dx(lc_rids)
            
            fT = fTs + fTw + fTc


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

As before 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_tdep_problem 3.5a_sz_tdep_problem.ipynb

[NbConvertApp] Converting notebook 3.5a_sz_tdep_problem.ipynb to script

[NbConvertApp] Writing 4394 bytes to ../../python/sz_problems/sz_tdep_problem.py
Contents