## Convection / flow modelling (with underworld)


  - Intro to underworld - quick overview of the basics
  - Solvers: Global flow from local source of buoyancy or bc. (Why MG)
  - Analytic solutions in Underworld and how to use them to estimate / test 
  - Scaling & Rayleigh number
  - Rayleigh number and Nusselt number
  - Meshes / interpolation ... checkpointing
  - Tracers / strain markers, analysis
  - Thermal / chemical convection and more 


### The basics

Underworld is configured to solve a couple Stokes / advection-diffusion problem on a finite element mesh. 

The basic tools we need are:
    
   - The mesh itself
   - Some unknowns living on the mesh
   - Some boundary conditions
   - A driving term for the flow
   - Constitutive terms (viscosity, elasticity, plasticity)
   - A way to solve the equations
   - Some means of analysing the result
   - A means of timestepping
     
### Equations

The dimensionless (indicated by primes) Navier-Stokes equation can be written

$$
	\frac{1}{\textrm{Pr}} \frac{D \mathbf{v}'  }{Dt' }  =
	\acute{\nabla}^2  \mathbf{v}'  -  \acute{\nabla} p' + \textrm{ Ra} T' \hat{\mathbf{z}}
$$

We will work in the infinite Prandtl number limit in which intertial terms vanish, and we will also assume incompressibility

$$
\nabla^2  \mathbf{v}  -  {\nabla} p = -\textrm{ Ra} T \hat{\mathbf{z}}
$$
$$
\nabla \cdot \mathbf{v} = 0
$$

The buoyancy force in the Stokes equation includes (but is not limited to) thermal variations which are themselves coupled to the solution of the Stokes equation through the thermal advection term.

$$
			\frac{\partial T}{\partial t} + \mathbf{v}\cdot \nabla T -
					\kappa \nabla^2 T = \frac{H}{C_p}
$$

   
### Hands on  - quasi-static problems
   
The simplest problem we can set up is a boundary driven flow


  - [Lid Driven Flow Notebook](020-LidDrivenFlow.ipynb) 
     - Exercise: [Point Driving Velocity](021-Exercise-PointSource.ipynb) 


The next thing to try is a buoyancy driven flow 


  - [Thermally-driven flow](030-ThermalBuoyancy.ipynb)
  

Which works pretty much as expected - by adding an appropriate function to the buoyancy term in the Stokes equation template.


### Hands on - time-stepping

Let's look at how to make a time-dependent problem. We start by looking at 
some passive tracers swirling around in a lid-driven flow.

   - [Lid Driven Flow with Tracers](040-LidDrivenFlowPlusTracers.ipynb) 
      - Exercise (take away): Kellogg & Turcotte, 1990 
      - Exercise: Streak plots    

Nice, but now we should let the flow be buoyancy driven and allow the buoyancy to drive the flow.

   - [Rayleigh Bénard Convection](050-RayleighBenardConvection.ipynb)
   - [Rayleigh Bénard Convection - pt2 ](051-RayleighBenardConvection-pt2.ipynb)
       - Exercise (take away): Nu - Ra and Vrms - Ra plots (again)
   - [Rayleigh Bénard Convection - pt3](052-RayleighBenardConvection-pt3.ipynb)
       - Exercise (take away): Arrhenius rheology and non-linear viscosity


## What have we not covered ?

There are many aspects of the `underworld` code and of modelling in general that we have not dealt with here, or have only touched upon in a general way. Some of these issues apply to all the different numerical methods and tools you might come across.

   - **Where can I find out more ?**  _Python codes include 'docstrings' which can be quite helpful and you should always take a look. uw.help() grabs the docstrings for you and leaves out some of the developer-specific complexities._

   - **Treat the code like a black box, how do you know when to trust the output ?** _In `underworld`, we treat this question very seriously and we bundle a whole collection of analytic Stokes flow problems in with the code ( `uw.help(uw.fn.analytic)` )_
   
   - **What are the trade-offs between resolution, accuracy, choice of solvers, speed, parallel performance etc etc ?** _You probably want to find the answer(s) to that question out for yourself, I just want to make you aware that it is one of the key questions in modeling. Having analytic or approximate solutions to hand is very helpful in this case. Also seeing if the results **converge** to a specific answer as you put more resources into the solution._ 

   - **Particles and material history.** _One of the reasons `underworld` was written was to look at problems of complex fluids - ones where the material history is relevant to the constitutive behaviour (e.g. anisotropy, viscoelasticity, strain-softening ...). There are examples of these in the `underworld` user guide and tutorials._
   
   - **Parallel ?** _Most of the time, the python scripts will run fine in parallel, some of the more convoluted things we saw in the examples are the result of guiding you to parallel safe methods (e.g. the `uw.fn.view.maxmin` which works in parallel when numpy.max() does not). Also, the whole `gLucifer` visualisation framework is designed for use in parallel. Have a look at the examples, they all work in parallel._
   
   - **Saving data.** _You can save and read back all field and swarm variables and this all works in parallel. The file formats are hdf5 and these can be read later (in serial if you can't process your results easily in parallel). This is important as you may / will want to checkpoint and restart jobs that take a long time to run. **Warning:** you will have to make sure your scripts handle whatever branches arise from restarting as we cannot snapshot your scripts and save them._
   
   - **Workflow.** _There are many, many python packages available which will just work with underworld. For example, to allow you to set up complicated pattens of particles using projections, mapping data to material properties or comparing observations and model outputs. Many are already built-in to `scipy`. You might end up with some issues in parallel if you are not careful (e.g. `scipy`'s fft routines), but perhaps you can create a workflow in which this processing is handled in serial, saves hdf5 information, and is then used to initialise a parallel script._




```python
import underworld as uw
from underworld import function as fn
import glucifer
import numpy as np
```


```python
print help(uw)
```

    Help on package underworld:
    
    NAME
        underworld
    
    FILE
        /Users/rfarrington/Underworld/Underworld2-development/underworld/__init__.py
    
    DESCRIPTION
        Underworld2 is a python-friendly version of the Underworld geodynamics
        code which provides a programmable and flexible front end to all the
        functionality of the code running in a parallel HPC environment. This
        gives signficant advantages to the user, with access to the power of
        python libraries for setup of complex problems, analysis at runtime,
        problem steering, and coupling of multiple problems.
        
        Underworld2 is integrated with the literate programming environment of
        the jupyter notebook system for tutorials and as a teaching tool for
        solid Earth geoscience.
        
        Underworld is an open-source, particle-in-cell finite element code
        tuned for large-scale geodynamics simulations. The numerical algorithms
        allow the tracking of history information through the high-strain
        deformation associated with fluid flow (for example, transport of the
        stress tensor in a viscoelastic, convecting medium, or the advection of
        fine-scale damage parameters by the large-scale flow). The finite
        element mesh can be static or dynamic, but it is not contrained to move
        in lock-step with the evolving geometry of the fluid. This hybrid approach
        is very well suited to complex fluids which is how the solid Earth behaves
        on a geological timescale.
    
    PACKAGE CONTENTS
        _net (package)
        _stgermain
        _uwid
        conditions (package)
        container (package)
        function (package)
        mesh (package)
        swarm (package)
        systems (package)
        utils (package)
    
    FUNCTIONS
        barrier()
            Creates an MPI barrier. All processes wait here for others to catch up.
        
        matplotlib_inline()
            This function simply enables Jupyter Notebook inlined matplotlib results.
            This function should be called at the start of your notebooks as a 
            replacement for the Jupyter Notebook *%matplotlib inline* magic. It provides
            the same functionality, however it allows notebooks to be converted to
            python without having to explicitly remove these calls.
        
        nProcs()
            Returns the number of processes being utilised by the simulation.
            
            Returns
            -------
            unsigned
                Number of processors.
        
        rank()
            Returns the rank of the current process.
            
            Returns
            -------
            unsigned
                Rank of current process.
    
    DATA
        __version__ = '2.4.0-dev'
    
    VERSION
        2.4.0-dev
    
    
    None



```python
## Tab completion:  fn.
## gives you a helpful list with all the available functions
## but it leaves out the _hidden_ ones which the following does not

print dir(fn.math)
```

    ['_Function', '__builtins__', '__doc__', '__file__', '__name__', '__package__', '_cfn', '_uw', 'abs', 'acos', 'acosh', 'asin', 'asinh', 'atan', 'atanh', 'cos', 'cosh', 'dot', 'erf', 'erfc', 'exp', 'log', 'log10', 'log2', 'pow', 'sin', 'sinh', 'sqrt', 'tan', 'tanh']



```python
print help(fn.analytic.SolCx)
```

    Help on class SolCx in module underworld.function.analytic:
    
    class SolCx(_SolBase)
     |  The boundary conditions are free-slip everywhere on a unit domain.
     |  There is a viscosity jump in the x direction at :math:`x=xc`.
     |  The flow is driven by the following density perturbation:
     |  
     |  .. math::
     |      \rho = \sin (nz \pi z) \cos (nx \pi x)
     |  
     |  Parameters
     |  ----------
     |  viscosityA : float
     |      Viscosty of region A.
     |  viscosityB : float
     |      Viscosty of region B.
     |  xc : float
     |      Location for viscosity jump.
     |  nx : unsigned
     |      Wavenumber in x axis.
     |  
     |  Method resolution order:
     |      SolCx
     |      _SolBase
     |      __builtin__.object
     |  
     |  Methods defined here:
     |  
     |  __init__(self, viscosityA=1.0, viscosityB=2.0, xc=0.25, nx=1, *args, **kwargs)
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors inherited from _SolBase:
     |  
     |  __dict__
     |      dictionary for instance variables (if defined)
     |  
     |  __weakref__
     |      list of weak references to the object (if defined)
     |  
     |  fn_bodyforce
     |  
     |  fn_pressure
     |  
     |  fn_strainRate
     |  
     |  fn_stress
     |  
     |  fn_velocity
     |  
     |  fn_viscosity
    
    None



```python
print help(uw.systems.Stokes)
```

    Help on class Stokes in module underworld.systems._stokes:
    
    class Stokes(underworld._stgermain.StgCompoundComponent)
     |  This class provides functionality for a discrete representation
     |  of the Stokes flow equations.
     |  
     |  Specifically, the class uses a mixed finite element method to
     |  construct a system of linear equations which may then be solved
     |  using an object of the underworld.system.Solver class.
     |  
     |  The underlying element types are determined by the supporting
     |  mesh used for the 'velocityField' and 'pressureField' parameters.
     |  
     |  The strong form of the given boundary value problem, for :math:`f`,
     |  :math:`g` and :math:`h` given, is
     |  
     |  .. math::
     |      \begin{align}
     |      \sigma_{ij,j} + f_i =& \: 0  & \text{ in }  \Omega \\
     |      u_{k,k} + \frac{p}{\lambda} =& \: H  & \text{ in }  \Omega \\
     |      u_i =& \: g_i & \text{ on }  \Gamma_{g_i} \\
     |      \sigma_{ij}n_j =& \: h_i & \text{ on }  \Gamma_{h_i} \\
     |      \end{align}
     |  
     |  where,
     |  
     |  * :math:`\sigma_{i,j}` is the stress tensor
     |  * :math:`u_i` is the velocity,
     |  * :math:`p`   is the pressure,
     |  * :math:`f_i` is a body force,
     |  * :math:`\lambda` is a bulk viscosity,
     |  * :math:`H` is the compressible equation source term,
     |  * :math:`g_i` are the velocity boundary conditions (DirichletCondition)
     |  * :math:`h_i` are the traction boundary conditions (NeumannCondition).
     |  
     |  The problem boundary, :math:`\Gamma`,
     |  admits the decompositions :math:`\Gamma=\Gamma_{g_i}\cup\Gamma_{h_i}` where
     |  :math:`\emptyset=\Gamma_{g_i}\cap\Gamma_{h_i}`. The equivalent weak form is:
     |  
     |  .. math::
     |      \int_{\Omega} w_{(i,j)} \sigma_{ij} \, d \Omega = \int_{\Omega} w_i \, f_i \, d\Omega + \sum_{j=1}^{n_{sd}} \int_{\Gamma_{h_j}} w_i \, h_i \,  d \Gamma
     |  
     |  where we must find :math:`u` which satisfies the above for all :math:`w`
     |  in some variational space.
     |  
     |  Parameters
     |  ----------
     |  velocityField : underworld.mesh.MeshVariable
     |      Variable used to record system velocity.
     |  pressureField : underworld.mesh.MeshVariable
     |      Variable used to record system pressure.
     |  fn_viscosity : underworld.function.Function
     |      Function which reports a viscosity value.
     |      Function must return scalar float values.
     |  fn_bodyforce : underworld.function.Function, Default = None
     |      Function which reports a body force for the system.
     |      Function must return float values of identical dimensionality
     |      to the provided velocity variable.
     |  fn_lambda : Removed use, fn_one_on_lambda instead
     |  fn_minus_one_on_lambda: underworld.function.Function, Default = None
     |      Function which defines a non solenoidal velocity field via the relationship
     |      div(velocityField) = -fn_minus_one_on_lambda * pressurefield + fn_source
     |      When this is left as None a incompressible formulation of the stokes equation is formed, ie, div(velocityField) = 0.
     |      fn_minus_one_on_lambda is incompatible with the 'penalty' stokes solver, ensure a
     |      'penalty' equal to 0 is used when fn_minus_one_on_lambda is used. By default this is the case.
     |  fn_source : underworld.function.Function, Default = None
     |      Function which defines a non solenoidal velocity field via the relationship
     |      div(velocityField) = -fn_minus_one_on_lambda * pressurefield + fn_source.
     |      fn_minus_one_on_lambda is incompatible with the 'penalty' stokes solver, ensure
     |      the 'penalty' of 0, is used when fn_lambda is used. By default this is the case.
     |  voronoi_swarm : underworld.swarm.Swarm
     |      If a voronoi_swarm is provided, voronoi type numerical integration is
     |      utilised. The provided swarm is used as the basis for the voronoi
     |      integration. If no voronoi_swarm is provided, Gauss integration
     |      is used.
     |  conditions : underworld.conditions.SystemCondition
     |      Numerical conditions to impose on the system. This should be supplied as
     |      the condition itself, or a list object containing the conditions.
     |  
     |  Notes
     |  -----
     |  Constructor must be called by collectively all processes.
     |  
     |  Method resolution order:
     |      Stokes
     |      underworld._stgermain.StgCompoundComponent
     |      underworld._stgermain.StgClass
     |      underworld._stgermain.LeftOverParamsChecker
     |      __builtin__.object
     |  
     |  Methods defined here:
     |  
     |  __init__(self, velocityField, pressureField, fn_viscosity, fn_bodyforce=None, fn_one_on_lambda=None, fn_lambda=None, fn_source=None, voronoi_swarm=None, conditions=[], _removeBCs=True, _fn_viscosity2=None, _fn_director=None, _fn_stresshistory=None, **kwargs)
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors defined here:
     |  
     |  eqResiduals
     |      Returns the stokes flow equations' residuals from the latest solve. Residual calculations
     |      use the matrices and vectors of the discretised problem.
     |      The residuals correspond to the momentum equation and the continuity equation.
     |      
     |      Return
     |      ------
     |      (r1, r2) - 2 tuple of doubles
     |          r1 is the momentum equation residual
     |          r2 is the continuity equation residual
     |      
     |      Notes
     |      -----
     |      This method must be called collectively by all processes.
     |  
     |  fn_bodyforce
     |      The body force function. You may change this function directly via this
     |      property.
     |  
     |  fn_one_on_lambda
     |      A bulk viscosity parameter
     |  
     |  fn_source
     |      The volumetric source term function. You may change this function directly via this
     |      property.
     |  
     |  fn_viscosity
     |      The viscosity function. You may change this function directly via this
     |      property.
     |  
     |  ----------------------------------------------------------------------
     |  Data and other attributes defined here:
     |  
     |  __abstractmethods__ = frozenset([])
     |  
     |  ----------------------------------------------------------------------
     |  Static methods inherited from underworld._stgermain.StgCompoundComponent:
     |  
     |  __new__(cls, *args, **kwargs)
     |      Creates stgermain instances of underlying objects.
     |      
     |      The list of required instances is provided as a dictionary ('_objectsDict') on the class, 
     |      with the entry key being the object name, and entry value being the object type.
     |      
     |      A second class data member ('_selfObjectName') provides the object name 
     |      which should be considered the object 'self'.
     |      
     |      For example:
     |      class Drawing(_stgermain.StgCompoundComponent):
     |          _selfObjectName = "_dr"  # child should set this
     |          _objectsDict = { "_cm": "lucColourMap",
     |                           "_dr": None            }
     |                           
     |      Note that this example is for an abstract class ('Drawing'), so it does not define 
     |      a type for '_dr', but instead defers to child classes to define the type of '_dr'.
     |      
     |      Args:
     |          objectsDictOverrule (dict): If provided, this will overwrite any objects 
     |                                      in the class level _objectDicts. Useful for
     |                                      programatically modifying behaviour.
     |          
     |      
     |      Returns:
     |          New created instance of child class.
     |  
     |  ----------------------------------------------------------------------
     |  Data and other attributes inherited from underworld._stgermain.StgCompoundComponent:
     |  
     |  __metaclass__ = <class 'underworld._stgermain._SetupClass'>
     |      This metaclass allows us to invoke a _setup method after the the __init__ method.
     |      Borrowed from:
     |      http://stackoverflow.com/questions/22261763/in-python-is-it-possible-to-write-a-method-that-will-be-automatically-called-aft
     |  
     |  ----------------------------------------------------------------------
     |  Methods inherited from underworld._stgermain.StgClass:
     |  
     |  __del__(self)
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors inherited from underworld._stgermain.LeftOverParamsChecker:
     |  
     |  __dict__
     |      dictionary for instance variables (if defined)
     |  
     |  __weakref__
     |      list of weak references to the object (if defined)
    
    None



```python

```

6ebd2a21d78c480c847184636be5e9180ce7acff
ipynb
Jupyter Notebook
notebooks/tutorials/ConvectionTutorial/Notebooks/010-Introduction.ipynb
lmoresi/uw-docker-landing-page
8c26e0cf05cbb3fbda311fc0bda1356f23293fe8
MIT
2019-06-17T08:30:23.000Z
2019-06-17T08:30:23.000Z
notebooks/tutorials/ConvectionTutorial/Notebooks/010-Introduction.ipynb
lmoresi/uw-docker-landing-page
8c26e0cf05cbb3fbda311fc0bda1356f23293fe8
MIT
notebooks/tutorials/ConvectionTutorial/Notebooks/010-Introduction.ipynb
lmoresi/uw-docker-landing-page
8c26e0cf05cbb3fbda311fc0bda1356f23293fe8
MIT
2020-06-23T13:02:27.000Z
2020-06-23T13:02:27.000Z
Qwen/Qwen-72B
1. YES
2. YES