Class soilSPL#
- class eroder.soilSPL.soilSPL(*args, **kwargs)[source]#
The class computes river incision expressed using a stream power formulation function of river discharge and slope also accounting for soil production.
A non-linear diffusion of soil based on soil thickness is also implemented in this class.
If the user has turned-on the sedimentation capability, this class will solve implicitly the stream power formulation accounting for a sediment transport/deposition term (Yuan et al, 2019).
Methods
For river-transported sediments reaching the marine realm, this function computes the related marine deposition diffusion.
Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediments in modulating erosion and deposition rate, considering non-linear slope dependency and accounting for soil production.
updateSoilThickness([deposition])Updates soil thickness through time from the increment in
self.tmp.Initialise
__init__(*args, **kwargs)Initialisation of soilSPL class.
Public Methods
updateSoilThickness([deposition])Updates soil thickness through time from the increment in
self.tmp.Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediments in modulating erosion and deposition rate, considering non-linear slope dependency and accounting for soil production.
For river-transported sediments reaching the marine realm, this function computes the related marine deposition diffusion.
Private Methods
_form_residual_soil(snes, h, F)The nonlinear system (SNES) at each time step is solved iteratively by assessing the residual of the SPL equation accounting for erosion, deposition (transport-limited) and soil production.
_monitorsoil(snes, its, norm)Non-linear SPL with soil production solver convergence evaluation.
_build_soil_snes([primary])Construct (and return) a soil-SPL SNES solver and its residual vector.
_soilErodibility(PA, surfK)(Re)build the dt-scaled stream-power erodibility coefficients
Kbr(bedrock) andK_soil(soil layer) for the CURRENTself.dtandself.hOldArray._soilSNESsolve(x, guess)Solve the soil-SPL SNES into
xfrom initial guessguess: the primary solver, then the complementary quasi-Newton / ngmres fallback on failure (same primary→fallback net as the flow KSP)._adaptiveSubstepSoil(x, PA, surfK)Adaptive sub-stepping for the soil-SPL solve when the full-Δt SNES diverges.
Solves the non-linear stream power law for the transport limited and soil case.
This function computes erosion deposition rates in metres per year and associated soil evolution.
_evalFunctionSoil(ts, t, x, xdot, f)The non-linear system for soil diffusion is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (
rosw)._evalJacobianSoil(ts, t, x, xdot, a, A, B)The non-linear system for soil diffusion is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (
rosw)._evalSolutionSoil(t, x)Evaluate the initial solution of the SNES system.
Public functions#
- soilSPL.updateSoilThickness(deposition=True)[source]#
Updates soil thickness through time from the increment in
self.tmp.Lsoilis a subaerial regolith cover, so the increment is not retained at subaqueous nodes (marine or ponded lake) — see_subaqueousMask— and is frozen under ice (_iceFrozenMask). These gates run in every mode, so callers invoke this unconditionally.depositiondistinguishes the two kinds of increment:deposition=True (lake/pit, marine): in regolith mode the deposit is stratigraphy, not soil, so the increment is skipped (the gates still run — a cell newly ponded by this step’s deposition is re-zeroed consistently with lumped mode). In lumped mode it is added.
deposition=False (soil creep / hillslope transport): always added — creep moves the weathering regolith itself, in both modes.
- soilSPL.erodepSPLsoil()[source]#
Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediments in modulating erosion and deposition rate, considering non-linear slope dependency and accounting for soil production.
It calls the private function _getEroDepRateSoil described above. Once erosion/deposition rates have been calculated, the function computes local thicknesses and soil evolution for the considered time step and update local elevation and cumulative erosion, deposition values.
- soilSPL.diffuseSoil()[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 \times (1.0 - e^{-h_s/H_0} \nabla h \right)\]It calls the following private functions:
_evalFunctionSoil
_evalJacobianSoil
_evalSolutionSoil
Note
PETSc SNES and time stepping TS approaches are used to solve the non-linear equation above over the considered time step.
Private functions#
- soilSPL._form_residual_soil(snes, h, F)[source]#
The nonlinear system (SNES) at each time step is solved iteratively by assessing the residual of the SPL equation accounting for erosion, deposition (transport-limited) and soil production.
- soilSPL._monitorsoil(snes, its, norm)[source]#
Non-linear SPL with soil production solver convergence evaluation.
- soilSPL._build_soil_snes(primary=True)[source]#
Construct (and return) a soil-SPL SNES solver and its residual vector.
Two configurations are available:
primary (
primary=True) – a Nonlinear GMRES accelerator (ngmres) right-preconditioned by Nonlinear Richardson (nrichardson). The Richardson sweep is what actually applies the Krylov solve (cg) and the multigrid preconditioner (HYPRE BoomerAMG by default, orself.soil_pc); a barengmresignores the KSP/PC entirely, which is why the previous setup stalled toSNES_DIVERGED_MAX_ITon the stiff soil-production residual.fallback (
primary=False) – a limited-memory quasi-Newton solver (qn, L-BFGS) with a matrix-free critical-point line search. It builds an approximate Jacobian from secant updates (no analytic Jacobian required) and is markedly more robust for stiff problems; it runs only when the primary solver fails to converge.
Each SNES gets its own PETSc options prefix so per-solver options (e.g. the line-search type) do not leak into the model’s other SNES/KSP objects.
- Parameters:
primary – select the primary (True) or fallback (False) solver.
- Returns:
the configured
(SNES, residual Vec)pair.
- soilSPL._soilErodibility(PA, surfK)[source]#
(Re)build the dt-scaled stream-power erodibility coefficients
Kbr(bedrock) andK_soil(soil layer) for the CURRENTself.dtandself.hOldArray. The erosion limiter follows the current elevation (drop to the receiver), so this is recomputed per sub-step by the adaptive sub-stepper as well as once for the full step.
- soilSPL._soilSNESsolve(x, guess)[source]#
Solve the soil-SPL SNES into
xfrom initial guessguess: the primary solver, then the complementary quasi-Newton / ngmres fallback on failure (same primary→fallback net as the flow KSP). Returns the final converged reason (>= 0converged,< 0both failed). Prints only the informational “fallback recovered” line; the caller decides what to do on an outright failure (sub-step, then discard).
- soilSPL._adaptiveSubstepSoil(x, PA, surfK)[source]#
Adaptive sub-stepping for the soil-SPL solve when the full-Δt SNES diverges. The divergence is a stiffness problem: a de-armoured node that has just captured a large discharge demands an enormous single-step incision (
Kbr·SⁿwithKbr ∝ Δt·A^m), so the slope term overshoots and the residual grows. BecauseKbr ∝ Δt, splitting the step intoNincrements ofΔt/Ncuts the demanded per-solve incision byNuntil the residual is smooth enough to converge; summing the increments recovers the FULL step’s erosion (unlike simply discarding).Escalates
Nover_soilSubstepN(default 4, 8, 16); accepts the firstNwhose every sub-step converges. If none does, the step is discarded (revert to the prior elevation) as a last resort. On success the step-start bookkeeping (hOld/hOldArray/soilH/dt) is restored and one final residual evaluation rebuildsnsoilH/_soilDepoGrowthconsistently with the full step (used by the soil update). Returns the final reason (0sub-stepped OK,-1discarded).
- soilSPL._solveSoil()[source]#
Solves the non-linear stream power law for the transport limited and soil case. This calls the following private function:
_form_residual_soil
Note
PETSc SNES approach is used to solve the nonlinear equation without forming an analytic Jacobian. The primary solver is a Nonlinear GMRES accelerator (
ngmres) right-preconditioned by Nonlinear Richardson (nrichardson), whose Krylov solve uses a Preconditioned Conjugate Gradient (cg) method with a multi-grid preconditioner (HYPRE BoomerAMG by default). If the primary solver stalls on the stiff soil-production residual, a limited-memory quasi-Newton fallback (qn, L-BFGS) with a critical-point line search is used. The iteration budget, tolerances and preconditioner are configurable through the YAMLsoilblock (maxIter,rtol,atol,pcType).
- soilSPL._getEroDepRateSoil()[source]#
This function computes erosion deposition rates in metres per year and associated soil evolution. This is done on the filled elevation.
The approach is based on BasicHySa governing equations from Terrainbento (as described in Appendix B20 from Barnhart et al. (2019)).
Note
The approach uses a continuous layer of soil-alluvium, which influences both hillslope and river-induced erosion. It relies on the SPACE algorithm of Shobe et al. (2017).
- soilSPL._evalFunctionSoil(ts, t, x, xdot, f)[source]#
The non-linear system for soil diffusion is solved iteratively using PETSc time stepping and SNES solution and is based on Rosenbrock W-scheme (
rosw).Here again, we evaluate the residual function on a DMPlex for an implicit time-stepping method.