Running goSPL#
Once the input file is ready, a goSPL simulation can be launched either from the command line or from a Python script — both run the same forward model.
Command line#
Installing goSPL provides a gospl console command:
# serial
gospl -i input.yml
# parallel (recommended) — choose any number of processes
mpirun -np 8 gospl -i input.yml -v
Options:
-i/--input(required) — the YAML input file.-v/--verbose— print per-step model progress (phase-grouped). The first line reports the goSPL version.--log— write the PETSc solver log summary.--profile— record a per-phase wall-clock profile (profile.json).--version— print the goSPL version and exit.
--i / --v are accepted as aliases of --input / --verbose. The
equivalent python -m gospl -i input.yml -v is available without a
console-script install. Note that fast mode is not a flag — set
domain: fast: true in the input file.
Python script#
The command above is a thin wrapper around the Model
class, so a run can equally be driven from Python (handy for parameter sweeps,
notebooks or coupling):
from gospl.model import Model as sim
# Read the input file (filename, verbose, showlog)
model = sim(args.input, args.verbose, args.log)
# Run the forward model
model.runProcesses()
# Clean up the PETSc objects
model.destroy()
Launch it under MPI the same way as any mpi4py program:
mpirun -np 8 python run_model.py -i input.yml
The number of processes used to run the model is independent of the number used later for post-processing.
Post-processing & utility commands#
Installing goSPL also provides several command-line tools (each is equally
python -m gospl.<module> if you haven’t reinstalled to register the
console scripts). They read the model output in the h5 directory; the
output fields page describes what each variable means.
Stratigraphic volume for ParaView — gospl-strata-volume#
Stacks the recorded layers into a 3-D wedge (triangular-prism) volume and
writes an XDMF (per-partition HDF5 + one .xdmf) for ParaView:
gospl-strata-volume --h5dir myrun/h5 # --field basic (default)
gospl-strata-volume --h5dir myrun/h5 --field lithology # or provenance
mpirun -np 4 gospl-strata-volume --h5dir myrun/h5 --field provenance \
--tout 100000 # ParaView time in yr
Key options: --outdir/--out (default strata), --field
basic|lithology|provenance, --tout/--tstart (label time in years),
--steps, --first-layer. Runs serially or under mpirun (partitions
split across ranks, independent of the run’s processor count).
Per-cell fields attached to every wedge: thickness, elevation
(wedge mid-height) and layer (the recorded stratigraphic-layer index)
always. --field basic (the default) adds just the per-layer
porosity and needs nothing beyond a plain stratigraphy run (no dual
lithology, no provenance) — use it for a single-lithology run. --field
lithology adds coarse/fine thickness,
fine fraction and per-fraction porosity; --field provenance adds the
per-class volume fraction (src_classN), the dominant source class and
the per-layer porosity (plus phiFine for a dual-lithology run). An
eroded / pinched-out layer (zero thickness, hence no source) inherits the
dominant source and composition of the cell directly below it in the
column, so it never renders as a no-source cell.
When the run used the groundwater / duricrust feature, every
field mode also attaches a per-layer induration (0–1) cell field — the
archived diagenetic crust degree stratDuri — so a cross-section shows buried
and exhumed duricrusts directly (no extra flag needed; absent for non-groundwater
runs).
Publication sections, wells & Wheeler — gospl-section#
Vector (PDF/SVG) stratigraphic figures via matplotlib:
# vertical cross-section coloured by lithology, with surface + basement lines
gospl-section --h5dir myrun/h5 --mesh input/mesh.npz:v:c --kind cross \
--along x --color-by lithology --out section.pdf
gospl-section ... --kind slice --z -50 # map at elevation z
gospl-section ... --kind well --xy 4.2e5,3.1e5 # synthetic borehole
gospl-section ... --kind wheeler --color-by thickness \
--strat-dt 1e4 # chronostratigraphic chart
# depositional facies from the deposition water depth, exaggerated x10,
# custom figure size, an interface line every 2 layers
gospl-section ... --kind cross --color-by facies --vexag 10 \
--figsize 12,3 --layer-lines 2
--color-by ∈ {deposition, thickness, lithology, coarse, porosity, age,
provenance, facies}; --along x|y or --path x0,y0;x1,y1;.... The
cross-section y-axis shows true elevation; --vexag applies a real
vertical exaggeration (data aspect, labels stay true); --figsize W,H sets
the figure size; --layer-lines N overlays a thin interface line every N
layers. ``facies`` classifies each layer by the water depth at
deposition (sea_level − stratZ): fluvial/deltaic plain (subaerial),
shoreface (0–20 m), distal offshore (20–50 m), upper slope (50–75 m), lower
slope (>75 m) — all tunable via --facies-depths (bin edges) and
--facies-colors; --color-by facies works for both cross and
wheeler (the Wheeler shoreline overlay is the subaerial↔marine boundary,
consistent with the facies). --figsize applies to every --kind. The
time axis of the Wheeler / age colouring is the real simulation
time: the per-layer interval is derived from the step’s stratal display time
(the .xmf Time) so the surface layer maps to that time —
--strat-dt YR overrides it. The sea-level datum (section line / Wheeler
shoreline trajectory / facies depth reference) is read from the simulation
for the step by default (the step’s .xmf sea constant);
--sea-level FLOAT overrides it. --legend-loc positions the legend
(e.g. 'upper left'), --title-fontsize sets the title size, and
--xlim MIN,MAX / --ylim MIN,MAX clip the cross-section / Wheeler to a
distance / elevation (or time) window (pinning a limit leaves the aspect auto,
so --vexag then has no effect). --figsize is honoured in the saved file
by default (pass
--tight to crop to content instead); for the cross-section, --vexag
exaggerates via the data aspect without overriding --figsize. Horizontal
distance is labelled in km and time in ky (the data stay in m / yr). The
synthetic well draws its colour bar horizontally at the base (few ticks,
small font, so it stays readable for a tall narrow figsize like 1,8);
well_panel() draws several wells on one
figure with a shared colour scale and a single colour bar (titled well 1,
well 2, … by default, or pass labels).
The same functions are importable for
notebooks (with extra keyword args — figsize, layer_lines,
facies_depths, facies_colors, facies_labels, legend_loc,
title_fontsize, xlim, ylim, well cbar_orientation / vmin /
vmax / colorbar / ylim, panel labels / vmin / vmax /
ylim)
(cross_section(),
horizontal_slice(),
synthetic_well(),
wheeler(),
well_panel()).
Regular-grid NetCDF for PyGMT / ArcGIS — gospl-grid#
Rasterises the surface to a CF-NetCDF grid with drainage basins, chi and drainage area; auto lon/lat for global meshes:
gospl-grid --h5dir myrun/h5 --mesh input/mesh.npz:v:c --out surface.nc \
--spacing 1000 --mn 0.5
Base level for catchments/chi defaults to the run’s sea level (--base-level
overrides); --latlim crops the polar caps (default 89.9). The sea level
used is written into the NetCDF as the global attribute sea_level and a
scalar sea_level variable (and returned in the result dict as
base_level), so the grid is self-describing. The priority-flood-filled
elevation is also written (filled). Every NetCDF variable carries its
units and a long_name definition (CF-style), so the file is
self-describing. The companion notebook API extracts /
plots per-basin river longitudinal profiles
(basin_rivers(),
plot_long_profile(),
plot_basin_map()). plot_long_profile plots
the raw elevation by default (which='elev' — keeps real lakes /
depressions + interpolation roughness as small peaks); which='filled' plots
the hydrologically-conditioned elevation, which is strictly monotonic
upstream. plot_basin_map overlays the sea-level coastline (the
elev == sea_level contour; defaults to the run’s sea level) and takes a
figsize.
Per-basin outflow fluxes — gospl-catchment#
For every drainage basin of a gridded surface, extract the cell of maximum
water discharge and the cell of maximum sediment load — i.e. each basin’s
river-mouth (outflow) point and its flux. It reads the gospl-grid
NetCDF directly (it uses the same variable names — FA, sedLoad,
basin, lon/lat — falling back to the legacy flowDischarge /
sedimentLoad / basinID names, so old files still work). The water flux
defaults to FA; pass --flow-var fillFA to use the depression-filled
accumulation instead (routes the trunk discharge through lakes):
gospl-catchment -i index.csv -o flowsed
The index CSV has two columns time,netcdf (one row per step):
time,netcdf
1,results/surface1.nc
5,results/surface5.nc
10,results/surface10.nc
For each step it writes flowsed/flow{time}.csv and flowsed/sed{time}.csv
(columns basin,lon,lat,val; val in m³/yr). Basins with <= --min-cells
cells (default 10) are skipped. The per-basin maximum is a single vectorised
lexsort (grouped arg-max), so a global 0.1° grid is processed in ~1 s per
step serially — the old MPI fan-out is no longer needed. From a notebook:
catchment_flux() (batch → CSVs) and
basin_outflow() (one file → two DataFrames).
When the run used the Level-B geochemistry, the gridded surface also carries
the dissolved-solute flux (riverSolute — the routed river dissolved load — and
its per-species fields riverSolute_<name>). gospl-catchment then also
writes flowsed/solute{time}.csv — for each basin, the solute outlet (cell
of maximum total solute flux) with columns basin,lon,lat,val plus one column
per species (e.g. carbonate, silica), each the species’ flux at that
outlet (they sum to val). The total field defaults to riverSolute, falling
back to the raw seepage export soluteflux (--solute-var overrides). From a
notebook, basin_solute_flux() returns that
DataFrame directly, and basin_outflow() adds a
"solute" entry when the field is present.
Sediment provenance — gospl-provenance#
Source-to-sink attribution (per-basin / per-pixel source mix, transport
distance, optional Cu-fertility layer) → HDF5 + XDMF + CSV. See
provenance for the inputs and --cu-weights.
ELA maps from paleo-temperature — gospl-ela#
A pre-processing helper that derives the glacier-geometry maps
(hela/hice) for the ice block from a per-vertex temperature map by
inverting the lapse rate:
gospl-ela --temperature t2m.npz --t-key t2m --lapse 0.0065 --t-ela -5 \
--elevation mesh.npz --z-key z --band 800 --out glaciers_0Ma
Then list the output in the ice.glaciers time series (hela: ['glaciers_0Ma',
'hela'] etc.; see Surface processes parameters). --reference surface|sealevel selects
whether the temperature map is at the surface or reduced to sea level.
Note
The console commands (gospl, gospl-strata-volume, gospl-section,
gospl-grid, gospl-catchment, gospl-provenance, gospl-ela)
appear after installing
goSPL. Until then, use the equivalent python -m gospl.<module> form. The
analysis tools need the optional extras: pip install gospl[analysis]
(numba, geopandas, netCDF4, matplotlib).
From a Jupyter notebook (import the API)#
Every tool is also importable, so it can be driven and the figures shown inline in a notebook (no CLI needed):
# --- stratigraphic sections / wells / Wheeler (matplotlib, inline) ---
from gospl.analyse.stratasection import (
load_strata, cross_section, horizontal_slice, synthetic_well, wheeler)
data = load_strata("myrun/h5", "input/mesh.npz", step=20)
cross_section(data, kind="x", color_by="lithology", vexag=20, sea_level=0)
synthetic_well(data, 4.2e5, 3.1e5, color_by="porosity")
wheeler(data, kind="x", color_by="thickness", sea_level=0, dt=1e4)
horizontal_slice(data, z=-50.0, color_by="age")
# --- gridded NetCDF + per-basin river profiles (PyGMT/ArcGIS + matplotlib) ---
from gospl.analyse.gridexport import (
grid_export, to_netcdf, basin_rivers, plot_long_profile, plot_basin_map)
g = grid_export("myrun/h5", "input/mesh.npz", step=20, spacing=1000.)
to_netcdf(g, "surface.nc") # for PyGMT / ArcGIS
riv = basin_rivers(g, basin_id=g["main_basin"], area_threshold=5e6)
plot_long_profile(riv); plot_basin_map(g, riv)
# --- ELA maps from temperature (pre-processing) ---
import numpy as np
from gospl.tools.ela_from_temperature import derive_ela
hela, hice = derive_ela(np.load("t2m.npz")["t2m"], 0.0065, -5.0, 800,
reference="surface", elevation=np.load("mesh.npz")["z"])
Each plotting call returns the matplotlib Axes/Figure (or PyGMT-ready
NetCDF) so you can tune it for a paper. The wedge volume (gospl-strata-volume)
is written for ParaView rather than inline rendering.