Class GridProcess

class tools.addprocess.GridProcess

When running goSPL on a 2D grid (i.e. not a global simulation), this class defines two processes solved directly on the unstructured mesh:

1. Flexural isostasy: it allows to compute isostatic deflections of Earth’s lithosphere with uniform or non-uniform flexural rigidity. Evolving surface loads are defined from erosion/deposition values associated to modelled surface processes.

2. Orographic rain: it accounts for change in rainfall patterns associated to change in topography. The orographic precipitation function is based on the Smith & Barstad (2004) linear model, solved directly on the unstructured mesh in parallel as two steady advection–relaxation equations (cloud water → hydrometeors → precipitation). The stratified airflow (mountain-wave) term of the original spectral solution is dropped, so the windward/lee rain-shadow is retained while no regular grid or FFT is required.

For global simulation, the library pyshtools provides a framework to estimate global-scale flexural isostasy.

Methods

applyFlexure()	This function computes the flexural isostasy equilibrium based on topographic change.
cptOrography()

Initialise

__init__()

Initialisation of the gridProcess class.

Public Methods

applyFlexure()	This function computes the flexural isostasy equilibrium based on topographic change.
cptOrography()

Private Methods

_buildOroMat(diag, offcoeff[, dirichlet])	Assemble a PETSc operator on the DMPlex from a per-row diagonal and the upwind off-diagonal coefficients returned by advecupwind (one column per neighbour, FVmesh_ngbID giving the neighbour indices).
_oroSolve(matrix, rhs, sol)	Solve one orographic advection-relaxation system with a dedicated (cached) GMRES + block-Jacobi KSP.
_updateTe()	Finds current elastic thickness map for the considered time interval.
_buildDHGrid()	Build the Driscoll-Healy (DH2, shape N x 2N) target grid for global flexure and precompute the interpolation weights linking the unstructured spherical mesh (self.mCoords) and the DH grid.
_unstr2dh(field)	Inverse-distance interpolation of a 1-D unstructured field defined at self.mCoords onto the DH2 grid (shape (dh_N, 2*dh_N)).
_dh2unstr(field_dh)	Bilinear interpolation of a DH2-grid array back to the unstructured mesh.
_cmptFlexGlobal(erodep_dh, te_dh[, ...])	Spectral thin-elastic-shell flexure solve on the precomputed Driscoll-Healy (DH2) grid.
_cmptFlexFEM(dzLocal)	Parallel finite-volume biharmonic flexure for FLAT (2D) models (flexure: method: 'fem' — the flat-model flexure solver).
_buildFlexFEM(Te)	Assemble (and cache) the flat-model FV biharmonic operator A = Lm·diag(D)·Lm + Δρg·I with its boundary conditions, and set up the reusable KSP.

Public functions

GridProcess.applyFlexure()

This function computes the flexural isostasy equilibrium based on topographic change. It is a simple routine that accounts for flexural isostatic rebound associated with erosional loading/unloading.

The function takes an initial (at time t) and final topography (at time t+Dt) (i.e. before and after erosion/deposition) and returns a corrected final topography that includes the effect of erosional/depositional unloading/loading. It uses a spectral method to solve the bi-harmonic equation governing the bending/flexure of a thin elastic plate floating on an inviscid fluid (the asthenosphere).

\[D (d^4 w / d^4 x ) + \Delta \rho g w = q\]

where \(D\) is the flexural rigidity, \(w\) is vertical deflection of the plate, \(q\) is the applied surface load, and \(\Delta \rho = \rho_m − \rho_f\) is the density of the mantle minus the density of the infilling material.

GridProcess.cptOrography()

Linear Theory of Orographic Precipitation following Smith & Barstad (2004), solved directly on the unstructured mesh in parallel.

The model is cast as two steady, vertically-integrated advection-relaxation equations for the cloud-water density \(q_c\) and the hydrometeor density \(q_s\), advected by a uniform wind \(\mathbf{v}\):

\[(\mathbf{v}\cdot\nabla + 1/\tau_c)\,q_c = C_w\,\mathbf{v}\cdot\nabla h, \qquad (\mathbf{v}\cdot\nabla + 1/\tau_f)\,q_s = q_c/\tau_c,\]

with the precipitation rate \(P = q_s/\tau_f\). The source is the terrain-forced uplift \(C_w\,\mathbf{v}\cdot\nabla h\): positive (condensation) on windward slopes, negative (evaporation) on lee slopes, which together with the downwind advection of \(q_c, q_s\) produces the windward-wet / lee-dry rain shadow. The elevation in the source is clamped to sea level, so submarine bathymetry produces no orographic forcing (the airflow follows the flat sea surface, not the seafloor).

In Fourier space these equations recover the Smith-Barstad transfer function with the stratified mountain-wave term \((1-i\,h_w m)\) set to one. That airflow term is therefore dropped (and the parameters nm, hw and latitude are inert); everything else is identical. The advection operator is the first-order upwind finite-volume scheme on the Voronoi mesh, so the operators depend only on the (uniform, constant) wind and are assembled once and cached.

Note

Please refer to the original manuscript of Smith and Barstad (2004) to understand the model physics and limitations.

Common bounds:

• precip_base : 0-10 [mm/h]

• precip_min : 0.001 - 1 [mm/h]

• conv_time : 200-2000 [s]

• fall_time : 200-2000 [s]

• cw : 0.001-0.02 [kg/m^3]

• rainfall_frequency : 0.1 - 24 [number of storms of 1 hour duration per day]

Private functions

GridProcess._buildOroMat(diag, offcoeff, dirichlet=False)

Assemble a PETSc operator on the DMPlex from a per-row diagonal and the upwind off-diagonal coefficients returned by advecupwind (one column per neighbour, FVmesh_ngbID giving the neighbour indices). This is the same incremental local-CSR assembly used for the tectonic advection matrix.

Parameters:

• diag – diagonal entries (length lpoints)

• offcoeff – off-diagonal coefficients, shape (lpoints, maxnb)

• dirichlet – when True, the (non-cyclic) domain edge rows are replaced by an identity row so the precipitation tracers are pinned to zero there (clean inflow / zero-padding equivalent).

Returns:

assembled PETSc matrix

GridProcess._oroSolve(matrix, rhs, sol)

Solve one orographic advection-relaxation system with a dedicated (cached) GMRES + block-Jacobi KSP. The operators are non-symmetric but diagonally-dominant M-matrices, so this converges quickly; the solution is warm-started from the previous step.

GridProcess._updateTe()

Finds current elastic thickness map for the considered time interval.

GridProcess._buildDHGrid()

Build the Driscoll-Healy (DH2, shape N x 2N) target grid for global flexure and precompute the interpolation weights linking the unstructured spherical mesh (self.mCoords) and the DH grid.

Forward direction (mesh -> DH): 3-D inverse-distance weighting via scipy.spatial.cKDTree on the existing Cartesian coordinates, which avoids polar / dateline singularities.

Backward direction (DH -> mesh): bilinear interpolation with precomputed integer indices and fractional offsets. The DH grid is implicitly extended at call time with a south-pole row and a wrap-around longitude column so every mesh point sits inside a 2 x 2 stencil.

Because the unstructured mesh is fixed for a given simulation, weights are built once and reused for every flexure step.

GridProcess._unstr2dh(field)

Inverse-distance interpolation of a 1-D unstructured field defined at self.mCoords onto the DH2 grid (shape (dh_N, 2*dh_N)).

GridProcess._dh2unstr(field_dh)

Bilinear interpolation of a DH2-grid array back to the unstructured mesh. The grid is padded at call time with a south-pole row (mean of the southernmost DH row) and a wrap-around longitude column.

GridProcess._cmptFlexGlobal(erodep_dh, te_dh, rho_infill=0.0, max_iter=50, flex_tol=0.0005, relax=1.0, anderson_depth=5)

Spectral thin-elastic-shell flexure solve on the precomputed Driscoll-Healy (DH2) grid.

Parameters:

erodep_dhnp.ndarray, shape (dh_N, 2*dh_N)

Erosion(-) / deposition(+) thickness in metres on the DH2 grid. Produced by _unstr2dh() from the unstructured mesh field.

te_dhfloat or np.ndarray of shape (dh_N, 2*dh_N)

Elastic thickness in metres. Scalar -> constant-Te spectral solve. Array (same grid as erodep_dh) -> iterative varying-Te solve.

rho_infillfloat

Density (kg/m^3) of the material filling the flexural moat (0 = air, 1030 = sea water).

max_iter, flex_tol, relax, anderson_depthint, float, float, int

Picard / Anderson iteration controls for the varying-Te branch. relax < 1 damps the update (useful for strong Te contrasts).

Returns:

np.ndarray, shape (dh_N, 2*dh_N)

Flexural deflection in metres on the DH2 grid. Sign: negative = down (subsidence under deposition; rebound gives positive w under erosion).

GridProcess._cmptFlexFEM(dzLocal)

Parallel finite-volume biharmonic flexure for FLAT (2D) models (flexure: method: 'fem' — the flat-model flexure solver).

Solves the plate equation with a single PETSc solve directly on the DMPlex: fully parallel, no gather-to-root, no regular grid, and varying elastic thickness in one linear solve (no iteration over the rigidity contrast). It replaced the former serial 'FD' and 'FFT' solvers, which gathered the load to one rank, solved on a regular grid, and interpolated back.

Solves the thin-plate-on-elastic-foundation equation with a spatially varying rigidity \(D(x)=E\,T_e(x)^3/[12(1-\nu^2)]\)

\[\nabla^2\!\big(D\,\nabla^2 w\big) + \Delta\rho\,g\,w = q, \qquad q = \rho_s\,g\,\mathrm{(erodep)}\]

Discretised with the cached finite-volume negative-Laplacian Lm (\(= -\nabla^2\), the hillslope stencil) applied twice, this is the single-field SPD-structured system

\[\big[\,L_m\,\mathrm{diag}(D)\,L_m + \Delta\rho g\,I\,\big]\,w = q,\]

solved with GMRES + GAMG (Lm is row-area-normalised, hence slightly non-symmetric, so GMRES rather than CG). The elastic foundation \(\Delta\rho g\,I\) (a positive diagonal) makes the operator well-posed — it pins the absolute deflection, so there is no nullspace and no mean removal (unlike the closed-sphere spectral solver). The solver defaults to a cached DIRECT factorisation — serial PETSc LU, parallel MUMPS — reused every step (only the RHS changes); it is options-prefixed (flexfem_) so an iterative Krylov method can be requested for meshes too large to factorise.

When this solver replaced the former regular-grid solver it agreed with it on a flat mesh where the deflection decays inside the domain: correlation 0.998 (natural BC), 0.9996 (clamped). Where the flexural wavelength approaches the domain size the two boundary discretisations differ by ~10 %.

Note

Boundary conditions (per side, from flex_bcN/S/E/W):

• 0Slope0Shear and Mirror — the natural zero-flux FV boundary (w'=0, w'''=0; a thin plate’s Mirror is 0Slope0Shear), imposed for free by the FV Laplacian — no modification.

• 0Displacement0Slope (clamped) — pin w=0 (Dirichlet zeroRows on that side’s nodes); the natural w'=0 from the inner FV Laplacian completes the clamp.

Sides map geographically. 0Moment0Shear and Periodic are not implemented.

The operator A and its KSP are CACHED and reused across steps: the geometry (Lm) is fixed and, for constant or piecewise-constant Te, A only changes when the Te interval advances. Each step then re-uses the factorisation / preconditioner and only the RHS changes — so after the one-off setup a step costs a back-substitution (serial) or a few warm Krylov iterations (parallel). Serial runs use a direct LU (fast on small meshes); parallel runs use GMRES+GAMG. The KSP is options-prefixed (flexfem_) to override.

Parameters:

dzLocal – local (lpoints) elevation change = surface load thickness (m).

Returns:

local (lpoints) flexural deflection (m), same sign convention as _cmptFlexGlobal (negative = subsidence under deposition).

GridProcess._buildFlexFEM(Te)

Assemble (and cache) the flat-model FV biharmonic operator A = Lm·diag(D)·Lm + Δρg·I with its boundary conditions, and set up the reusable KSP. Called from _cmptFlexFEM() on the first flexure step and whenever the elastic-thickness interval advances; the result is reused for every step in between (only the RHS changes).