Contents

Landscape evolution

Three input files are provided in the recipe folder.

  • inputStrat.yml: records the stratigraphy but does not consider multiple lithologies

  • inputSedLay.yml: records the stratigraphy for multiple lithologies and uses the first sedimentary layers file defined in the Inputs definition notebook (sedlay25km.npz)

  • inputStrat.yml: records the stratigraphy for multiple lithologies and uses the second sedimentary layers file defined in the Inputs definition notebook (surflay25km.npz)

Run gospl

Running gospl is done by calling the runModel.py script with the name of the input file (below inputSedLay.yml) 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 inputSedLay.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.

Caution

This function will take approximately 5 minutes to complete with one CPU…

# On a single processor uncomment the following line...
%run script/runModel.py -i inputSedLay.yml

# In parallel...
#!mpirun -np 5 python3 script/runModel.py -i input.yml

--- Initialisation Phase (3.16 seconds)

+++ Output Simulation Time: 0.00 years

--- Computational Step                       (28.00 seconds)

+++ Output Simulation Time: 5000.00 years

--- Computational Step                       (23.37 seconds)

+++ Output Simulation Time: 10000.00 years

--- Computational Step                       (22.11 seconds)

+++ Output Simulation Time: 15000.00 years

--- Computational Step                       (18.77 seconds)

+++ Output Simulation Time: 20000.00 years

--- Computational Step                       (17.84 seconds)

+++ Output Simulation Time: 25000.00 years

--- Computational Step                       (16.71 seconds)

+++
+++ Total run time (130.00 seconds)
+++

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 outGoM).

Amongst the temporal variables outputed by gospl for this specific simulation, you will find:

  • surface elevation elev.

  • cumulative erosion & deposition values erodep.

  • flow accumulation flowAcc.

  • river sediment load sedLoad.

  • precipitation maps based on forcing conditions rain.

Loading libraries

gospl outputs are hdf5 files which are produced for each time steps and each partition (when run in parallel).

Using the readOutput.py script, we will first read these files and build for each recorded time the corresponding vtk file.

The function requires several arguments:

  • path: the path to the input file

  • filename: the name of the input file

  • step: the step you wish to output (here set to 5 corresponding to the last output based on the input parameters: start time 0 year, end time 50 thousand years with an output every 10 thousand years)

  • nbstep: the number of time steps to plot (useful if one want to output a netdcf file containing all time steps (done in the following section).

  • uplift_forcing: set to False as we are not considering any tectonic forcing

Here we will build the last time step (i.e. step set to 5):

import os
from script import readOutput as rout
from script import stratal as strat

step = 5
uplift_forcing = False

# Reading the final output generated by gospl
vtkgrid = rout.readOutput(path='./', filename='inputSedLay.yml', 
                          step=step, nbstep=None, 
                          uplift=uplift_forcing)

Post-processing

Now that the outputs from our simulation have been loaded on disk we will export them as a vtk grid using the exportVTK function.

The function computes the level of filling requires to have a depression-less DEM (it can be seen as different lakes levels). It also labels each catchment basin.

The created vtk file contains the following information:

  • elev: elevation,

  • erodep: cumulative erosion/deposition,

  • rain: rainfall,

  • FA: log of the flow accumulation,

  • SL: sediment load,

  • fill: filling for depression-less DEM, and

  • basin: each drainage basin label.

# Exporting the final output as a VTK mesh
out_path = 'export'
if not os.path.exists(out_path):
    os.makedirs(out_path)

    
vtkfile = os.path.join(out_path, 'GoMstrat'+str(step)+'.vtk')
vtkgrid.exportVTK(vtkfile, globNgh=True)

Writing VTK file export/GoMstrat5.vtk

Interactive plotting

We can now visualise the VTK output in the notebook directly using the pyvista library.

Using the top left widget you will be able to chose amongst the different outputed variables from gospl, define a colorscale and manually set the colorbar range for each of these variables.

import matplotlib
import numpy as np
import pyvista as pv

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'

mesh = pv.read(vtkfile)
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()
Inputs definition Extract stratigraphy