Note
Dual-lithology example¶
This example will run in approximately 20 minutes.
To run the notebook several post and pre-processing librairies are required:
panel
pyvista
pyevtk
jigsaw
To install these dependencies read the documentation in the user guide.
In addition, to these libraries several Python functions have been defined to ease the creation of gospl
input files and the visualisation of the ouputs. These functions are located in the script
folder in the same directory as this notebook:
buildMesh.py
creates the initial mesh using thejigsaw
libraryrunModel.py
runsgospl
based on the provided input filereadOutput.py
builds aVTK
file from model outputsstratal.py
extract the recorded stratigraphy in a specific region
[1]:
import numpy as np
import pyvista as pv
from script import stratal as strat
from script import readOutput as rout
import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1 import make_axes_locatable
label_size = 7
matplotlib.rcParams['xtick.labelsize'] = label_size
matplotlib.rcParams['ytick.labelsize'] = label_size
matplotlib.rc('font', size=6)
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
1. Input files: topography and rainfall¶
First we have a data
folder containing relevant inputs for the simulation.
First, we will use a 0.1 degree resolution
netcdf
topography grid at present day (freely available paleoDEM files can be obtained from different sources, an example being the gPlates webportal), andSecond, we will use a rainfall map based on current day precipitation, the dataset is again a
netcdf
file at 1.0 degree.
NOTE The folder is organised with specific sub-directory that are used later on during the mesh creation. You will need to follow this structure if you want to use your own dataset. Specifically the following folder names will be search: paleomap and precipitation.
2. Initial unstructured mesh generation¶
We start by creating the unstructured mesh using jigsaw
library. It is ran from the terminal using buildMesh.py
script.
python3 script/buildMesh.py -t=0 -d=data -s=5000,1000,100
This function takes 3 arguments:
t
the time interval in Ma of the starting simulation time (here 0 Ma as we start at present day),d
the folder containing the input conditions. The code assumes that the paleo-elevations are located in the folder data under paleomap and arenetCDF
files of the formXXMa.nc
withXX
the specified time. It also assumes that the displacement maps are located in velocity and are off the formvelocity_ XXMa.xy
. Lastly for the paleo-precipitation, the assumed file is under precipitation and arenetCDF
files of the formXXMa.nc
as well.s
is the space conditions for thejigsaw
algorithm and consists of 3 values: the spacing in km for the mesh in the deep ocean (<=-1000 m), the spacing in km across shelf margin (>=-1000 m and < 0m) and the spacing in km in the continental domain. Here, the resolution will vary from 5000 km to 1000 km and 100 km respectively. This is set at pretty low resolution to ensure a quicker runtime. We often use the following values-s=100,30,15
in our run.
This function will take approximtely 5 minutes to complete…
[2]:
%run script/buildMesh.py -t=0 -d=data -s=5000,1000,100
Read scotese map (0.01 seconds)
Writing Jigsaw topo input .msh file (15.27 seconds)
Load topography grid (1.92 seconds)
Build space function (6.85 seconds)
Perform triangulation (22.71 seconds)
Get unstructured mesh (0.98 seconds)
Creation goSPL files took (48.55 seconds)
From the previous script, you will get a new folder created called input0
that contains the mesh information and the associated discretised elevation and rainfall values. Specifically you will have:
input0/mesh0.vtk
: the VTK file containing the elevation at specified resolution,input0/0Ma.npz
: a Numpy compressed file with the mesh coordinates, cells, each vertice neighbours indices and the nodes elevation,input0/rain0Ma.npz
: a Numpy compressed file with the precipitation values for each mesh vertice.
Before going further, you can check the mesh validity by loading the created VTK
file (input0/mesh0.vtk
).
[3]:
mesh = pv.read('input0/mesh0.vtk')
elev = mesh.get_array(name='value')
earthRadius = 6.371e6
scale = 20.
factor = 1.+ (elev/earthRadius)*scale
mesh.points[:, 0] *= factor
mesh.points[:, 1] *= factor
mesh.points[:, 2] *= factor
contour = mesh.contour([0])
plot = pv.PlotterITK()
plot.add_mesh(mesh, scalars="value")
plot.add_mesh(contour, color="black", opacity=1.)
plot.show()
3. Creating the initial stratigraphic mesh¶
In addition to the previous input files, we will specify an initial stratigraphic mesh. Here we impose a simple uniform 50 km thick layer for the entire globe composed of a 30% of fine and 70% of coarse sediment.
NOTE More complex stratigraphy can be defined with multiple layers varying spatially in thickness and composition.
In addition, to the layer thicknesses and compositions, gospl
requires the porosity of each sediment type present in a given layer to be specified.
Here we follow the approach presented in the technical guide of the documentation, and the porosity is considered to varies with depth z.
[4]:
# Loading gospl mesh
loadMesh = np.load('input0/0Ma.npz')
gCoords = loadMesh["v"]
gZ = loadMesh["z"]
# Define layers variables
H = np.zeros((len(gZ),2)) # thickness
Z = H.copy() # elevation
Fperc = H.copy() # fine fraction
Fphi = H.copy() # fine porosity
Sphi = H.copy() # coarse porosity
H[:,1] = 50.0e3 # 50 km thick
Fperc[:,1] = 0.3 # 30% of fines
Z[:,0] = gZ - 50.e3 # elevation at the layer centre (base layer is assumed to have a 0 thickness)
Z[:,1] = gZ - 25.e3 # elevation at the layer centre (top layer centre elevation is at 25 km depth)
phis = 0.49 # Coarse sediment surface porosity
phif = 0.63 # Fine sediment surface porosity
z0s = 3700.0 # e-folding fine sediment thickness for porosity reduction 3700 m
z0f = 1960.0 # e-folding fine sediment thickness for porosity reduction 1960 m
# Compute porosity based on above equation
Fphi1 = phif * np.exp(-25.e3/z0f)
Sphi1 = phis * np.exp(-25.e3/z0s)
Fphi[:,1] = Fphi1
Sphi[:,1] = Sphi1
# Save the stratigraphic mesh as a Numpy file...
np.savez_compressed('input0/sed0Ma', strataH=H, strataF=Fperc, strataZ=Z, phiF=Fphi, phiS=Sphi)
4. Run gospl¶
Running gospl
is done by calling the runModel.py
script with the name of the input file as argument.
The Python script takes the following arguments:
-i XXXX.yml
specifying the input file name (required)-l True/False
for outputing PETSc log during runtime (default is set to False)-v True/False
for verbosing option during runtime (default is set to False)
You can open the input.yml
file to look at the parameters that are setup for this model. A complete list of the gospl
input variables is available in the user guide documentation.
This function will take approximately 5 minutes to complete with one CPU…
[5]:
# On a single processor...
%run script/runModel.py -i input.yml
# In parallel...
#!mpirun -np 4 python3 script/runModel.py -i input.yml
--- Initialisation Phase (0.85 seconds)
+++ Output Simulation Time: 0.00 years
--- Computational Step (2.28 seconds)
--- Computational Step (1.90 seconds)
--- Computational Step (2.15 seconds)
--- Computational Step (2.18 seconds)
--- Computational Step (1.93 seconds)
--- Computational Step (1.96 seconds)
--- Computational Step (2.29 seconds)
--- Computational Step (2.30 seconds)
--- Computational Step (2.18 seconds)
--- Computational Step (2.34 seconds)
+++ Output Simulation Time: 100000.00 years
--- Computational Step (2.50 seconds)
--- Computational Step (2.15 seconds)
--- Computational Step (2.20 seconds)
--- Computational Step (2.24 seconds)
--- Computational Step (2.16 seconds)
--- Computational Step (2.22 seconds)
--- Computational Step (2.26 seconds)
--- Computational Step (2.29 seconds)
--- Computational Step (2.44 seconds)
--- Computational Step (2.32 seconds)
+++ Output Simulation Time: 200000.00 years
--- Computational Step (2.71 seconds)
--- Computational Step (2.30 seconds)
--- Computational Step (2.32 seconds)
--- Computational Step (2.32 seconds)
--- Computational Step (2.32 seconds)
--- Computational Step (2.23 seconds)
--- Computational Step (2.32 seconds)
--- Computational Step (2.42 seconds)
--- Computational Step (2.34 seconds)
--- Computational Step (2.28 seconds)
+++ Output Simulation Time: 300000.00 years
--- Computational Step (2.89 seconds)
--- Computational Step (2.43 seconds)
--- Computational Step (2.51 seconds)
--- Computational Step (2.40 seconds)
--- Computational Step (2.46 seconds)
--- Computational Step (2.33 seconds)
--- Computational Step (2.50 seconds)
--- Computational Step (2.38 seconds)
--- Computational Step (2.46 seconds)
--- Computational Step (2.37 seconds)
+++ Output Simulation Time: 400000.00 years
--- Computational Step (3.14 seconds)
--- Computational Step (2.34 seconds)
--- Computational Step (2.60 seconds)
--- Computational Step (2.37 seconds)
--- Computational Step (2.55 seconds)
--- Computational Step (2.39 seconds)
--- Computational Step (2.47 seconds)
--- Computational Step (2.49 seconds)
--- Computational Step (2.50 seconds)
--- Computational Step (2.46 seconds)
+++ Output Simulation Time: 500000.00 years
--- Computational Step (3.30 seconds)
--- Computational Step (2.54 seconds)
--- Computational Step (2.56 seconds)
--- Computational Step (2.54 seconds)
--- Computational Step (2.56 seconds)
--- Computational Step (2.71 seconds)
--- Computational Step (2.71 seconds)
--- Computational Step (2.78 seconds)
--- Computational Step (2.81 seconds)
--- Computational Step (2.82 seconds)
+++ Output Simulation Time: 600000.00 years
--- Computational Step (3.73 seconds)
--- Computational Step (2.85 seconds)
--- Computational Step (2.83 seconds)
--- Computational Step (2.76 seconds)
--- Computational Step (2.80 seconds)
--- Computational Step (2.83 seconds)
--- Computational Step (2.94 seconds)
--- Computational Step (2.80 seconds)
--- Computational Step (2.99 seconds)
--- Computational Step (2.58 seconds)
+++ Output Simulation Time: 700000.00 years
--- Computational Step (4.07 seconds)
--- Computational Step (2.93 seconds)
--- Computational Step (3.06 seconds)
--- Computational Step (2.69 seconds)
--- Computational Step (3.09 seconds)
--- Computational Step (2.70 seconds)
--- Computational Step (3.11 seconds)
--- Computational Step (2.96 seconds)
--- Computational Step (3.02 seconds)
--- Computational Step (2.89 seconds)
+++ Output Simulation Time: 800000.00 years
--- Computational Step (4.13 seconds)
--- Computational Step (2.91 seconds)
--- Computational Step (3.04 seconds)
--- Computational Step (2.72 seconds)
--- Computational Step (2.94 seconds)
--- Computational Step (2.63 seconds)
--- Computational Step (2.94 seconds)
--- Computational Step (2.64 seconds)
--- Computational Step (2.96 seconds)
--- Computational Step (2.65 seconds)
+++ Output Simulation Time: 900000.00 years
--- Computational Step (4.18 seconds)
--- Computational Step (2.86 seconds)
--- Computational Step (2.94 seconds)
--- Computational Step (2.89 seconds)
--- Computational Step (3.04 seconds)
--- Computational Step (2.88 seconds)
--- Computational Step (2.96 seconds)
--- Computational Step (2.90 seconds)
--- Computational Step (3.21 seconds)
--- Computational Step (2.93 seconds)
+++ Output Simulation Time: 1000000.00 years
--- Computational Step (4.53 seconds)
+++
+++ Total run time (270.56 seconds)
+++
5. Visualisation in a notebook environment¶
The preferred way for visualising the model output is via Paraview
by loading the time series file called gospl.xdmf
available in the output folder (here called dual-lithologies
).
Amongst the temporal variables outputed by gospl
you will find:
surface elevation elev.
cumulative erosion & deposition values erodep.
flow accumulation flowAcc before pit filling.
flow accumulation fillAcc for depressionless surface.
river sediment load sedLoad.
fine sediment load sedLoadf when dual lithologies are accounted for.
uplift subsidence values if vertical tectonic forcing is considered uplift.
horizontal displacement values when considered hdisp.
precipitation maps based on forcing conditions rain.
Several filters, rendering and calculation can be done with Paraview
but are beyond the scope of this example.
Here you will use the readOutput.py
functions available in the script
folder to visualise directly the model output in the notebook at the final time step.
The function requires several arguments:
path
: the path to the input filefilename
: the name of the input filestep
: the step you wish to output (here set to 10 corresponding to the last output based on the input parameters: start time 0 year, end time 1 million years with an output every 0.1 million years)back
: set toFalse
as the simulation is not a backward forward modeluplift
: set toFalse
as we are not considering any tectonic forcing
[6]:
# Reading the final output generated by gospl
output = rout.readOutput(path='./', filename='input.yml',step=10, back=False, uplift=False)
# Exporting the final output as a VTK mesh
output.exportVTK('step10.vtk')
Writing VTK file step10.vtk
We can now visualise the VTK
output in the notebook directly:
[7]:
mesh = pv.read('step10.vtk')
elev = mesh.get_array(name='elev')
earthRadius = 6.371e6
scale = 20.
factor = 1.+ (elev/earthRadius)*scale
mesh.points[:, 0] *= factor
mesh.points[:, 1] *= factor
mesh.points[:, 2] *= factor
contour = mesh.contour([0])
plot = pv.PlotterITK()
plot.add_mesh(mesh, scalars="elev")
plot.add_mesh(contour, color="black", opacity=1.)
plot.show()
6. Extract stratigraphy¶
Finally, we will look at the recorded stratigraphy. The stratigraphic layer are recorded in gospl
as HDF5
files stored in the output folder as stratal.XX.pX.h5
where XX
is the output step and X
the processor number.
The following information are stored:
elevation at time of deposition, considered to be to the current elevation for the top stratigraphic layer stratZ.
thickness of each stratigrapic layer stratH accounting for both erosion & deposition events.
proportion of fine sediment stratF contains in each stratigraphic layer.
porosity of coarse sediment phiS in each stratigraphic layer computed at center of each layer.
porosity of fine sediment phiF in each stratigraphic layer computed at center of each layer.
We will use the stratal.py
function to extract the information above. It requires the following arguments: 1. path
: the path to the input file 2. filename
: the name of the input file 3. layer
: the stratal file you wish to output
[8]:
# Load the function and specify the input file
strati = strat.stratal(path='./', filename='input.yml', layer=10)
# Read the stratigraphic dataset
strati.readStratalData()
# Interpolate the spherical dataset on a lon/lat regular grid
# by specifying the desired resolution and interpolation neighbours
strati.buildLonLatMesh(res=0.1, nghb=3)
Maximum layer number: 10
Current number of sedimentary: 52
Start building regular stratigraphic arrays
Percentage of arrays built : [####################] 100.0% DONE
We can visualise the create maps directly by doing…
[9]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
ax.imshow(np.flipud(strati.zi[-1,:,:]), extent=(-180, 180, -90, 90), vmin=-8000, vmax=8000, cmap=cm.bwr)
ax.set(xlabel='Longitude', ylabel='Latitude', yticks=np.arange(-90,120,30), xticks=np.arange(-180,180,30))
ax.minorticks_on()
We will now extract the stratigraphic layer for a specific region by using the writeMesh
function which takes the following argument:
vtkfile
the output name of theVTK
stratigraphic mesh to createlats
latitude of the lower left and upper right corner of the region (specified between -90 and 90 degree)lons
longitude of the lower left and upper right corner of the region (specified between -180 and 180 degree)sigma
the standard deviation for Gaussian kernel as defined in SciPygaussian_filter
function.
The function returns the domain length in meters along the X and Y borders.
Here we chose a region around the Ganges–Brahmaputra–Meghna delta.
This function will take approximately 2 minutes to complete.
[10]:
length = strati.writeMesh(vtkfile='strati2D',
lats=[9,24],
lons=[82,97],
sigma=2.)
The function will build a VTK
structured mesh containing the stratigraphic record for the region.
Here we will set the slice at the center of the domain…
We can viusualise the stratigraphic layers in the notebook:
[11]:
mesh = pv.read('strati2D.vts')
mesh.set_active_scalars('layID')
threshold = mesh.threshold([2,52])
surface = mesh.threshold([50,52])
# Position cross-section at the center of the region
slices = threshold.slice_orthogonal(x=length[0]/2, y=length[-1]/2, z=-10)
scale_factor = 5
slices[0].points[:, -1] *= scale_factor
slices[1].points[:, -1] *= scale_factor
contours0 = slices[0].contour(np.linspace(0, 51, 52))
contours1 = slices[1].contour(np.linspace(0, 51, 52))
surface.points[:, -1] *= scale_factor
p = pv.PlotterITK()
p.add_mesh(surface, scalars="dep elev", opacity=0.25)
p.add_mesh(slices[0], scalars="percfine")
p.add_mesh(slices[1], scalars="percfine")
p.add_mesh(contours0, color="black", opacity=1.)
p.add_mesh(contours1, color="black", opacity=1.)
p.show()
[ ]: