Class UnstMesh#

class mesher.unstructuredmesh.UnstMesh[source]#

This class defines the 2D or spherical mesh characteristics and builds a PETSc DMPlex that encapsulates this unstructured mesh, with interfaces for both topology and geometry. The PETSc DMPlex is used for parallel redistribution and for load balancing.

Note

goSPL is built around a Finite-Volume method (FVM) for representing and evaluating partial differential equations. It requires the definition of several mesh variables such as:

the number of neighbours surrounding every node,

the cell area defined using Voronoi area,

the length of the edge connecting every nodes, and

the length of the Voronoi faces shared by each node with his neighbours.

In addition to mesh defintions, the class declares several functions related to forcing conditions (e.g. paleo-precipitation maps, tectonic (vertical and horizontal) displacements, stratigraphic layers…). These functions are defined within the UnstMesh class as they rely heavily on the mesh structure.

Important

The grid (2D or spherical) requires locally-orthogonal Voronoi/Delaunay staggering, or an unstructured C-grid type numerical formulation as described in Engwirda 2017

Finally a function to clean all PETSc variables is defined and called at the end of a simulation.

Methods

`applyForces`()	Finds the different values for climatic, tectonic and sea-level forcing that will be applied at any given time interval during the simulation.
`applyTectonics`()	Finds the different values for tectonic forces that will be applied at any given time interval during the simulation.
`destroy_DMPlex`()	Destroys PETSc DMPlex objects and associated PETSc local/global Vectors and Matrices at the end of the simulation.

Initialise

__init__()

The initialisation of UnstMesh class calls the private function _buildMesh.

Public Methods

`applyForces`()	Finds the different values for climatic, tectonic and sea-level forcing that will be applied at any given time interval during the simulation.
`applyTectonics`()	Finds the different values for tectonic forces that will be applied at any given time interval during the simulation.
`destroy_DMPlex`()	Destroys PETSc DMPlex objects and associated PETSc local/global Vectors and Matrices at the end of the simulation.

Private Methods

`_meshfrom_cell_list`(dim, cells, coords)	Creates a DMPlex from a list of cells and coordinates.
`_meshStructure`()	Defines the mesh structure and the associated voronoi parameters used in the Finite Volume method.
`_generateVTKmesh`(points, cells)	A global VTK mesh is generated to compute the distance between mesh vertices and coastlines position.
`_buildMesh`()	This function is at the core of the UnstMesh class.
`_gatherGlobalOnRoot`(local_full)	Assemble a full-mesh (`mpoints`) array on rank 0 from the owned nodes of a local (`lpoints`) array.
`_set_DMPlex_boundary_points`(label)	In case of a 2D mesh, this function finds the points that join the edges that have been marked as "boundary" faces in the directed acyclic graph (DAG) then sets them as boundaries.
`_get_boundary`([label])	In case of a 2D mesh, this function finds the nodes on the boundary from the DM.
`_readErosionDeposition`()	Reads existing cumulative erosion depostion from a previous experiment as defined in the YAML input file following the nperodep key.
`_updateRain`()	Finds current rain values for the considered time interval and computes the volume of water available for runoff on each vertex.
`_updateIce`()	Refresh the glacier-geometry fields (terminus / ELA / ice-cap altitude) for the current time from the `_iceTimeSeries` built in the input parser — the ice analogue of `_updateRain`.
`_resolveIceField`(field)	Materialise one glacier-geometry field (a `(scalar, map_spec)` pair, exactly one non-None) to a full-mesh array.
`_updateEvap`()	Finds current evaporation values for the considered time interval and stores them as a per-node m/yr numpy array on self.evapVal.
`_updateEroFactor`()	Finds current erodibility factor values for the considered time interval.

Public functions#

UnstMesh.applyForces()[source]#: Finds the different values for climatic, tectonic and sea-level forcing that will be applied at any given time interval during the simulation.

UnstMesh.applyTectonics()[source]#: Finds the different values for tectonic forces that will be applied at any given time interval during the simulation.

UnstMesh.destroy_DMPlex()[source]#: Destroys PETSc DMPlex objects and associated PETSc local/global Vectors and Matrices at the end of the simulation.

Private functions#

UnstMesh._meshfrom_cell_list(dim, cells, coords)[source]#

Creates a DMPlex from a list of cells and coordinates.

Note

PETSc DMPlex requires to be initialised on one processor before load balancing.

Parameters:

dim – topological dimension of the mesh
cells – vertices of each cell
coords – coordinates of each vertex

UnstMesh._meshStructure()[source]#

Defines the mesh structure and the associated voronoi parameters used in the Finite Volume method.

Important

The mesh structure is built locally on a single partition of the global mesh.

Once the voronoi definitions have been obtained a call to the fortran subroutine definetin is performed to order each node and the dual mesh components, it records:

all cells surrounding a given vertice,
all edges connected to a given vertice,
the triangulation edge lengths,
the voronoi edge lengths.

UnstMesh._generateVTKmesh(points, cells)[source]#: A global VTK mesh is generated to compute the distance between mesh vertices and coastlines position.

Note

The distance to the coastline for every marine vertices is used to define a maximum shelf slope during deposition.

The coastline contours are efficiently obtained from VTK contouring function.

UnstMesh._buildMesh()[source]#

This function is at the core of the UnstMesh class. It encapsulates both mesh construction (triangulation and voronoi representation for the Finite Volume discretisation), PETSc DMPlex distribution and several PETSc vectors allocation.

The function relies on several private functions from the class:

_generateVTKmesh
_meshfrom_cell_list
_meshStructure
_readErosionDeposition
readStratLayers

Note

It is worth mentionning that partitioning and field distribution from global to local PETSc DMPlex takes a lot of time for large mesh.

UnstMesh._gatherGlobalOnRoot(local_full)[source]#

Assemble a full-mesh (mpoints) array on rank 0 from the owned nodes of a local (lpoints) array. Returns the global array on rank 0 and None on every other rank.

This replaces the np.zeros(self.mpoints) + sentinel; arr[self.locIDs] = local; Allreduce(MPI.MAX) idiom for fields that are only consumed on rank 0 (the serial flexure / orography solves). With it, non-root ranks never allocate an mpoints-sized array, and the data moved is a single Gatherv of the owned values (~``mpoints`` total) instead of an mpoints-sized Allreduce on every rank — both the per-rank memory and the communication cost of these phases stop growing with the global mesh size on every process.

The owned nodes (self.glIDs) — not self.locIDs — are gathered: each global vertex is owned by exactly one rank, so they tile 0..mpoints-1 with no overlap, and the result is identical to the MAX reduction because a ghost node always carries the same value as its owner.

Parameters:: local_full – an lpoints (local, with ghosts) array.
Returns:: the assembled mpoints array on rank 0, None elsewhere.

UnstMesh._set_DMPlex_boundary_points(label)[source]#: In case of a 2D mesh, this function finds the points that join the edges that have been marked as “boundary” faces in the directed acyclic graph (DAG) then sets them as boundaries.

UnstMesh._get_boundary(label='boundary')[source]#: In case of a 2D mesh, this function finds the nodes on the boundary from the DM.

UnstMesh._readErosionDeposition()[source]#

Reads existing cumulative erosion depostion from a previous experiment as defined in the YAML input file following the nperodep key.

This functionality can be used when restarting from a previous simulation in which the mesh has been modified either to account for horizontal advection or to refine/coarsen a specific region during a given time period.

UnstMesh._updateRain()[source]#: Finds current rain values for the considered time interval and computes the volume of water available for runoff on each vertex.

Note

It is worth noting that the precipitation maps are considered as runoff water. If one wants to account for evaporation and infiltration you will need to modify the precipitation maps accordingly as a pre-processing step.

UnstMesh._updateIce()[source]#

Refresh the glacier-geometry fields (terminus / ELA / ice-cap altitude) for the current time from the _iceTimeSeries built in the input parser — the ice analogue of _updateRain.

Each interval supplies, per field, either a uniform scalar or a per-vertex [file, key] map; both are materialised here to full-mesh arrays (elaMesh / iceMesh / termMesh) which iceAccumulation indexes by locIDs. Maps are (re)loaded only when the active interval changes (step changes, like the precipitation maps).

UnstMesh._resolveIceField(field)[source]#: Materialise one glacier-geometry field (a (scalar, map_spec) pair, exactly one non-None) to a full-mesh array.

UnstMesh._updateEvap()[source]#

Finds current evaporation values for the considered time interval and stores them as a per-node m/yr numpy array on self.evapVal.

Mirrors _updateRain but only supports two source types: a eUniform scalar (m/yr) or an eMap`+`eKey pair pointing to a per-node array in an .npz file. Elevation-banded evaporation is intentionally not supported in v1 (see DESIGN_EVAPORATION.md D4).

Note

Evaporation is treated as a within-step sink at two points in the flow solver: (i) subtracted from rainfall before the IDA solve (channel runoff), and (ii) subtracted from per-pit inflow inside _distributeDownstream before the spillover decision (lake-surface evap). Lakes that cannot sustain themselves against the local evap budget simply do not form. The same cell-area assumption is used on land and over lake surfaces; users who need a stronger lake-surface rate should encode the spatial pattern directly in eMap.

UnstMesh._updateEroFactor()[source]#: Finds current erodibility factor values for the considered time interval.

Note

It is worth noting that the erodibility factor is an indice representing different lithological classes (see Moosdorf et al., 2018).