Class hillSLP#

class sed.hillslope.hillSLP(*args, **kwargs)[source]#

This class encapsulates all the functions related to hillslope (soil creep) processes (both linear and non-linear cases are implemented).

Methods

getHillslope()

This function is the main entry point to compute linear and non-linear hillslope processes.

Initialise

__init__(*args, **kwargs)

The initialisation of hillSLP class consists in the declaration of several PETSc vectors.

Public Methods

getHillslope()

This function is the main entry point to compute linear and non-linear hillslope processes.

Private Methods

_hillSlope([smooth])

This function computes hillslope using a linear diffusion law commonly referred to as soil creep:

_buildDiffMat(Kd)

Assemble the implicit linear-diffusion operator \((I - L)\) for the finite-volume Laplacian, given per-node coefficients Kd (m², i.e. a diffusivity × time, or the mesh-size-based marine-smoothing coefficient N · cell_area).

_assembleDiffMat(coeffs)

Assemble a PETSc finite-volume stencil matrix from a per-node coefficient array coeffs of shape (lpoints, 1 + maxnb): column 0 is the diagonal, columns 1..maxnb are the off-diagonal entries for the maxnb neighbours in self.FVmesh_ngbID.

_assembleDiffMatCSR(coeffs)

Single-pass CSR assembly of the FV stencil matrix from coeffs (lpoints, 1+maxnb) (col 0 = diagonal, cols 1..maxnb = neighbour entries for self.FVmesh_ngbID).

_makeDiffusionKSP(prefix)

Create a cached KSP for a linear diffusion / smoothing solve: PETSc Richardson + block-Jacobi (matching flowplex._solve_KSP), nonzero initial guess, with the per-rank ILU pivots shifted so the block-Jacobi PCSetUp cannot fail on a degenerate / ocean-only partition.

_solveSmooth(rhs, sol)

Solve one implicit step of the cached marine-smoothing diffusion operator (see _hillSlope smooth=2).

_diff_nl_monitor(snes, its, norm)

Non-linear diffusion solver convergence evaluation.

_form_residual_nl_hillslope(snes, h, F)

The nonlinear system (SNES) at each time step is solved iteratively by assessing the residual of the hillslope equation.

_hillSlopeNL()

This function computes hillslope using a non-linear diffusion law:

_evalFunctionMardDiff(ts, t, x, xdot, f)

The non-linear system for freshly-deposited marine sediments diffusion at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (rosw).

_evalJacobianMardDiff(ts, t, x, xdot, a, A, B)

The non-linear system for freshly-deposited marine sediments diffusion at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (rosw).

_evalSolutionMardDiff(t, x)

Evaluate the initial solution of the SNES system.

_diffuseOcean(dh)

For river-transported sediments reaching the marine realm, this function computes the related marine deposition diffusion.

_diffuseImplicit(dh, mask, Cd_val[, label])

Solve the non-linear thickness diffusion equation

_diffuseImplicitPicard(dh, mask, Cd_val[, label])

Lagged-diffusivity (Picard) alternative to _diffuseImplicit for the non-linear thickness diffusion \(\partial_t h = \nabla\cdot(C_d(h) \nabla h)\) over one model step.

Public functions#

hillSLP.getHillslope()[source]#

This function is the main entry point to compute linear and non-linear hillslope processes.

Private functions#

hillSLP._hillSlope(smooth=0)[source]#

This function computes hillslope using a linear diffusion law commonly referred to as soil creep:

\[\frac{\partial z}{\partial t}= \kappa_{D} \nabla^2 z\]

in which \(\kappa_{D}\) is the diffusion coefficient and can be defined with different values for the marine and land environments (set with hillslopeKa and hillslopeKm in the YAML input file). It encapsulates, in a simple formulation, processes operating on superficial sedimentary layers. Main controls on variations of \(\kappa_{D}\) include substrate, lithology, soil depth, climate and biological activity.

Note

The hillslope processes in goSPL are considered to be happening at the same rate for coarse and fine sediment sizes.

Parameters:

smooth – smoothing mode. 0 (default) computes linear soil-creep hillslope and updates the elevation in place. 1 returns a one-off smoothed copy of the ice-surface proxy used to derive glacial flow directions. 2 returns a smoothed bathymetry used only to derive coherent marine flow directions in _matOcean — it never alters the elevation. The smooth=2 operator uses a mesh-size- based, timestep-independent diffusivity (see MARINE_SMOOTH_N_*) and is cached + reused across steps, rebuilt only when the coastline moves.

Note

This method uses scratch Vecs self.tmp / self.tmpL (and, for smooth=1, self.tmp1 as the input field). smooth=0 also uses self.hOld.

hillSLP._buildDiffMat(Kd)[source]#

Assemble the implicit linear-diffusion operator \((I - L)\) for the finite-volume Laplacian, given per-node coefficients Kd (m², i.e. a diffusivity × time, or the mesh-size-based marine-smoothing coefficient N · cell_area).

Shared by the linear hillslope (smooth=0), the ice-surface smoothing (smooth=1) and the cached marine-smoothing operator (smooth=2).

Parameters:

Kd – per-node diffusion coefficient array (lpoints, m²).

Returns:

the assembled PETSc diffusion matrix (caller owns it; the smooth=0/1 callers destroy it after the solve, the smooth=2 caller caches it).

hillSLP._assembleDiffMat(coeffs)[source]#

Assemble a PETSc finite-volume stencil matrix from a per-node coefficient array coeffs of shape (lpoints, 1 + maxnb): column 0 is the diagonal, columns 1..maxnb are the off-diagonal entries for the maxnb neighbours in self.FVmesh_ngbID. Used by _buildDiffMat (sethillslopecoeff output) and the lagged- diffusivity marine solver (jacobiancoeff output, scaled to I + dt·L).

Parameters:

coeffs – (lpoints, 1+maxnb) diagonal + neighbour coefficients.

Returns:

the assembled PETSc matrix (caller owns / destroys it).

hillSLP._assembleDiffMatCSR(coeffs)[source]#

Single-pass CSR assembly of the FV stencil matrix from coeffs (lpoints, 1+maxnb) (col 0 = diagonal, cols 1..maxnb = neighbour entries for self.FVmesh_ngbID). Equivalent to _assembleDiffMat (diag + per-neighbour axpy) but builds the matrix in ONE setValuesLocalCSR pass — ADD_VALUES sums shared columns exactly, mirroring seaplex._matOcean (#450). Much cheaper than the 12-axpy loop, used in the Picard hot loop (_diffuseImplicitPicard), which rebuilds the operator every iteration.

_buildDiffMat deliberately keeps the 12-axpy _assembleDiffMat so the bit-faithful cached smoother / linear-hillslope operators (#457 / #459) are not perturbed by a different floating-point summation order.

Parameters:

coeffs – (lpoints, 1+maxnb) diagonal + neighbour coefficients.

Returns:

the assembled PETSc matrix (caller owns / destroys it).

hillSLP._makeDiffusionKSP(prefix)[source]#

Create a cached KSP for a linear diffusion / smoothing solve: PETSc Richardson + block-Jacobi (matching flowplex._solve_KSP), nonzero initial guess, with the per-rank ILU pivots shifted so the block-Jacobi PCSetUp cannot fail on a degenerate / ocean-only partition. The options prefix scopes that shift (and any user override) to this solver only. Shared by the cached marine smoother (_solveSmooth) and the cached linear hillslope (_hillSlope smooth=0).

Parameters:

prefix – PETSc options prefix for this KSP (e.g. "marsmooth_").

Returns:

the configured PETSc KSP (caller caches and owns it).

hillSLP._solveSmooth(rhs, sol)[source]#

Solve one implicit step of the cached marine-smoothing diffusion operator (see _hillSlope smooth=2).

The operator (self._smoothMat) is rebuilt only when the coastline moves, so the block-Jacobi/ILU preconditioner is factorised once and reused across every step with an unchanged land/sea mask — that factorisation (not the 1-iteration Richardson solve) was the dominant per-step cost of the previous build-from-scratch approach.

Parameters:

  • rhs – PETSc global Vec — the field to smooth (read).

  • sol – PETSc global Vec — the smoothed result (written; scratch self.tmp at the call site).

hillSLP._diff_nl_monitor(snes, its, norm)[source]#

Non-linear diffusion solver convergence evaluation.

hillSLP._form_residual_nl_hillslope(snes, h, F)[source]#

The nonlinear system (SNES) at each time step is solved iteratively by assessing the residual of the hillslope equation.

hillSLP._hillSlopeNL()[source]#

This function computes hillslope using a non-linear diffusion law:

\[\frac{\partial h}{\partial t}= \nabla \cdot \left( C_d(h) \nabla h \right)\]

Two flavors of non-linear diffusion are possible based on user defined parameters:

  1. a non-critical hillslope model following the work of Wang et al. (2024).

  2. a non-linear depth-dependent creep law as described in Barnhart et al. (2019) (section 3.4.5).

Note

In this implementation of the SNES, we do not form the Jacobian and PETSc calculates it based on the residual function. Here, a Nonlinear Generalized Minimum Residual method is used ngmres, a Preconditioned Conjugate Gradient cg method is defined for the KSP and the preconditioner allows for multi-grid methods based on the HYPRE BoomerAMG approach.

hillSLP._evalFunctionMardDiff(ts, t, x, xdot, f)[source]#

The non-linear system for freshly-deposited marine sediments diffusion at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (rosw).

Here we evaluate the residual function on a DMPlex for an implicit time-stepping method.

hillSLP._evalJacobianMardDiff(ts, t, x, xdot, a, A, B)[source]#

The non-linear system for freshly-deposited marine sediments diffusion at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (rosw).

Here, we define the Jacobian matrix A and the preconditioner matrix B on a DMPlex.

hillSLP._evalSolutionMardDiff(t, x)[source]#

Evaluate the initial solution of the SNES system.

hillSLP._diffuseOcean(dh)[source]#

For river-transported sediments reaching the marine realm, this function computes the related marine deposition diffusion. It is based on a non-linear diffusion approach.

\[\frac{\partial h}{\partial t}= \nabla \cdot \left( C_d(h) \nabla h \right)\]

It calls the following private functions:

  • _evalFunctionMardDiff

  • _evalJacobianMardDiff

  • _evalSolutionMardDiff

Note

PETSc SNES and time stepping TS approaches are used to solve the non-linear equation above over the considered time step.

Parameters:

dh – Numpy Array of incoming marine depositional thicknesses

hillSLP._diffuseImplicit(dh, mask, Cd_val, label='diffusion')[source]#

Solve the non-linear thickness diffusion equation

\[\partial_t h = \nabla \cdot (C_d(h) \nabla h)\]

on the cells indicated by mask for one global timestep self.dt. Cells outside mask get a zero diffusion coefficient and therefore do not move. The cached PETSc TS (rosw) solver is reused across calls; self.Cd and self.minDiff_vec are swapped in-place so the same solver can be applied to marine sediments, lake sediments, or any other sub-domain.

Parameters:

  • dh – per-node initial deposit thickness (m)

  • mask – bool array (lpoints) or integer index array; cells where diffusion is active.

  • Cd_val – scalar non-linear diffusion coefficient (m^2/yr)

  • label – log label (e.g. “marine” or “lake”) for the timing print

Returns:

ndepo - smoothed deposit thickness (m), zeroed where negative

hillSLP._diffuseImplicitPicard(dh, mask, Cd_val, label='diffusion')[source]#

Lagged-diffusivity (Picard) alternative to _diffuseImplicit for the non-linear thickness diffusion \(\partial_t h = \nabla\cdot(C_d(h) \nabla h)\) over one model step.

Instead of the adaptive non-linear TS, this takes self.picardSub backward-Euler sub-steps of size dt/picardSub; within each sub-step it freezes the diffusivity \(C_d(h^k)\) and solves the resulting linear system \((I + \Delta t_{sub} L(C_d^k)) h^{k+1} = h^{start}\) with self.picardIts Picard updates. The operator L is built from jacobiancoeff with a zero derivative term (Kp=0), so it is the exact linear \(\nabla\cdot(C_d\nabla)\) operator consistent with the TS residual fctcoeff — including the “no flux across a zero- diffusivity face” marine-mask gating.

Each solve is linear (a cached richardson+bjacobi KSP, no SNES) and, because \(C_d\) is frozen, smooth — so it avoids the error-estimate rejections the adaptive TS hits at the \(C_d\) kink (dh<0.1). picardSub is the accuracy/speed knob. Same signature / return as _diffuseImplicit.

Parameters:

  • dh – per-node initial deposit thickness (m)

  • mask – cells where diffusion is active

  • Cd_val – scalar non-linear diffusion coefficient (m^2/yr)

  • label – log label

Returns:

ndepo - smoothed deposit thickness (m), zeroed where negative