Class FAMesh¶
- class flow.flowplex.FAMesh(*args, **kwargs)[source]¶
This class calculates drainage area in an implicit, iterative manner using PETSc solvers. It accounts for multiple flow direction paths (SFD to MFD) based on user input declaration.
Note
The class follows the parallel approach described in Richardson et al., 2014 where the iterative nature of the computational algorithms used to solve the linear system creates the possibility of accelerating the solution by providing an initial guess.
For drainage computation, the class uses a depression-less surface and computes river incision expressed using a stream power formulation function of river discharge and slope.
Initialise
__init__(*args, **kwargs)The initialisation of FAMesh class consists in the declaration of PETSc vectors and matrices.
Public Methods
This function is the main entry point for flow accumulation computation.
This function defines the flow direction matrices.
River incision is based on a standard form of the stream power law assuming detachment-limited behaviour.
Private Methods
_buildFlowDirection(h[, down])This function builds from neighbouring slopes the flow directions.
_distributeDownstream(pitVol, FA, hl, step)In cases where rivers flow in depressions, they might fill the sink completely and overspill or remain within the depression, forming a lake.
This function computes erosion rates in metres per year.
_matrix_build([nnz])Creates a sparse PETSc matrix.
_matrix_build_diag(V[, nnz])Builds a PETSc diagonal matrix based on a given array V
_solve_KSP(guess, matrix, vector1, vector2)PETSc scalable linear equations solvers (KSP) component provides Krylov subspace iterative method and a preconditioner.
Public functions¶
- FAMesh.flowAccumulation()[source]¶
This function is the main entry point for flow accumulation computation.
Note
Flow accumulation (FA) calculations are a core component of landscape evolution models as they are often used as proxy to estimate flow discharge, sediment load, river width, bedrock erosion, and sediment deposition. Until recently, conventional FA algorithms were serial and limited to small spatial problems.
gospl model computes the flow discharge from FA and the net precipitation rate using a parallel implicit drainage area (IDA) method proposed by Richardson et al., 2014 but adapted to unstructured grids.
It calls the following private functions:
_buildFlowDirection
_solve_KSP
_distributeDownstream
- FAMesh.matrixFlow()[source]¶
This function defines the flow direction matrices.
Note
The matrix is built incrementally looping through the number of flow direction paths defined by the user. It proceeds by assembling a local Compressed Sparse Row (CSR) matrix to a global PETSc matrix.
When setting up the flow matrix in PETSc, we preallocate the non-zero entries of the matrix before starting filling in the values. Using PETSc sparse matrix storage scheme has the advantage that matrix-vector multiplication is extremely fast.
The matrix coefficients consist of weights (comprised between 0 and 1) calculated based on the number of downslope neighbours and proportional to the slope.
- Parameters
flow – boolean to compute matrix for either downstream water or sediment transport
- FAMesh.riverIncision()[source]¶
River incision is based on a standard form of the stream power law assuming detachment-limited behaviour.
It calls the private function _getErosionRate described above. Once erosion rates have been calculated, the function computes local eroded thicknesses for the considered time step and update local elevation and cumulative erosion, deposition values.
If multiple lithologies are considered, the stratigraphic record is updated based on eroded thicknesses.
Important
The approach assumes that the volume of rock eroded using the stream power law accounts for both the solid and void phase.
Private functions¶
- FAMesh._buildFlowDirection(h, down=True)[source]¶
This function builds from neighbouring slopes the flow directions. It calls a fortran subroutine that locally computes for each vertice:
the indices of receivers (downstream) nodes depending on the desired number of flow directions (SFD to MFD).
the distances to the receivers based on mesh resolution.
the associated weights calculated based on the number of receivers and proportional to the slope.
- Parameters
h – elevation numpy array
- FAMesh._distributeDownstream(pitVol, FA, hl, step)[source]¶
In cases where rivers flow in depressions, they might fill the sink completely and overspill or remain within the depression, forming a lake. This function computes the excess of water (if any) able to flow dowstream.
Important
The excess water is then added to the downstream flow accumulation (FA) and used to estimate rivers’ erosion.
- Parameters
pitVol – volume of depressions
FA – excess flow accumulation array
hl – current elevation array
step – downstream distribution step
- Returns
pitVol, excess (updated volume in each depression and boolean set to True is excess flow remains to be distributed)
- FAMesh._getErosionRate()[source]¶
This function computes erosion rates in metres per year. This is done on the filled elevation. We use the filled-limited elevation to ensure that erosion is not going to be underestimate by small depressions which are likely to be filled (either by sediments or water) during a single time step.
The simplest law to simulate fluvial incision is based on the detachment-limited stream power law, in which erosion rate depends on drainage area \(A\), net precipitation \(P\) and local slope \(S\) and takes the form:
\[E = − \kappa P^d (PA)^m S^n\]\(\kappa\) is a dimensional coefficient describing the erodibility of the channel bed as a function of rock strength, bed roughness and climate, \(d\), \(m\) and \(n\) are dimensionless positive constants.
Default formulation assumes \(d = 0\), \(m = 0.5\) and \(n = 1\). The precipitation exponent \(d\) allows for representation of climate-dependent chemical weathering of river bed across non-uniform rainfall.
Important
In gospl, the coefficients m and n are fixed and the only variables that the user can tune are the coefficient d and the erodibility \(\kappa\). E is in m/yr and therefore the erodibility dimension is \((m yr)^{−0.5}\).
The erosion rate is solved by an implicit time integration method, the matrix system is based on the receiver distributions and is assembled from local Compressed Sparse Row (CSR) matrices into a global PETSc matrix. The PETSc scalable linear equations solvers (KSP) is used with both an iterative method and a preconditioner and erosion rate solution is obtained using PETSc Richardson solver (richardson) with block Jacobian preconditioning (bjacobi).
- FAMesh._matrix_build(nnz=(1, 1))[source]¶
Creates a sparse PETSc matrix.
Note
To achieve good performance during matrix assembly, the function preallocates the matrix storage by setting the array nnz.
- Parameters
nnz – array containing the number of nonzeros in the various rows
- Returns
sparse PETSc matrix
- FAMesh._matrix_build_diag(V, nnz=(1, 1))[source]¶
Builds a PETSc diagonal matrix based on a given array V
- Parameters
V – diagonal data array
nnz – array containing the number of nonzero blocks
- Returns
sparse PETSc matrix
- FAMesh._solve_KSP(guess, matrix, vector1, vector2)[source]¶
PETSc scalable linear equations solvers (KSP) component provides Krylov subspace iterative method and a preconditioner. Here, flow accumulation solution is obtained using PETSc Richardson solver (richardson) with block Jacobian preconditioning (bjacobi).
Note
The solver choice was made based on the convergence results from Richardson et al. (2014) but can be changed if better solver and preconditioner combinations are found.
Using such iterative method allows for an initial guess to be provided. When this initial guess is close to the solution, the number of iterations required for convergence dramatically decreases. Here the flow discharge solution from previous time step can be passed as an initial guess to the solver as discharge often exhibits little change between successive time intervals.
- Parameters
guess – Boolean specifying if the iterative KSP solver initial guess is nonzero (when provided it corresponds to the previous flow discharge values)
matrix – PETSc sparse matrix used by the KSP solver composed of diagonal terms set to unity (identity matrix) and off-diagonal terms (weights between 0 and 1). The weights are calculated based on the number of downslope neighbours (based on the chosen number of flow direction directions) and are proportional to the slope.
vector1 – PETSc vector corresponding to the local volume of water available for runoff during a given time step (e.g. voronoi area times local precipitation rate)
vector2 – PETSc vector corresponding to the unknown flow discharge values
- Returns
vector2 PETSc vector of the new flow discharge values