Solving PDEs in Geosciences Using Python

1
Solving PDEs in Geosciences Using Python

• Lutz Gross
• ESSCC _at_ The University of Queensland
• l.gross_at_uq.edu.au

2
Coworkers

• Ken Steube
• Artak Amirbekyan
• Peter Hornby (CSIRO EM)
• Hans Muhlhaus

3
Introduction to escript
4
Model of the Earths Mantel
5
Fault System Dynamics
6
Input Data
Topography bathymetry
geological map

• fault data base

7
Problem Characteristics

• Partial differential equation (PDE)
• time-depended
• non-linear
• complex
• coupled
• 3D (FEM or similar)
• Input parameters
• constants
• output from other PDE
• data bases, typically linked to geometry

8
escript Objectives

• Environment for PDE-based modeling
• Works for complex geometries
• Efficient for large scale 3D problems (?
  parallelized)
• Easy to use for non-software-engineers Supports
  code reuse
• access to analysis capabilities (e.g.
  visualization)

9
escript extension of Python

• Python programming language
• Interpreting
• Interactive
• Object-oriented
• Extendable
• see www.python.org

10
Approach Abstraction
11
escript Modules

• Data Handling escript
• PDEs escript.linearPDEs
• Geometries pycad
• Visualization pyvisi
• Models modelframe

12
Architecture
Interactive/script
GUI
Escript
modelframe
Model
escript
linearPDEs
pyvisi
pycad
finley
vtk
libescript.so
gmsh
13
escript

• Handels data on a Domain (FEM mesh)
• With spatial distribution
• stored on sample points
• tensors up to order 4
• data parallel OpenMP MPI mixed mode

14
Smoothness

• Smoothness defined by FunctionSpace attribute

15
Data handling

• Various internal data representation are used
• Constant
• Piecewise constant/Tagged
• Expanded
• Appropriate translation in expressions

16
escript (cont.)

• Transfer data to and from libraries
• e.g. coefficients and solution of PDEs
• via API (representation neutral)
• dump/load of data
• Archiving via NetCDF

17
escript.linearPDE
Interface to a general linear PDE for
u Implemented by finley solver
18
Finley Library

• C-library to solve general, linear PDE
• Finite Elements
• Unstructured, isoparametric meshes
• Contact Elements
• Implements escript.linearPDE
• Parallelized for OpenMP MPI mixed

19
Temperature Diffusion Model
20
Implementation

• def temperature(mydomain,Tinit,kappa,rhocp,Q)
• myPDELinearPDE(mydomain)
• mypde.setValue(AkappaKronecker(), \
  Drhocp/dt) TTinit
• t0.
• while tltt_end
• mypde.setValue(YQrhocp/dtT)
  TmyPDE.getSolution()
• tdt

21
Interested ?

• Documentation and other stuff
• http//shake200.uq.edu.au/twiki/bin/view/ESSCC/Esy
  sUser
• download Version 1
• supports OpenMP
• OSL v3 licence
• MPI version later in 2008

22
Application to Geosciences
23
escript Applications

• Convection of the Earths Mantel
• With complex rheology !
• Subduction zones
• Volcanic Lava Flow
• Reactive Transport
• Crustal Fault Dynamics
• Tsunami

24
Mantel Convection (Muhlhaus)
velocity vi, pressure p, temperature T
25
Convection Results
26
(No Transcript)
27
Level Set Method
Fgt0
Flt0
F0
28
Level Set Method (cont.)
Reinitialization (ever 1-10 time steps)
Linearization
29
Saddle Point Problem
30
Uzawa Iteration
Eliminate velocity
Apply preconditioner (Elman et al. 1996 for
constant ?)
31
Uzawa Iteration (cont.)
Implemented as
Needs PCG acceleration for variable ?
32
escript.LinearPDE class
2nd order system of linear PDE on a Domain O
33
escript implementation
from escript import from escript.linearPDEs
import LinearPDE def solve(dom,f,eta,v,p,stop) p
de1LinearPDE(dom) pde1.setValue(Aetaltsome
tensorgt,Yf) pde2LinearPDE(dom) pde2.setReducti
onOn() LBB condition ! pde2.setValue(D1./eta)
while not stop(p) pde1.setValue(X-etasymmetr
ic(grad(v)) \ -pkronecker(dom)) vpde1.sol
ve() pde2.setValue(Y-div(v)) ppde2.solve()
return v,p
34
Temperature Advection-Diffusion

• Explicit Taylor-Galerkin
• Watch courant condition
• Easy implementation in escript

35
Remarks

• Fails for vi,i?0
• Tends to be over-diffusive
• Requires fine meshes pushes dt
• solver on the C/C level ?

36
New escript PDE interface
FEM
Flux controlled transport solver in C
37
Algebraic up-winding
HOP
LOP (extrema diminishing)
38
Flux Control (Kuzmin et al. 2001)

• Criteria
• Mass conservation fij-fji and aijaij
• 0 aij 1
• if aij1 ? HOP
• if aij0 ? LOP
• Where local min/max of u could decrease/grow
  further gt aij0
• Leads to a weakly non-linear problem
• Use sparse matrix-vector product technology

39
Pure Advection Problem
40
Results (Ra105, A22)
41
Summary

• escript is an environment for PDE-based numerical
  modelling
• Geosciences in escript
• Stokes type problems
• Transport problems
• Level Set
• And still a lot of things to do