Function veloren_world::sim::erosion::erode

source · 
fn erode(
    map_size_lg: MapSizeLg,
    h: &mut [f64],
    b: &mut [f64],
    wh: &mut [f64],
    max_uplift: f32,
    max_g: f32,
    kdsed: f64,
    _seed: &RandomField,
    rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
    uplift: impl Fn(usize) -> f32 + Sync,
    n_f: impl Fn(usize) -> f32 + Sync,
    m_f: impl Fn(usize) -> f32 + Sync,
    kf: impl Fn(usize) -> f64 + Sync,
    kd: impl Fn(usize) -> f64,
    g: impl Fn(usize) -> f32 + Sync,
    epsilon_0: impl Fn(usize) -> f32 + Sync,
    alpha: impl Fn(usize) -> f32 + Sync,
    is_ocean: impl Fn(usize) -> bool + Sync,
    height_scale: impl Fn(f32) -> f64 + Sync,
    k_da_scale: impl Fn(f64) -> f64,
    threadpool: &ThreadPool,
)
Expand description

Erode all chunks by amount.

Our equation is:

dh(p) / dt = uplift(p)−k * A(p)^m * slope(p)^n

where A(p) is the drainage area at p, m and n are constants (we choose m = 0.4 and n = 1), and k is a constant. We choose

k = 2.244 * uplift.max() / (desired_max_height)

since this tends to produce mountains of max height desired_max_height; and we set desired_max_height = 1.0 to reflect limitations of mountain scale.

This algorithm does this in four steps:

  1. Sort the nodes in h by height (so the lowest node by altitude is first in the list, and the highest node by altitude is last).

  2. Iterate through the list in reverse. For each node, we compute its drainage area as the sum of the drainage areas of its “children” nodes (i.e. the nodes with directed edges to this node). To do this efficiently, we start with the “leaves” (the highest nodes), which have no neighbors higher than them, hence no directed edges to them. We add their area to themselves, and then to all neighbors that they flow into (their “ancestors” in the flow graph); currently, this just means the node immediately downhill of this node. As we go lower, we know that all our “children” already had their areas computed, which means that we can repeat the process in order to derive all the final areas.

  3. Now, iterate through the list in order. Whether we used the filling method to compute a “filled” version of each depression, or used the lake connection algorithm described in [1], each node is guaranteed to have zero or one drainage edges out, representing the direction of water flow for that node. For nodes i with zero drainage edges out (boundary nodes and lake bottoms) we set the slope to 0 (so the change in altitude is uplift(i)).

    For nodes with at least one drainage edge out, we take advantage of the fact that we are computing new heights in order and rewrite our equation as (letting j = downhill[i], A[i] be the computed area of point i, p(i) be the x-y position of point i, flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1):

    h[i](t + dt) = hi + δt * (uplift[i] + flux(i) * h[j](t + δt)) / (1 + flux(i) * δt).

Since we compute heights in ascending order by height, and j is downhill from i, h[j] will always be the new h[j](t + δt), while h[i] will still not have been computed yet, so we only need to visit each node once.

Afterwards, we also apply a hillslope diffusion process using an ADI (alternating direction implicit) method:

https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90

We also borrow the implementation for sediment transport from

https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90

The approximate equation for soil production function (predicting the rate at which bedrock turns into soil, increasing the distance between the basement and altitude) is taken from equation (11) from [2]. This (among numerous other sources) also includes at least one prediction that hillslope diffusion should be nonlinear, which we sort of attempt to approximate.

[1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al.. Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion. Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175. ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩

[2] William E. Dietrich, Dino G. Bellugi, Leonard S. Sklar, Jonathan D. Stock Geomorphic Transport Laws for Predicting Landscape Form and Dynamics. Prediction in Geomorphology, Geophysical Monograph 135. Copyright 2003 by the American Geophysical Union 10.1029/135GM09