veloren_world/sim/
erosion.rs

1use super::{diffusion, downhill, uphill};
2use crate::{config::CONFIG, util::RandomField};
3use common::{
4    terrain::{
5        MapSizeLg, NEIGHBOR_DELTA, TerrainChunkSize, neighbors, uniform_idx_as_vec2,
6        vec2_as_uniform_idx,
7    },
8    vol::RectVolSize,
9};
10use tracing::{debug, error, info, warn};
11// use faster::*;
12use itertools::izip;
13use noise::NoiseFn;
14use num::{Float, Zero};
15use ordered_float::{FloatCore, NotNan};
16use rayon::prelude::*;
17use std::{
18    cmp::{Ordering, Reverse},
19    collections::BinaryHeap,
20    fmt, mem,
21    time::Instant,
22};
23use vek::*;
24
25pub type Alt = f64;
26pub type Compute = f64;
27pub type Computex8 = [Compute; 8];
28
29/* code used by sharp in future
30/// Compute the water flux at all chunks, given a list of chunk indices sorted
31/// by increasing height.
32pub fn get_drainage(
33    map_size_lg: MapSizeLg,
34    newh: &[u32],
35    downhill: &[isize],
36    _boundary_len: usize,
37) -> Box<[f32]> {
38    // FIXME: Make the below work.  For now, we just use constant flux.
39    // Initially, flux is determined by rainfall.  We currently treat this as the
40    // same as humidity, so we just use humidity as a proxy.  The total flux
41    // across the whole map is normalize to 1.0, and we expect the average flux
42    // to be 0.5.  To figure out how far from normal any given chunk is, we use
43    // its logit. NOTE: If there are no non-boundary chunks, we just set
44    // base_flux to 1.0; this should still work fine because in that case
45    // there's no erosion anyway. let base_flux = 1.0 / ((map_size_lg.chunks_len())
46    // as f32);
47    let base_flux = 1.0;
48    let mut flux = vec![base_flux; map_size_lg.chunks_len()].into_boxed_slice();
49    newh.iter().rev().for_each(|&chunk_idx| {
50        let chunk_idx = chunk_idx as usize;
51        let downhill_idx = downhill[chunk_idx];
52        if downhill_idx >= 0 {
53            flux[downhill_idx as usize] += flux[chunk_idx];
54        }
55    });
56    flux
57}
58*/
59
60/// Compute the water flux at all chunks for multiple receivers, given a list of
61/// chunk indices sorted by increasing height and weights for each receiver.
62pub fn get_multi_drainage(
63    map_size_lg: MapSizeLg,
64    mstack: &[u32],
65    mrec: &[u8],
66    mwrec: &[Computex8],
67    _boundary_len: usize,
68) -> Box<[Compute]> {
69    // FIXME: Make the below work.  For now, we just use constant flux.
70    // Initially, flux is determined by rainfall.  We currently treat this as the
71    // same as humidity, so we just use humidity as a proxy.  The total flux
72    // across the whole map is normalize to 1.0, and we expect the average flux
73    // to be 0.5.  To figure out how far from normal any given chunk is, we use
74    // its logit. NOTE: If there are no non-boundary chunks, we just set
75    // base_flux to 1.0; this should still work fine because in that case
76    // there's no erosion anyway.
77    let base_area = 1.0;
78    let mut area = vec![base_area; map_size_lg.chunks_len()].into_boxed_slice();
79    mstack.iter().for_each(|&ij| {
80        let ij = ij as usize;
81        let wrec_ij = &mwrec[ij];
82        let area_ij = area[ij];
83        mrec_downhill(map_size_lg, mrec, ij).for_each(|(k, ijr)| {
84            area[ijr] += area_ij * wrec_ij[k];
85        });
86    });
87    area
88    /*
89      a=dx*dy*precip
90      do ij=1,nn
91        ijk=mstack(ij)
92        do k =1,mnrec(ijk)
93          a(mrec(k,ijk))=a(mrec(k,ijk))+a(ijk)*mwrec(k,ijk)
94        enddo
95      enddo
96    */
97}
98
99/// Kind of water on this tile.
100#[derive(Clone, Copy, Debug, PartialEq)]
101pub enum RiverKind {
102    Ocean,
103    Lake {
104        /// In addition to a downhill node (pointing to, eventually, the bottom
105        /// of the lake), each lake also has a "pass" that identifies
106        /// the direction out of which water should flow from this lake
107        /// if it is minimally flooded.  While some lakes may be too full for
108        /// this to be the actual pass used by their enclosing lake, we
109        /// still use this as a way to make sure that lake connections
110        /// to rivers flow in the correct direction.
111        neighbor_pass_pos: Vec2<i32>,
112    },
113    /// River should be maximal.
114    River {
115        /// Dimensions of the river's cross-sectional area, as m × m (rivers are
116        /// approximated as an open rectangular prism in the direction
117        /// of the velocity vector).
118        cross_section: Vec2<f32>,
119    },
120}
121
122impl RiverKind {
123    pub fn is_ocean(&self) -> bool { matches!(*self, RiverKind::Ocean) }
124
125    pub fn is_river(&self) -> bool { matches!(*self, RiverKind::River { .. }) }
126
127    pub fn is_lake(&self) -> bool { matches!(*self, RiverKind::Lake { .. }) }
128}
129
130impl PartialOrd for RiverKind {
131    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
132        match (*self, *other) {
133            (RiverKind::Ocean, RiverKind::Ocean) => Some(Ordering::Equal),
134            (RiverKind::Ocean, _) => Some(Ordering::Less),
135            (_, RiverKind::Ocean) => Some(Ordering::Greater),
136            (RiverKind::Lake { .. }, RiverKind::Lake { .. }) => None,
137            (RiverKind::Lake { .. }, _) => Some(Ordering::Less),
138            (_, RiverKind::Lake { .. }) => Some(Ordering::Greater),
139            (RiverKind::River { .. }, RiverKind::River { .. }) => None,
140        }
141    }
142}
143
144/// From velocity and cross_section we can calculate the volumetric flow rate,
145/// as the cross-sectional area times the velocity.
146///
147/// TODO: we might choose to include a curve for the river, as long as it didn't
148/// allow it to cross more than one neighboring chunk away.  For now we defer
149/// this to rendering time.
150///
151/// NOTE: This structure is 57 (or more likely 64) bytes, which is kind of big.
152#[derive(Clone, Debug, Default)]
153pub struct RiverData {
154    /// A velocity vector (in m / minute, i.e. voxels / second from a game
155    /// perspective).
156    ///
157    /// TODO: To represent this in a better-packed way, use u8s instead (as
158    /// "f8s").
159    pub(crate) velocity: Vec3<f32>,
160    /// The computed derivative for the segment of river starting at this chunk
161    /// (and flowing downhill).  Should be 0 at endpoints.  For rivers with
162    /// more than one incoming segment, we weight the derivatives by flux
163    /// (cross-sectional area times velocity) which is correlated
164    /// with mass / second; treating the derivative as "velocity" with respect
165    /// to length along the river, we treat the weighted sum of incoming
166    /// splines as the "momentum", and can divide it by the total incoming
167    /// mass as a way to find the velocity of the center of mass.  We can
168    /// then use this derivative to find a "tangent" for the incoming river
169    /// segment at this point, and as the linear part of the interpolating
170    /// spline at this point.
171    ///
172    /// Note that we aren't going to have completely smooth curves here anyway,
173    /// so we will probably end up applying a dampening factor as well
174    /// (maybe based on the length?) to prevent extremely wild oscillations.
175    pub(crate) spline_derivative: Vec2<f32>,
176    /// If this chunk is part of a river, this should be true.  We can't just
177    /// compute this from the cross section because once a river becomes
178    /// visible, we want it to stay visible until it reaches its sink.
179    pub river_kind: Option<RiverKind>,
180    /// We also have a second record for recording any rivers in nearby chunks
181    /// that manage to intersect this chunk, though this is unlikely to
182    /// happen in current gameplay.  This is because river areas are allowed
183    /// to cross arbitrarily many chunk boundaries, if they are wide enough.
184    /// In such cases we may choose to render the rivers as particularly deep in
185    /// those places.
186    pub(crate) neighbor_rivers: Vec<u32>,
187}
188
189impl RiverData {
190    pub fn is_ocean(&self) -> bool {
191        self.river_kind
192            .as_ref()
193            .map(RiverKind::is_ocean)
194            .unwrap_or(false)
195    }
196
197    pub fn is_river(&self) -> bool {
198        self.river_kind
199            .as_ref()
200            .map(RiverKind::is_river)
201            .unwrap_or(false)
202    }
203
204    pub fn is_lake(&self) -> bool {
205        self.river_kind
206            .as_ref()
207            .map(RiverKind::is_lake)
208            .unwrap_or(false)
209    }
210
211    pub fn near_river(&self) -> bool { self.is_river() || !self.neighbor_rivers.is_empty() }
212
213    pub fn near_water(&self) -> bool { self.near_river() || self.is_lake() || self.is_ocean() }
214}
215
216/// Draw rivers and assign them heights, widths, and velocities.  Take some
217/// liberties with the constant factors etc. in order to make it more likely
218/// that we draw rivers at all.
219pub fn get_rivers<F: fmt::Debug + Float + Into<f64>, G: Float + Into<f64>>(
220    map_size_lg: MapSizeLg,
221    continent_scale_hack: f64,
222    newh: &[u32],
223    water_alt: &[F],
224    downhill: &[isize],
225    indirection: &[i32],
226    drainage: &[G],
227) -> Box<[RiverData]> {
228    // For continuity-preserving quadratic spline interpolation, we (appear to) need
229    // to build up the derivatives from the top down.  Fortunately this
230    // computation seems tractable.
231
232    let mut rivers = vec![RiverData::default(); map_size_lg.chunks_len()].into_boxed_slice();
233    let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
234    // (Roughly) area of a chunk, times minutes per second.
235    // NOTE: Clearly, this should "actually" be 1/60 (or maybe 1/64, if you want to
236    // retain powers of 2).
237    //
238    // But since we want rivers to form more often than they do in real life, we use
239    // this as a way to control the frequency of river formation.  As grid_scale
240    // increases, mins_per_sec should decrease, until it hits 1 / 60 or 1/ 64.
241    // For example, if grid_scale is multiplied by 4, mins_per_sec should be
242    // multiplied by 1 / (4.0 * 4.0).
243    let mins_per_sec = 1.0 / (continent_scale_hack * continent_scale_hack)/*1.0 / 16.0*//*1.0 / 64.0*/;
244    let chunk_area_factor = neighbor_coef.x * neighbor_coef.y * mins_per_sec;
245    // NOTE: This technically makes us discontinuous, so we should be cautious about
246    // using this.
247    let derivative_divisor = 1.0;
248    newh.iter().rev().for_each(|&chunk_idx| {
249        let chunk_idx = chunk_idx as usize;
250        let downhill_idx = downhill[chunk_idx];
251        if downhill_idx < 0 {
252            // We are in the ocean.
253            debug_assert!(downhill_idx == -2);
254            rivers[chunk_idx].river_kind = Some(RiverKind::Ocean);
255            return;
256        }
257        let downhill_idx = downhill_idx as usize;
258        let downhill_pos = uniform_idx_as_vec2(map_size_lg, downhill_idx);
259        let dxy = (downhill_pos - uniform_idx_as_vec2(map_size_lg, chunk_idx)).map(|e| e as f64);
260        let neighbor_dim = neighbor_coef * dxy;
261        // First, we calculate the river's volumetric flow rate.
262        let chunk_drainage = drainage[chunk_idx].into();
263        // Volumetric flow rate is just the total drainage area to this chunk, times
264        // rainfall height per chunk per minute, times minutes per second
265        // (needed in order to use this as a m³ volume).
266        // TODO: consider having different rainfall rates (and including this
267        // information in the computation of drainage).
268        let volumetric_flow_rate =
269            chunk_drainage * chunk_area_factor * CONFIG.rainfall_chunk_rate as f64;
270        let downhill_drainage = drainage[downhill_idx].into();
271
272        // We know the drainage to the downhill node is just chunk_drainage - 1.0 (the
273        // amount of rainfall this chunk is said to get), so we don't need to
274        // explicitly remember the incoming mass.  How do we find a slope for
275        // endpoints where there is no incoming data? Currently, we just assume
276        // it's set to 0.0. TODO: Fix this when we add differing amounts of
277        // rainfall.
278        let incoming_drainage = downhill_drainage - 1.0;
279        let get_river_spline_derivative =
280            |neighbor_dim: Vec2<f64>, spline_derivative: Vec2<f32>| {
281                // "Velocity of center of mass" of splines of incoming flows.
282                let river_prev_slope = spline_derivative.map(|e| e as f64);
283                // NOTE: We need to make sure the slope doesn't get *too* crazy.
284                // ((dpx - cx) - 4 * MAX).abs() = bx
285                // NOTE: This will fail if the distance between chunks in any direction
286                // is exactly TerrainChunkSize::RECT * 4.0, but hopefully this should not be
287                // possible. NOTE: This isn't measuring actual distance, you can
288                // go farther on diagonals.
289                let max_deriv = neighbor_dim - neighbor_coef * 2.0 * 2.0.sqrt();
290                let extra_divisor = river_prev_slope
291                    .map2(max_deriv, |e, f| (e / f).abs())
292                    .reduce_partial_max();
293                // Set up the river's spline derivative.  For each incoming river at pos with
294                // river_spline_derivative bx, we can compute our interpolated slope as:
295                //   d_x = 2 * (chunk_pos - pos - bx) + bx
296                //       = 2 * (chunk_pos - pos) - bx
297                //
298                // which is exactly twice what was weighted by uphill nodes to get our
299                // river_spline_derivative in the first place.
300                //
301                // NOTE: this probably implies that the distance shouldn't be normalized, since
302                // the distances aren't actually equal between x and y... we'll
303                // see what happens.
304                (if extra_divisor > 1.0 {
305                    river_prev_slope / extra_divisor
306                } else {
307                    river_prev_slope
308                })
309                .map(|e| e as f32)
310            };
311
312        let river = &rivers[chunk_idx];
313        let river_spline_derivative =
314            get_river_spline_derivative(neighbor_dim, river.spline_derivative);
315
316        let indirection_idx = indirection[chunk_idx];
317        // Find the lake we are flowing into.
318        let lake_idx = if indirection_idx < 0 {
319            // If we're a lake bottom, our own indirection is negative.
320            let pass_idx = (-indirection_idx) as usize;
321            // NOTE: Must exist since this lake had a downhill in the first place.
322            let neighbor_pass_idx = downhill[pass_idx] as usize/*downhill_idx*/;
323            let lake_neighbor_pass = &mut rivers[neighbor_pass_idx];
324            // We definitely shouldn't have encountered this yet!
325            debug_assert!(lake_neighbor_pass.velocity == Vec3::zero());
326            // TODO: Rethink making the lake neighbor pass always a river or lake, no matter
327            // how much incoming water there is?  Sometimes it looks weird
328            // having a river emerge from a tiny pool.
329            lake_neighbor_pass.river_kind = Some(RiverKind::River {
330                cross_section: Vec2::default(),
331            });
332            chunk_idx
333        } else {
334            indirection_idx as usize
335        };
336
337        // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
338        // pushed towards the point identified by pass_idx).
339        let pass_idx = if downhill[lake_idx] < 0 {
340            // Flows into nothing, so this lake is its own pass.
341            lake_idx
342        } else {
343            (-indirection[lake_idx]) as usize
344        };
345
346        // Add our spline derivative to the downhill river (weighted by the chunk's
347        // drainage). NOTE: Don't add the spline derivative to the lake side of
348        // the pass for our own lake, because we don't want to preserve weird
349        // curvature from before we hit the lake in the outflowing river (this
350        // will not apply to one-chunk lakes, which are their own pass).
351        if pass_idx != downhill_idx {
352            // TODO: consider utilizing height difference component of flux as well;
353            // currently we just discard it in figuring out the spline's slope.
354            let downhill_river = &mut rivers[downhill_idx];
355            let weighted_flow = (neighbor_dim * 2.0 - river_spline_derivative.map(|e| e as f64))
356                / derivative_divisor
357                * chunk_drainage
358                / incoming_drainage;
359            downhill_river.spline_derivative += weighted_flow.map(|e| e as f32);
360        }
361
362        let neighbor_pass_idx = downhill[pass_idx];
363        // Find our own water height.
364        let chunk_water_alt = water_alt[chunk_idx];
365        if neighbor_pass_idx >= 0 {
366            // We may be a river.  But we're not sure yet, since we still could be
367            // underwater.  Check the lake height and see if our own water height is within
368            // ε of it.
369            let lake_water_alt = water_alt[lake_idx];
370            if chunk_water_alt == lake_water_alt {
371                // Not a river.
372                // Check whether we we are the lake side of the pass.
373                // NOTE: Safe because this is a lake.
374                let (neighbor_pass_pos, river_spline_derivative) = if pass_idx == chunk_idx {
375                    // This is a pass, so set our flow direction to point to the neighbor pass
376                    // rather than downhill.
377                    // NOTE: Safe because neighbor_pass_idx >= 0.
378                    (
379                        uniform_idx_as_vec2(map_size_lg, downhill_idx),
380                        river_spline_derivative,
381                    )
382                } else {
383                    // Try pointing towards the lake side of the pass.
384                    (
385                        uniform_idx_as_vec2(map_size_lg, pass_idx),
386                        river_spline_derivative,
387                    )
388                };
389                let lake = &mut rivers[chunk_idx];
390                lake.spline_derivative = river_spline_derivative;
391                lake.river_kind = Some(RiverKind::Lake {
392                    neighbor_pass_pos: neighbor_pass_pos
393                        * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
394                });
395                return;
396            }
397        // Otherwise, we must be a river.
398        } else {
399            // We are flowing into the ocean.
400            debug_assert!(neighbor_pass_idx == -2);
401            // But we are not the ocean, so we must be a river.
402        }
403        // Now, we know we are a river *candidate*.  We still don't know whether we are
404        // actually a river, though.  There are two ways for that to happen:
405        // (i) We are already a river.
406        // (ii) Using the Gauckler–Manning–Strickler formula for cross-sectional
407        //      average velocity of water, we establish that the river can be
408        //      "big enough" to appear on the Veloren map.
409        //
410        // This is very imprecise, of course, and (ii) may (and almost certainly will)
411        // change over time.
412        //
413        // In both cases, we preemptively set our child to be a river, to make sure we
414        // have an unbroken stream.  Also in both cases, we go to the effort of
415        // computing an effective water velocity vector and cross-sectional
416        // dimensions, as well as figuring out the derivative of our
417        // interpolating spline (since this percolates through the whole river
418        // network).
419        let downhill_water_alt = water_alt[downhill_idx];
420        let neighbor_distance = neighbor_dim.magnitude();
421        let dz = (downhill_water_alt - chunk_water_alt).into();
422        let slope = dz.abs() / neighbor_distance;
423        if slope == 0.0 {
424            // This is not a river--how did this even happen?
425            let pass_idx = (-indirection_idx) as usize;
426            error!(
427                "Our chunk (and downhill, lake, pass, neighbor_pass): {:?} (to {:?}, in {:?} via \
428                 {:?} to {:?}), chunk water alt: {:?}, lake water alt: {:?}",
429                uniform_idx_as_vec2(map_size_lg, chunk_idx),
430                uniform_idx_as_vec2(map_size_lg, downhill_idx),
431                uniform_idx_as_vec2(map_size_lg, lake_idx),
432                uniform_idx_as_vec2(map_size_lg, pass_idx),
433                if neighbor_pass_idx >= 0 {
434                    Some(uniform_idx_as_vec2(map_size_lg, neighbor_pass_idx as usize))
435                } else {
436                    None
437                },
438                water_alt[chunk_idx],
439                water_alt[lake_idx]
440            );
441            panic!("Should this happen at all?");
442        }
443        let slope_sqrt = slope.sqrt();
444        // Now, we compute a quantity that is proportional to the velocity of the chunk,
445        // derived from the Manning formula, equal to
446        // volumetric_flow_rate / slope_sqrt * CONFIG.river_roughness.
447        let almost_velocity = volumetric_flow_rate / slope_sqrt * CONFIG.river_roughness as f64;
448        // From this, we can figure out the width of the chunk if we know the height.
449        // For now, we hardcode the height to 0.5, but it should almost
450        // certainly be much more complicated than this.
451        // let mut height = 0.5f32;
452        // We approximate the river as a rectangular prism.  Theoretically, we need to
453        // solve the following quintic equation to determine its width from its
454        // height:
455        //
456        // h^5 * w^5 = almost_velocity^3 * (w + 2 * h)^2.
457        //
458        // This is because one of the quantities in the Manning formula (the unknown) is
459        // R_h = (area of cross-section / h)^(2/3).
460        //
461        // Unfortunately, quintic equations do not in general have algebraic solutions,
462        // and it's not clear (to me anyway) that this one does in all cases.
463        //
464        // In practice, for high ratios of width to height, we can approximate the
465        // rectangular prism's perimeter as equal to its width, so R_h as equal
466        // to height.  This greatly simplifies the calculation.  For simplicity,
467        // we do this even for low ratios of width to height--I found that for
468        // most real rivers, at least big ones, this approximation is
469        // "good enough."  We don't need to be *that* realistic :P
470        //
471        // NOTE: Derived from a paper on estimating river width.
472        let mut width = 5.0
473            * (CONFIG.river_width_to_depth as f64
474                * (CONFIG.river_width_to_depth as f64 + 2.0).powf(2.0 / 3.0))
475            .powf(3.0 / 8.0)
476            * volumetric_flow_rate.powf(3.0 / 8.0)
477            * slope.powf(-3.0 / 16.0)
478            * (CONFIG.river_roughness as f64).powf(3.0 / 8.0);
479        width = width.max(0.0);
480
481        let mut height = if width == 0.0 {
482            CONFIG.river_min_height as f64
483        } else {
484            (almost_velocity / width).powf(3.0 / 5.0)
485        };
486
487        // We can now weight the river's drainage by its direction, which we use to help
488        // improve the slope of the downhill node.
489        let river_direction = Vec3::new(neighbor_dim.x, neighbor_dim.y, dz.signum() * dz);
490
491        // Now, we can check whether this is "really" a river.
492        // Currently, we just check that width and height are at least 0.5 and
493        // CONFIG.river_min_height.
494        let river = &rivers[chunk_idx];
495        let is_river = river.is_river() || width >= 0.5 && height >= CONFIG.river_min_height as f64;
496        let downhill_river = &mut rivers[downhill_idx];
497
498        if is_river {
499            // Provisionally make the downhill chunk a river as well.
500            downhill_river.river_kind = Some(RiverKind::River {
501                cross_section: Vec2::default(),
502            });
503
504            // Additionally, if the cross-sectional area for this river exceeds the max
505            // river width, the river is overflowing the two chunks adjacent to
506            // it, which we'd prefer to avoid since only its two immediate
507            // neighbors (orthogonal to the downhill direction) are guaranteed
508            // uphill of it. Solving this properly most likely requires
509            // modifying the erosion model to take channel width into account,
510            // which is a formidable task that likely requires rethinking the
511            // current grid-based erosion model (or at least, requires tracking some
512            // edges that aren't implied by the grid graph).  For now, we will solve this
513            // problem by making the river deeper when it hits the max width,
514            // until it consumes all the available energy in this part of the
515            // river.
516            let max_width = TerrainChunkSize::RECT_SIZE.x as f64 * CONFIG.river_max_width as f64;
517            if width > max_width {
518                width = max_width;
519                height = (almost_velocity / width).powf(3.0 / 5.0);
520            }
521        }
522        // Now we can compute the river's approximate velocity magnitude as well, as
523        let velocity_magnitude =
524            1.0 / CONFIG.river_roughness as f64 * height.powf(2.0 / 3.0) * slope_sqrt;
525
526        // Set up the river's cross-sectional area.
527        let cross_section = Vec2::new(width as f32, height as f32);
528        // Set up the river's velocity vector.
529        let mut velocity = river_direction;
530        velocity.normalize();
531        velocity *= velocity_magnitude;
532
533        let river = &mut rivers[chunk_idx];
534        // NOTE: Not trying to do this more cleverly because we want to keep the river's
535        // neighbors. TODO: Actually put something in the neighbors.
536        river.velocity = velocity.map(|e| e as f32);
537        river.spline_derivative = river_spline_derivative;
538        river.river_kind = if is_river {
539            Some(RiverKind::River { cross_section })
540        } else {
541            None
542        };
543    });
544    rivers
545}
546
547/// Precompute the maximum slope at all points.
548///
549/// TODO: See if allocating in advance is worthwhile.
550fn get_max_slope(
551    map_size_lg: MapSizeLg,
552    h: &[Alt],
553    rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
554    height_scale: impl Fn(usize) -> Alt + Sync,
555) -> Box<[f64]> {
556    let min_max_angle = (15.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
557    let max_max_angle = (60.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
558    let max_angle_range = max_max_angle - min_max_angle;
559    h.par_iter()
560        .enumerate()
561        .map(|(posi, &z)| {
562            let wposf = uniform_idx_as_vec2(map_size_lg, posi).map(|e| e as f64)
563                * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
564            let height_scale = height_scale(posi);
565            let wposz = z / height_scale;
566            // Normalized to be between 6 and and 54 degrees.
567            let div_factor = (2.0 * TerrainChunkSize::RECT_SIZE.x as f64) / 8.0;
568            let rock_strength = rock_strength_nz.get([wposf.x, wposf.y, wposz * div_factor]);
569            let rock_strength = rock_strength.clamp(-1.0, 1.0) * 0.5 + 0.5;
570            // Logistic regression.  Make sure x ∈ (0, 1).
571            let logit = |x: f64| x.ln() - (-x).ln_1p();
572            // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
573            let logistic_2_base = 3.0f64.sqrt() * std::f64::consts::FRAC_2_PI;
574            // Assumes μ = 0, σ = 1
575            let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
576
577            // We do log-odds against center, so that our log odds are 0 when x = 0.25,
578            // lower when x is lower, and higher when x is higher.
579            //
580            // (NOTE: below sea level, we invert it).
581            //
582            // TODO: Make all this stuff configurable... but honestly, it's so complicated
583            // that I'm not sure anyone would be able to usefully tweak them on
584            // a per-map basis?  Plus it's just a hacky heuristic anyway.
585            let center = 0.4;
586            let dmin = center - 0.05;
587            let dmax = center + 0.05;
588            let log_odds = |x: f64| logit(x) - logit(center);
589            let rock_strength = logistic_cdf(
590                1.0 * logit(rock_strength.clamp(1e-7, 1.0f64 - 1e-7))
591                    + 1.0
592                        * log_odds(
593                            (wposz / CONFIG.mountain_scale as f64)
594                                .abs()
595                                .clamp(dmin, dmax),
596                        ),
597            );
598            // NOTE: If you want to disable varying rock strength entirely, uncomment  this
599            // line. let max_slope = 3.0.sqrt() / 3.0;
600            rock_strength * max_angle_range + min_max_angle //max_slope
601        })
602        .collect::<Vec<_>>()
603        .into_boxed_slice()
604}
605
606// simd alternative
607#[cfg(not(feature = "simd"))]
608#[derive(Copy, Clone)]
609struct M32(u32);
610#[cfg(not(feature = "simd"))]
611impl M32 {
612    #[inline]
613    fn splat(x: bool) -> Self { if x { Self(u32::MAX) } else { Self(u32::MIN) } }
614
615    #[inline]
616    fn any(&self) -> bool { self.0 != 0 }
617}
618#[cfg(feature = "simd")]
619type M32 = std::simd::Mask<i32, 1>;
620
621/// Erode all chunks by amount.
622///
623/// Our equation is:
624///
625///   dh(p) / dt = uplift(p)−k * A(p)^m * slope(p)^n
626///
627///   where A(p) is the drainage area at p, m and n are constants
628///   (we choose m = 0.4 and n = 1), and k is a constant.  We choose
629///
630///   k = 2.244 * uplift.max() / (desired_max_height)
631///
632/// since this tends to produce mountains of max height desired_max_height;
633/// and we set desired_max_height = 1.0 to reflect limitations of mountain
634/// scale.
635///
636/// This algorithm does this in four steps:
637///
638/// 1. Sort the nodes in h by height (so the lowest node by altitude is first in
639///    the list, and the highest node by altitude is last).
640/// 2. Iterate through the list in *reverse.*  For each node, we compute its
641///    drainage area as the sum of the drainage areas of its "children" nodes
642///    (i.e. the nodes with directed edges to this node). To do this
643///    efficiently, we start with the "leaves" (the highest nodes), which have
644///    no neighbors higher than them, hence no directed edges to them. We add
645///    their area to themselves, and then to all neighbors that they flow into
646///    (their "ancestors" in the flow graph); currently, this just means the
647///    node immediately downhill of this node. As we go lower, we know that all
648///    our "children" already had their areas computed, which means that we can
649///    repeat the process in order to derive all the final areas.
650/// 3. Now, iterate through the list in *order.*  Whether we used the filling
651///    method to compute a "filled" version of each depression, or used the lake
652///    connection algorithm described in [1], each node is guaranteed to have
653///    zero or one drainage edges out, representing the direction of water flow
654///    for that node. For nodes i with zero drainage edges out (boundary nodes
655///    and lake bottoms) we set the slope to 0 (so the change in altitude is
656///    uplift(i)).
657///
658///    For nodes with at least one drainage edge out, we take
659///    advantage of the fact that we are computing new heights in order and
660///    rewrite our equation as (letting j = downhill[i], A[i] be the computed
661///    area of point i, p(i) be the x-y position of point i, flux(i) = k *
662///    A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1):
663///
664///    h[i](t + dt) = h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt)) / (1 +
665///    flux(i) * δt).
666///
667/// Since we compute heights in ascending order by height, and j is downhill
668/// from i, h[j] will    always be the *new* h[j](t + δt), while h[i] will still
669/// not have been computed yet, so we    only need to visit each node once.
670///
671/// Afterwards, we also apply a hillslope diffusion process using an ADI
672/// (alternating direction implicit) method:
673///
674/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
675///
676/// We also borrow the implementation for sediment transport from
677///
678/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90
679///
680/// The approximate equation for soil production function (predicting the rate
681/// at which bedrock turns into soil, increasing the distance between the
682/// basement and altitude) is taken from equation (11) from [2]. This (among
683/// numerous other sources) also includes at least one prediction that hillslope
684/// diffusion should be nonlinear, which we sort of attempt to approximate.
685///
686/// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani,
687///     Bedrich Benes, Eric Galin, et al..
688///     Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
689///     Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2),
690///     pp.165-175.     ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
691///
692/// [2] William E. Dietrich, Dino G. Bellugi, Leonard S. Sklar,
693///     Jonathan D. Stock
694///     Geomorphic Transport Laws for Predicting Landscape Form and Dynamics.
695///     Prediction in Geomorphology, Geophysical Monograph 135.
696///     Copyright 2003 by the American Geophysical Union
697///     10.1029/135GM09
698#[expect(clippy::too_many_arguments)]
699fn erode(
700    // Underlying map dimensions.
701    map_size_lg: MapSizeLg,
702    // Height above sea level of topsoil
703    h: &mut [Alt],
704    // Height above sea level of bedrock
705    b: &mut [Alt],
706    // Height above sea level of water
707    wh: &mut [Alt],
708    max_uplift: f32,
709    max_g: f32,
710    kdsed: f64,
711    _seed: &RandomField,
712    rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
713    uplift: impl Fn(usize) -> f32 + Sync,
714    n_f: impl Fn(usize) -> f32 + Sync,
715    m_f: impl Fn(usize) -> f32 + Sync,
716    kf: impl Fn(usize) -> f64 + Sync,
717    kd: impl Fn(usize) -> f64,
718    g: impl Fn(usize) -> f32 + Sync,
719    epsilon_0: impl Fn(usize) -> f32 + Sync,
720    alpha: impl Fn(usize) -> f32 + Sync,
721    is_ocean: impl Fn(usize) -> bool + Sync,
722    // scaling factors
723    height_scale: impl Fn(f32) -> Alt + Sync,
724    k_da_scale: impl Fn(f64) -> f64,
725    threadpool: &rayon::ThreadPool,
726) {
727    let compute_stats = true;
728    debug!("Done draining...");
729    // NOTE: To experimentally allow erosion to proceed below sea level, replace 0.0
730    // with -<Alt as Float>::infinity().
731    let min_erosion_height = 0.0; // -<Alt as Float>::infinity();
732
733    // NOTE: The value being divided by here sets the effective maximum uplift rate,
734    // as everything is scaled to it!
735    let dt = max_uplift as f64 / 1e-3;
736    debug!(?dt, "");
737    // Minimum sediment thickness before we treat erosion as sediment based.
738    let sediment_thickness = |_n| /*6.25e-5*/1.0e-4 * dt;
739    let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
740    let chunk_area = neighbor_coef.x * neighbor_coef.y;
741    let min_length = neighbor_coef.reduce_partial_min();
742    let max_stable = min_length * min_length / dt;
743
744    // Debris flow area coefficient (m^(-2q)).
745    let q = 0.2;
746    // NOTE: Set to 1.0 to make (assuming n = 1) the erosion equation linear during
747    // each stream power iteration.  This will result in significant speedups,
748    // at the cost of less interesting erosion behavior (linear vs. nonlinear).
749    let q_ = 1.5;
750    let k_da = 2.5 * k_da_scale(q);
751    let nx = usize::from(map_size_lg.chunks().x);
752    let ny = usize::from(map_size_lg.chunks().y);
753    let dx = TerrainChunkSize::RECT_SIZE.x as f64;
754    let dy = TerrainChunkSize::RECT_SIZE.y as f64;
755
756    // ```ignore
757    // ε₀ and α are part of the soil production approximate
758    // equation:
759    //
760    // -∂z_b / ∂t = ε₀ * e^(-αH)
761    //
762    // where
763    //    z_b is the elevation of the soil-bedrock interface (i.e. the basement),
764    //    ε₀ is the production rate of exposed bedrock (H = 0),
765    //    H is the soil thickness normal to the ground surface,
766    //    and α is a parameter (units of 1 / length).
767    //
768    // Note that normal depth at i, for us, will be interpreted as the soil depth vector,
769    //   sed_i = ((0, 0), h_i - b_i),
770    // projected onto the ground surface slope vector,
771    //   ground_surface_i = ((p_i - p_j), h_i - h_j),
772    // yielding the soil depth vector
773    //   H_i = sed_i - sed_i ⋅ ground_surface_i / (ground_surface_i ⋅ ground_surface_i) * ground_surface_i
774    //
775    //       = ((0, 0), h_i - b_i) -
776    //         (0 * ||p_i - p_j|| + (h_i - b_i) * (h_i - h_j)) / (||p_i - p_j||^2 + (h_i - h_j)^2)
777    //         * (p_i - p_j, h_i - h_j)
778    //       = ((0, 0), h_i - b_i) -
779    //         ((h_i - b_i) * (h_i - h_j)) / (||p_i - p_j||^2 + (h_i - h_j)^2)
780    //         * (p_i - p_j, h_i - h_j)
781    //       = (h_i - b_i) *
782    //         (((0, 0), 1) - (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2) * (p_i - p_j, h_i - h_j))
783    //   H_i_fact = (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2)
784    //   H_i = (h_i - b_i) * ((((0, 0), 1) - H_i_fact * (p_i - p_j, h_i - h_j)))
785    //       = (h_i - b_i) * (-H_i_fact * (p_i - p_j), 1 - H_i_fact * (h_i - h_j))
786    //   ||H_i|| = (h_i - b_i) * √(H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (h_i - h_j))^2)
787    //           = (h_i - b_i) * √(H_i_fact^2 * ||p_i - p_j||^2 + 1 - 2 * H_i_fact * (h_i - h_j) +
788    //                             H_i_fact^2 * (h_i - h_j)^2)
789    //           = (h_i - b_i) * √(H_i_fact^2 * (||p_i - p_j||^2 + (h_i - h_j)^2) +
790    //                             1 - 2 * H_i_fact * (h_i - h_j))
791    //           = (h_i - b_i) * √((h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2) * (||p_i - p_j||^2 + (h_i - h_j)^2) +
792    //                             1 - 2 * (h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2))
793    //           = (h_i - b_i) * √((h_i - h_j)^2 - 2(h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2) + 1)
794    //
795    // where j is i's receiver and ||p_i - p_j|| is the horizontal displacement between them.  The
796    // idea here is that we first compute the hypotenuse between the horizontal and vertical
797    // displacement of ground (getting the horizontal component of the triangle), and then this is
798    // taken as one of the non-hypotenuse sides of the triangle whose other non-hypotenuse side is
799    // the normal height H_i, while their square adds up to the vertical displacement (h_i - b_i).
800    // If h and b have different slopes, this may not work completely correctly, but this is
801    // probably fine as an approximation.
802    // ```
803
804    // Spatio-temporal variation in net precipitation rate ((m / year) / (m / year))
805    // (ratio of precipitation rate at chunk relative to mean precipitation rate
806    // p₀).
807    let p = 1.0;
808    // Dimensionless multiplier for stream power erosion constant when land becomes
809    // sediment.
810    let k_fs_mult_sed = 4.0;
811    // Dimensionless multiplier for G when land becomes sediment.
812    let g_fs_mult_sed = 1.0;
813    let ((dh, newh, maxh, mrec, mstack, mwrec, area), (mut max_slopes, h_t)) = threadpool.join(
814        || {
815            let mut dh = downhill(
816                map_size_lg,
817                |posi| h[posi],
818                |posi| is_ocean(posi) && h[posi] <= 0.0,
819            );
820            debug!("Computed downhill...");
821            let (boundary_len, _indirection, newh, maxh) =
822                get_lakes(map_size_lg, |posi| h[posi], &mut dh);
823            debug!("Got lakes...");
824            let (mrec, mstack, mwrec) = get_multi_rec(
825                map_size_lg,
826                |posi| h[posi],
827                &dh,
828                &newh,
829                wh,
830                nx,
831                ny,
832                dx as Compute,
833                dy as Compute,
834                maxh,
835                threadpool,
836            );
837            debug!("Got multiple receivers...");
838            // TODO: Figure out how to switch between single-receiver and multi-receiver
839            // drainage, as the former is much less computationally costly.
840            // let area = get_drainage(map_size_lg, &newh, &dh, boundary_len);
841            let area = get_multi_drainage(map_size_lg, &mstack, &mrec, &mwrec, boundary_len);
842            debug!("Got flux...");
843            (dh, newh, maxh, mrec, mstack, mwrec, area)
844        },
845        || {
846            threadpool.join(
847                || {
848                    let max_slope = get_max_slope(map_size_lg, h, rock_strength_nz, |posi| {
849                        height_scale(n_f(posi))
850                    });
851                    debug!("Got max slopes...");
852                    max_slope
853                },
854                || h.to_vec().into_boxed_slice(),
855            )
856        },
857    );
858
859    assert!(h.len() == dh.len() && dh.len() == area.len());
860
861    // max angle of slope depends on rock strength, which is computed from noise
862    // function. TODO: Make more principled.
863    let mid_slope = (30.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
864
865    type SimdType = f32;
866    type MaskType = M32;
867
868    // Precompute factors for Stream Power Law.
869    let czero = <SimdType as Zero>::zero();
870    let (k_fs_fact, k_df_fact) = threadpool.join(
871        || {
872            dh.par_iter()
873                .enumerate()
874                .map(|(posi, &posj)| {
875                    let mut k_tot = [czero; 8];
876                    if posj < 0 {
877                        // Egress with no outgoing flows, no stream power.
878                        k_tot
879                    } else {
880                        let old_b_i = b[posi];
881                        let sed = h_t[posi] - old_b_i;
882                        let n = n_f(posi);
883                        // Higher rock strength tends to lead to higher erodibility?
884                        let kd_factor = 1.0;
885                        let k_fs = kf(posi) / kd_factor;
886
887                        let k = if sed > sediment_thickness(n) {
888                            // Sediment
889                            k_fs_mult_sed * k_fs
890                        } else {
891                            // Bedrock
892                            k_fs
893                        } * dt;
894                        let n = n as f64;
895                        let m = m_f(posi) as f64;
896
897                        let mwrec_i = &mwrec[posi];
898                        mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
899                            let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
900                                - uniform_idx_as_vec2(map_size_lg, posj))
901                            .map(|e| e as f64);
902                            let neighbor_distance = (neighbor_coef * dxy).magnitude();
903                            let knew = (k * (p * chunk_area * (area[posi] * mwrec_i[kk])).powf(m)
904                                / neighbor_distance.powf(n))
905                                as SimdType;
906                            k_tot[kk] = knew;
907                        });
908                        k_tot
909                    }
910                })
911                .collect::<Vec<[SimdType; 8]>>()
912        },
913        || {
914            dh.par_iter()
915                .enumerate()
916                .map(|(posi, &posj)| {
917                    let mut k_tot = [czero; 8];
918                    let uplift_i = uplift(posi) as f64;
919                    debug_assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.0);
920                    if posj < 0 {
921                        // Egress with no outgoing flows, no stream power.
922                        k_tot
923                    } else {
924                        let area_i = area[posi];
925                        let max_slope = max_slopes[posi];
926                        let chunk_area_pow = chunk_area.powf(q);
927
928                        let old_b_i = b[posi];
929                        let sed = h_t[posi] - old_b_i;
930                        let n = n_f(posi);
931                        let g_i = if sed > sediment_thickness(n) {
932                            // Sediment
933                            (g_fs_mult_sed * g(posi)) as f64
934                        } else {
935                            // Bedrock
936                            g(posi) as f64
937                        };
938
939                        // Higher rock strength tends to lead to higher curvature?
940                        let kd_factor = (max_slope / mid_slope).powi(2);
941                        let k_da = k_da * kd_factor;
942
943                        let mwrec_i = &mwrec[posi];
944                        mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
945                            let mwrec_kk = mwrec_i[kk];
946
947                            #[rustfmt::skip]
948                            // Working equation:
949                            //   U = uplift per time
950                            //   D = sediment deposition per time
951                            //   E = fluvial erosion per time
952                            //   0 = U + D - E - k_df * (1 + k_da * (mrec_kk * A)^q) * (∂B/∂p)^(q_)
953                            //
954                            //   k_df = (U + D - E) / (1 + k_da * (mrec_kk * A)^q) / (∂B/∂p)^(q_)
955                            //
956                            // Want: ∂B/∂p = max slope at steady state, i.e.
957                            //     ∂B/∂p = max_slope
958                            // Then:
959                            //   k_df = (U + D - E) / (1 + k_da * (mrec_kk * A)^q) / max_slope^(q_)
960                            // Letting
961                            //   k = k_df * Δt
962                            // we get:
963                            //     k = (U + D - E)Δt / (1 + k_da * (mrec_kk * A)^q) / (ΔB)^(q_)
964                            //
965                            // Now ∂B/∂t under constant uplift, without debris flow (U + D - E), is
966                            //   U + D - E = U - E + G/(p̃A) * ∫_A((U - ∂h/∂t) * dA)
967                            //
968                            // Observing that at steady state ∂h/∂t should theoretically
969                            // be 0, we can simplify to:
970                            //   U + D = U + G/(p̃A) * ∫_A(U * dA)
971                            //
972                            // Upper bounding this at uplift = max_uplift/∂t for the whole prior
973                            // drainage area, and assuming we account for just mrec_kk of
974                            // the combined uplift and deposition, we get:
975                            //
976                            //   U + D ≤ mrec_kk * U + G/p̃ * max_uplift/∂t
977                            //   (U + D - E)Δt ≤ (mrec_kk * uplift_i + G/p̃ * mrec_kk * max_uplift - EΔt)
978                            //
979                            // therefore
980                            //   k * (1 + k_da * (mrec_kk * A)^q) * max_slope^(q_) ≤ (mrec_kk * (uplift_i + G/p̃ * max_uplift) - EΔt)
981                            // i.e.
982                            //   k ≤ (mrec_kk * (uplift_i + G/p̃ * max_uplift) - EΔt) / (1 + k_da * (mrec_kk * A)^q) / max_slope^q_
983                            //
984                            // (eliminating EΔt maintains the sign, but it's somewhat imprecise;
985                            //  we can address this later, e.g. by assigning a debris flow / fluvial erosion ratio).
986                            let chunk_neutral_area = 0.1e6; // 1 km^2 * (1000 m / km)^2 = 1e6 m^2
987                            let k = (mwrec_kk * (uplift_i + max_uplift as f64 * g_i / p))
988                                / (1.0 + k_da * (mwrec_kk * chunk_neutral_area).powf(q))
989                                / max_slope.powf(q_);
990                            let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
991                                - uniform_idx_as_vec2(map_size_lg, posj))
992                            .map(|e| e as f64);
993                            let neighbor_distance = (neighbor_coef * dxy).magnitude();
994
995                            let knew = (k
996                                * (1.0 + k_da * chunk_area_pow * (area_i * mwrec_kk).powf(q))
997                                / neighbor_distance.powf(q_))
998                                as SimdType;
999                            k_tot[kk] = knew;
1000                        });
1001                        k_tot
1002                    }
1003                })
1004                .collect::<Vec<[SimdType; 8]>>()
1005        },
1006    );
1007
1008    debug!("Computed stream power factors...");
1009
1010    let mut lake_water_volume: Box<[Compute]> =
1011        vec![0.0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1012    let mut elev: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1013    let mut h_p: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1014    let mut deltah: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1015
1016    // calculate the elevation / SPL, including sediment flux
1017    let tol1: Compute = 1.0e-4_f64 * (maxh as Compute + 1.0);
1018    let tol2: Compute = 1.0e-3_f64 * (max_uplift as Compute + 1.0);
1019    let tol = tol1.max(tol2);
1020    let mut err = 2.0 * tol;
1021
1022    // Some variables for tracking statistics, currently only for debugging
1023    // purposes.
1024    let mut minh = <Alt as Float>::infinity();
1025    let mut maxh = 0.0;
1026    let mut nland = 0usize;
1027    let mut ncorr = 0usize;
1028    let mut sums = 0.0;
1029    let mut sumh = 0.0;
1030    let mut sumsed = 0.0;
1031    let mut sumsed_land = 0.0;
1032    let mut ntherm = 0usize;
1033    // ln of product of actual slopes (only where actual is above critical).
1034    let mut prods_therm = 0.0;
1035    // ln of product of critical slopes (only where actual is above critical).
1036    let mut prodscrit_therm = 0.0;
1037    let avgz = |x, y: usize| if y == 0 { f64::INFINITY } else { x / y as f64 };
1038    let geomz = |x: f64, y: usize| {
1039        if y == 0 {
1040            f64::INFINITY
1041        } else {
1042            (x / y as f64).exp()
1043        }
1044    };
1045
1046    // Gauss-Seidel iteration
1047
1048    let mut lake_silt: Box<[Compute]> = vec![0.0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1049    let mut lake_sill = vec![-1isize; map_size_lg.chunks_len()].into_boxed_slice();
1050
1051    let mut n_gs_stream_power_law = 0;
1052    // NOTE: Increasing this can theoretically sometimes be necessary for
1053    // convergence, but in practice it is fairly unlikely that you should need
1054    // to do this (especially if you stick to g ∈ [0, 1]).
1055    let max_n_gs_stream_power_law = 99;
1056    let mut mstack_inv = vec![0usize; dh.len()];
1057    let mut h_t_stack = vec![Zero::zero(); dh.len()];
1058    let mut dh_stack = vec![0isize; dh.len()];
1059    let mut h_stack = vec![Zero::zero(); dh.len()];
1060    let mut b_stack = vec![Zero::zero(); dh.len()];
1061    let mut area_stack = vec![Zero::zero(); dh.len()];
1062    assert_eq!(mstack.len(), dh.len());
1063    assert_eq!(b.len(), dh.len());
1064    assert_eq!(h_t.len(), dh.len());
1065    let mstack_inv = &mut *mstack_inv;
1066    mstack.iter().enumerate().for_each(|(stacki, &posi)| {
1067        let posi = posi as usize;
1068        mstack_inv[posi] = stacki;
1069        dh_stack[stacki] = dh[posi];
1070        h_t_stack[stacki] = h_t[posi];
1071        h_stack[stacki] = h[posi];
1072        b_stack[stacki] = b[posi];
1073        area_stack[stacki] = area[posi];
1074    });
1075
1076    while err > tol && n_gs_stream_power_law < max_n_gs_stream_power_law {
1077        debug!("Stream power iteration #{:?}", n_gs_stream_power_law);
1078
1079        // Reset statistics in each loop.
1080        maxh = 0.0;
1081        minh = <Alt as Float>::infinity();
1082        nland = 0usize;
1083        ncorr = 0usize;
1084        sums = 0.0;
1085        sumh = 0.0;
1086        sumsed = 0.0;
1087        sumsed_land = 0.0;
1088        ntherm = 0usize;
1089        prods_therm = 0.0;
1090        prodscrit_therm = 0.0;
1091
1092        let start_time = Instant::now();
1093        // Keep track of how many iterations we've gone to test for convergence.
1094        n_gs_stream_power_law += 1;
1095
1096        threadpool.join(
1097            || {
1098                // guess/update the elevation at t+Δt (k)
1099                (&mut *h_p, &*h_stack)
1100                    .into_par_iter()
1101                    .for_each(|(h_p, h_)| {
1102                        *h_p = (*h_) as Compute;
1103                    });
1104            },
1105            || {
1106                // calculate erosion/deposition of sediment at each node
1107                (&*mstack, &mut *deltah, &*h_t_stack, &*h_stack)
1108                    .into_par_iter()
1109                    .for_each(|(&posi, deltah, &h_t_i, &h_i)| {
1110                        let posi = posi as usize;
1111                        let uplift_i = uplift(posi) as Alt;
1112                        let delta = (h_t_i + uplift_i - h_i) as Compute;
1113                        *deltah = delta;
1114                    });
1115            },
1116        );
1117        debug!(
1118            "(Done pre-computation, time={:?}ms).",
1119            start_time.elapsed().as_millis()
1120        );
1121        #[rustfmt::skip]
1122        // sum the erosion in stack order
1123        //
1124        // After:
1125        // deltah_i = Σ{j ∈ {i} ∪ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_i(t + Δt, k))
1126        let start_time = Instant::now();
1127        izip!(&*mstack, &*dh_stack, &h_t_stack, &*h_p)
1128            .enumerate()
1129            .for_each(|(stacki, (&posi, &posj, &h_t_i, &h_p_i))| {
1130                let posi = posi as usize;
1131                let deltah_i = deltah[stacki];
1132                if posj < 0 {
1133                    lake_silt[stacki] = deltah_i;
1134                } else {
1135                    let uplift_i = uplift(posi) as Alt;
1136                    let uphill_deltah_i = deltah_i - ((h_t_i + uplift_i) as Compute - h_p_i);
1137                    let lposi = lake_sill[stacki];
1138                    if lposi == stacki as isize {
1139                        if uphill_deltah_i <= 0.0 {
1140                            lake_silt[stacki] = 0.0;
1141                        } else {
1142                            lake_silt[stacki] = uphill_deltah_i;
1143                        }
1144                    }
1145                    let mwrec_i = &mwrec[posi];
1146                    mrec_downhill(map_size_lg, &mrec, posi).for_each(|(k, posj)| {
1147                        let stack_posj = mstack_inv[posj];
1148                        deltah[stack_posj] += deltah_i * mwrec_i[k];
1149                    });
1150                }
1151            });
1152        debug!(
1153            "(Done sediment transport computation, time={:?}ms).",
1154            start_time.elapsed().as_millis()
1155        );
1156        // ```ignore
1157        // do ij=nn,1,-1
1158        //   ijk=stack(ij)
1159        //   ijr=rec(ijk)
1160        //   if (ijr.ne.ijk) then
1161        //     dh(ijk)=dh(ijk)-(ht(ijk)-hp(ijk))
1162        //     if (lake_sill(ijk).eq.ijk) then
1163        //       if (dh(ijk).le.0.d0) then
1164        //         lake_sediment(ijk)=0.d0
1165        //       else
1166        //         lake_sediment(ijk)=dh(ijk)
1167        //       endif
1168        //     endif
1169        //     dh(ijk)=dh(ijk)+(ht(ijk)-hp(ijk))
1170        //     dh(ijr)=dh(ijr)+dh(ijk)
1171        //   else
1172        //     lake_sediment(ijk)=dh(ijk)
1173        //   endif
1174        // enddo
1175        // ```
1176
1177        let start_time = Instant::now();
1178        (
1179            &*mstack,
1180            &mut *elev,
1181            &*dh_stack,
1182            &*h_t_stack,
1183            &*area_stack,
1184            &*deltah,
1185            &*h_p,
1186            &*b_stack,
1187        )
1188            .into_par_iter()
1189            .for_each(
1190                |(&posi, elev, &dh_i, &h_t_i, &area_i, &deltah_i, &h_p_i, &b_i)| {
1191                    let posi = posi as usize;
1192
1193                    let uplift_i = uplift(posi) as Alt;
1194                    if dh_i < 0 {
1195                        *elev = (h_t_i + uplift_i) as Compute;
1196                    } else {
1197                        let old_h_after_uplift_i = (h_t_i + uplift_i) as Compute;
1198                        let area_i = area_i as Compute;
1199                        let uphill_silt_i = deltah_i - (old_h_after_uplift_i - h_p_i);
1200                        let old_b_i = b_i;
1201                        let sed = h_t_i - old_b_i;
1202                        let n = n_f(posi);
1203                        let g_i = if sed > sediment_thickness(n) {
1204                            (g_fs_mult_sed * g(posi)) as Compute
1205                        } else {
1206                            g(posi) as Compute
1207                        };
1208                        // Make sure deposition coefficient doesn't result in more deposition than
1209                        // there actually was material to deposit.  The
1210                        // current assumption is that as long as we
1211                        // are storing at most as much sediment as there actually was along the
1212                        // river, we are in the clear.
1213                        let g_i_ratio = g_i / (p * area_i);
1214                        // One side of nonlinear equation (23):
1215                        //
1216                        // h_i(t) + U_i * Δt + G / (p̃ * Ã_i) * Σ
1217                        // {j ∈ upstream_i(t)}(h_j(t, FINAL)
1218                        // + U_j * Δt - h_j(t + Δt, k))
1219                        //
1220                        // where
1221                        //
1222                        // Ã_i = A_i / (∆x∆y) = N_i,
1223                        // number of cells upstream of cell i.
1224                        *elev = old_h_after_uplift_i + uphill_silt_i * g_i_ratio;
1225                    }
1226                },
1227            );
1228        debug!(
1229            "(Done elevation estimation, time={:?}ms).",
1230            start_time.elapsed().as_millis()
1231        );
1232
1233        let start_time = Instant::now();
1234        // TODO: Consider taking advantage of multi-receiver flow here.
1235        // Iterate in ascending height order.
1236        let mut sum_err: Compute = 0.0_f64;
1237        izip!(&*mstack, &*elev, &*b_stack, &*h_t_stack, &*dh_stack, &*h_p)
1238            .enumerate()
1239            .rev()
1240            .for_each(|(stacki, (&posi, &elev_i, &b_i, &h_t_i, &dh_i, &h_p_i))| {
1241                let iteration_error = 0.0;
1242                let posi = posi as usize;
1243                let old_elev_i = elev_i;
1244                let old_b_i = b_i;
1245                let old_ht_i = h_t_i;
1246                let sed = old_ht_i - old_b_i;
1247
1248                let posj = dh_i;
1249                if posj < 0 {
1250                    if posj == -1 {
1251                        panic!("Disconnected lake!");
1252                    }
1253                    if h_t_i > 0.0 {
1254                        warn!("Ocean above zero?");
1255                    }
1256                    // Egress with no outgoing flows.
1257                    // wh for oceans is always at least min_erosion_height.
1258                    let uplift_i = uplift(posi) as Alt;
1259                    wh[posi] = FloatCore::max(min_erosion_height, h_t_i + uplift_i);
1260                    lake_sill[stacki] = posi as isize;
1261                    lake_water_volume[stacki] = 0.0;
1262                } else {
1263                    let posj = posj as usize;
1264
1265                    // Has an outgoing flow edge (posi, posj).
1266                    // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
1267                    // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 +
1268                    // flux(i) * δt). NOTE: posj has already been computed since
1269                    // it's downhill from us. Therefore, we can rely on wh being
1270                    // set to the water height for that node.
1271                    // let h_j = h[posj_stack] as f64;
1272                    let wh_j = wh[posj];
1273                    let old_h_i = h_stack[stacki];
1274                    let mut new_h_i = old_h_i;
1275
1276                    // Only perform erosion if we are above the water level of the previous node.
1277                    // NOTE: Can replace wh_j with h_j here (and a few other places) to allow
1278                    // erosion underwater, producing very different looking
1279                    // maps!
1280                    if old_elev_i > wh_j
1281                    /* h_j */
1282                    {
1283                        let dtherm = 0.0f64;
1284                        let h0 = old_elev_i + dtherm;
1285
1286                        // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 +
1287                        // ∂t * kp^mAi^m / ||pi - pj||)
1288                        let n = n_f(posi) as f64;
1289
1290                        // Fluvial erosion.
1291                        let k_df_fact = &k_df_fact[posi];
1292                        let k_fs_fact = &k_fs_fact[posi];
1293                        if (n - 1.0).abs() <= 1.0e-3 && (q_ - 1.0).abs() <= 1.0e-3 {
1294                            let mut f = h0;
1295                            let mut df = 1.0;
1296                            mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
1297                                let posj_stack = mstack_inv[posj];
1298                                let h_j = h_stack[posj_stack];
1299                                // This can happen in cases where receiver kk is neither uphill of
1300                                // nor downhill from posi's direct receiver.
1301                                // NOTE: Fastscape does h_t[posi] + uplift(posi) as f64 >= h_t[posj]
1302                                // + uplift(posj) as f64
1303                                // NOTE: We also considered using old_elev_i > wh[posj] here.
1304                                if old_elev_i > h_j {
1305                                    let elev_j = h_j;
1306                                    let fact = k_fs_fact[kk] as f64 + k_df_fact[kk] as f64;
1307                                    f += fact * elev_j;
1308                                    df += fact;
1309                                }
1310                            });
1311                            new_h_i = f / df;
1312                        } else {
1313                            // Local Newton-Raphson
1314                            // TODO: Work out how (if possible) to make this converge for tiny n.
1315                            let omega1 = 0.875f64 * n;
1316                            let omega2 = 0.875f64 / q_;
1317                            let omega = omega1.max(omega2);
1318                            let tolp = 1.0e-3;
1319                            let mut errp = 2.0 * tolp;
1320                            let mut rec_heights = [0.0; 8];
1321                            let mut mask = [MaskType::splat(false); 8];
1322                            mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
1323                                let posj_stack = mstack_inv[posj];
1324                                let h_j = h_stack[posj_stack];
1325                                // NOTE: Fastscape does h_t[posi] + uplift(posi) as f64 >= h_t[posj]
1326                                // + uplift(posj) as f64
1327                                // NOTE: We also considered using old_elev_i > wh[posj] here.
1328                                if old_elev_i > h_j {
1329                                    mask[kk] = MaskType::splat(true);
1330                                    rec_heights[kk] = h_j as SimdType;
1331                                }
1332                            });
1333                            while errp > tolp {
1334                                let mut f = new_h_i - h0;
1335                                let mut df = 1.0;
1336                                izip!(&mask, &rec_heights, k_fs_fact, k_df_fact).for_each(
1337                                    |(&mask_kk, &rec_heights_kk, &k_fs_fact_kk, &k_df_fact_kk)| {
1338                                        if mask_kk.any() {
1339                                            let h_j = rec_heights_kk;
1340                                            let elev_j = h_j;
1341                                            let dh =
1342                                                FloatCore::max(0.0, new_h_i as SimdType - elev_j);
1343                                            let powf = |a: SimdType, b| a.powf(b);
1344                                            let dh_fs_sample = k_fs_fact_kk as SimdType
1345                                                * powf(dh, n as SimdType - 1.0);
1346                                            let dh_df_sample = k_df_fact_kk as SimdType
1347                                                * powf(dh, q_ as SimdType - 1.0);
1348                                            // Want: h_i(t+Δt) = h0 - fact * (h_i(t+Δt) -
1349                                            // h_j(t+Δt))^n
1350                                            // Goal: h_i(t+Δt) - h0 + fact * (h_i(t+Δt) -
1351                                            // h_j(t+Δt))^n = 0
1352                                            f += ((dh_fs_sample + dh_df_sample) * dh) as f64;
1353                                            // ∂h_i(t+Δt)/∂n = 1 + fact * n * (h_i(t+Δt) -
1354                                            // h_j(t+Δt))^(n - 1)
1355                                            df += (n as SimdType * dh_fs_sample
1356                                                + q_ as SimdType * dh_df_sample)
1357                                                as f64;
1358                                        }
1359                                    },
1360                                );
1361                                // hn = h_i(t+Δt, k) - (h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) -
1362                                // h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)
1363                                let hn = new_h_i - f / df;
1364                                // errp = |(h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) -
1365                                // h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)|
1366                                errp = (hn - new_h_i).abs();
1367                                // h_i(t+∆t, k+1) = ...
1368                                new_h_i = new_h_i * (1.0 - omega) + hn * omega;
1369                            }
1370                            /* omega=0.875d0/n
1371                            tolp=1.d-3
1372                            errp=2.d0*tolp
1373                            h0=elev(ijk)
1374                            do while (errp.gt.tolp)
1375                              f=h(ijk)-h0
1376                              df=1.d0
1377                              if (ht(ijk).gt.ht(ijr)) then
1378                                fact = kfint(ijk)*dt*a(ijk)**m/length(ijk)**n
1379                                f=f+fact*max(0.d0,h(ijk)-h(ijr))**n
1380                                df=df+fact*n*max(0.d0,h(ijk)-h(ijr))**(n-1.d0)
1381                              endif
1382                              hn=h(ijk)-f/df
1383                              errp=abs(hn-h(ijk))
1384                              h(ijk)=h(ijk)*(1.d0-omega)+hn*omega
1385                            enddo */
1386                        }
1387
1388                        lake_sill[stacki] = posi as isize;
1389                        lake_water_volume[stacki] = 0.0;
1390
1391                        // If we dipped below the receiver's water level, set our height to the
1392                        // receiver's water level.
1393                        // NOTE: If we want erosion to proceed underwater, use h_j here instead of
1394                        // wh_j.
1395                        if new_h_i <= wh_j
1396                        /* h_j */
1397                        {
1398                            if compute_stats {
1399                                ncorr += 1;
1400                            }
1401                            // NOTE: Why wh_j?
1402                            // Because in the next round, if the old height is still wh_j or under,
1403                            // it will be set precisely equal to the
1404                            // estimated height, meaning it effectively
1405                            // "vanishes" and just deposits sediment to its receiver.
1406                            // (This is probably related to criteria for block Gauss-Seidel, etc.).
1407                            // NOTE: If we want erosion to proceed underwater, use h_j here instead
1408                            // of wh_j.
1409                            new_h_i = wh_j;
1410                        } else if compute_stats && new_h_i > 0.0 {
1411                            let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1412                                - uniform_idx_as_vec2(map_size_lg, posj))
1413                            .map(|e| e as f64);
1414                            let neighbor_distance = (neighbor_coef * dxy).magnitude();
1415                            let dz = (new_h_i - wh_j).max(0.0);
1416                            let mag_slope = dz / neighbor_distance;
1417
1418                            nland += 1;
1419                            sumsed_land += sed;
1420                            sumh += new_h_i;
1421                            sums += mag_slope;
1422                        }
1423                    } else {
1424                        new_h_i = old_elev_i;
1425                        let posj_stack = mstack_inv[posj];
1426                        let lposj = lake_sill[posj_stack];
1427                        lake_sill[stacki] = lposj;
1428                        if lposj >= 0 {
1429                            let lposj = lposj as usize;
1430                            lake_water_volume[lposj] += (wh_j - old_elev_i) as Compute;
1431                        }
1432                    }
1433                    // Set max_slope to this node's water height (max of receiver's water height and
1434                    // this node's height).
1435                    wh[posi] = wh_j.max(new_h_i) as Alt;
1436                    h_stack[stacki] = new_h_i as Alt;
1437                }
1438                if compute_stats {
1439                    sumsed += sed;
1440                    let h_i = h_stack[stacki];
1441                    if h_i > 0.0 {
1442                        minh = h_i.min(minh);
1443                    }
1444                    maxh = h_i.max(maxh);
1445                }
1446
1447                // Add sum of squares of errors.
1448                sum_err +=
1449                    (iteration_error + h_stack[stacki] as Compute - h_p_i as Compute).powi(2);
1450            });
1451        debug!(
1452            "(Done erosion computation, time={:?}ms)",
1453            start_time.elapsed().as_millis()
1454        );
1455
1456        err = (sum_err / mstack.len() as Compute).sqrt();
1457        debug!("(RMSE: {:?})", err);
1458        if max_g == 0.0 {
1459            err = 0.0;
1460        }
1461        if n_gs_stream_power_law == max_n_gs_stream_power_law {
1462            warn!(
1463                "Beware: Gauss-Seidel scheme not convergent: err={:?}, expected={:?}",
1464                err, tol
1465            );
1466        }
1467    }
1468
1469    (&*mstack_inv, &mut *h)
1470        .into_par_iter()
1471        .enumerate()
1472        .for_each(|(posi, (&stacki, h))| {
1473            assert_eq!(posi, mstack[stacki] as usize);
1474            *h = h_stack[stacki];
1475        });
1476
1477    // update the basement
1478    //
1479    // NOTE: Despite this not quite applying since basement order and height order
1480    // differ, we once again borrow the implicit FastScape stack order.  If this
1481    // becomes a problem we can easily compute a separate stack order just for
1482    // basement. TODO: Consider taking advantage of multi-receiver flow here.
1483    b.par_iter_mut()
1484        .zip_eq(h.par_iter())
1485        .enumerate()
1486        .for_each(|(posi, (b, &h_i))| {
1487            let old_b_i = *b;
1488            let uplift_i = uplift(posi) as Alt;
1489
1490            // First, add uplift...
1491            let mut new_b_i = old_b_i + uplift_i;
1492
1493            let posj = dh[posi];
1494            // Sediment height normal to bedrock.  NOTE: Currently we can actually have
1495            // sediment and bedrock slope at different heights, meaning there's
1496            // no uniform slope.  There are probably more correct ways to
1497            // account for this, such as averaging, integrating, or doing things
1498            // by mass / volume instead of height, but for now we use the time-honored
1499            // technique of ignoring the problem.
1500            let vertical_sed = (h_i - new_b_i).max(0.0);
1501            let h_normal = if posj < 0 {
1502                // Egress with no outgoing flows; for now, we assume this means normal and
1503                // vertical coincide.
1504                vertical_sed
1505            } else {
1506                let posj = posj as usize;
1507                let h_j = h[posj];
1508                let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1509                    - uniform_idx_as_vec2(map_size_lg, posj))
1510                .map(|e| e as f64);
1511                let neighbor_distance_squared = (neighbor_coef * dxy).magnitude_squared();
1512                let dh = h_i - h_j;
1513                // H_i_fact = (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2)
1514                let h_i_fact = dh / (neighbor_distance_squared + dh * dh);
1515                let h_i_vertical = 1.0 - h_i_fact * dh;
1516                // ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact *
1517                // (h_i - h_j))^2))
1518                vertical_sed
1519                    * (h_i_fact * h_i_fact * neighbor_distance_squared
1520                        + h_i_vertical * h_i_vertical)
1521                        .sqrt()
1522            };
1523            // Rate of sediment production: -∂z_b / ∂t = ε₀ * e^(-αH)
1524            let p_i = epsilon_0(posi) as f64 * dt * f64::exp(-alpha(posi) as f64 * h_normal);
1525
1526            new_b_i -= p_i as Alt;
1527
1528            // Clamp basement so it doesn't exceed height.
1529            new_b_i = new_b_i.min(h_i);
1530            *b = new_b_i;
1531        });
1532    debug!("Done updating basement and applying soil production...");
1533
1534    // update the height to reflect sediment flux.
1535    if max_g > 0.0 {
1536        //  If max_g = 0.0, lake_silt will be too high during the first iteration since
1537        // our  initial estimate for h is very poor; however, the elevation
1538        // estimate will have been  unaffected by g.
1539        (&mut *h, &*mstack_inv)
1540            .into_par_iter()
1541            .enumerate()
1542            .for_each(|(posi, (h, &stacki))| {
1543                let lposi = lake_sill[stacki];
1544                if lposi >= 0 {
1545                    let lposi = lposi as usize;
1546                    if lake_water_volume[lposi] > 0.0 {
1547                        // +max(0.d0,min(lake_sediment(lake_sill(ij)),
1548                        // lake_water_volume(lake_sill(ij))))/
1549                        // lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
1550                        *h += (FloatCore::max(0.0, lake_silt[stacki].min(lake_water_volume[lposi]))
1551                            / lake_water_volume[lposi]
1552                            * (wh[posi] - *h) as Compute) as Alt;
1553                    }
1554                }
1555            });
1556    }
1557    // do ij=1,nn
1558    //   if (lake_sill(ij).ne.0) then
1559    //     if (lake_water_volume(lake_sill(ij)).gt.0.d0) h(ij)=h(ij) &
1560    //     +max(0.d0,min(lake_sediment(lake_sill(ij)),
1561    // lake_water_volume(lake_sill(ij))))/ &
1562    //     lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
1563    //   endif
1564    // enddo
1565
1566    debug!(
1567        "Done applying stream power (max height: {:?}) (avg height: {:?}) (min height: {:?}) (avg \
1568         slope: {:?})\n        (above talus angle, geom. mean slope [actual/critical/ratio]: {:?} \
1569         / {:?} / {:?})\n        (old avg sediment thickness [all/land]: {:?} / {:?})\n        \
1570         (num land: {:?}) (num thermal: {:?}) (num corrected: {:?})",
1571        maxh,
1572        avgz(sumh, nland),
1573        minh,
1574        avgz(sums, nland),
1575        geomz(prods_therm, ntherm),
1576        geomz(prodscrit_therm, ntherm),
1577        geomz(prods_therm - prodscrit_therm, ntherm),
1578        avgz(sumsed, newh.len()),
1579        avgz(sumsed_land, nland),
1580        nland,
1581        ntherm,
1582        ncorr,
1583    );
1584
1585    // Apply thermal erosion.
1586    maxh = 0.0;
1587    minh = <Alt as Float>::infinity();
1588    sumh = 0.0;
1589    sums = 0.0;
1590    sumsed = 0.0;
1591    sumsed_land = 0.0;
1592    nland = 0usize;
1593    ncorr = 0usize;
1594    ntherm = 0usize;
1595    prods_therm = 0.0;
1596    prodscrit_therm = 0.0;
1597    newh.iter().for_each(|&posi| {
1598        let posi = posi as usize;
1599        let old_h_i = h[posi];
1600        let old_b_i = b[posi];
1601        let sed = old_h_i - old_b_i;
1602
1603        let max_slope = max_slopes[posi];
1604        let n = n_f(posi);
1605        max_slopes[posi] = if sed > sediment_thickness(n) && kdsed > 0.0 {
1606            // Sediment
1607            kdsed
1608        } else {
1609            // Bedrock
1610            kd(posi)
1611        };
1612
1613        let posj = dh[posi];
1614        if posj < 0 {
1615            // Egress with no outgoing flows.
1616            if posj == -1 {
1617                panic!("Disconnected lake!");
1618            }
1619            // wh for oceans is always at least min_erosion_height.
1620            wh[posi] = FloatCore::max(min_erosion_height, old_h_i as Alt);
1621        } else {
1622            let posj = posj as usize;
1623            // Find the water height for this chunk's receiver; we only apply thermal
1624            // erosion for chunks above water.
1625            let wh_j = wh[posj];
1626            let mut new_h_i = old_h_i;
1627            if wh_j < old_h_i {
1628                // NOTE: Currently assuming that talus angle is not eroded once the substance is
1629                // totally submerged in water, and that talus angle if part of the substance is
1630                // in water is 0 (or the same as the dry part, if this is set to wh_j), but
1631                // actually that's probably not true.
1632                let old_h_j = h[posj];
1633                let h_j = old_h_j;
1634                // Test the slope.
1635                // Hacky version of thermal erosion: only consider lowest neighbor, don't
1636                // redistribute uplift to other neighbors.
1637                let (posk, h_k) = (posj, h_j);
1638                let (posk, h_k) = if h_k < h_j { (posk, h_k) } else { (posj, h_j) };
1639                let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1640                    - uniform_idx_as_vec2(map_size_lg, posk))
1641                .map(|e| e as f64);
1642                let neighbor_distance = (neighbor_coef * dxy).magnitude();
1643                let dz = (new_h_i - h_k).max(0.0);
1644                let mag_slope = dz / neighbor_distance;
1645                if mag_slope >= max_slope {
1646                    let dtherm = 0.0;
1647                    new_h_i -= dtherm;
1648                    if new_h_i <= wh_j {
1649                        if compute_stats {
1650                            ncorr += 1;
1651                        }
1652                    } else if compute_stats && new_h_i > 0.0 {
1653                        let dz = (new_h_i - h_j).max(0.0);
1654                        let slope = dz / neighbor_distance;
1655                        sums += slope;
1656                        nland += 1;
1657                        sumh += new_h_i;
1658                        sumsed_land += sed;
1659                    }
1660                    if compute_stats {
1661                        ntherm += 1;
1662                        prodscrit_therm += max_slope.ln();
1663                        prods_therm += mag_slope.ln();
1664                    }
1665                } else {
1666                    // Poorly emulating nonlinear hillslope transport as described by
1667                    // http://eps.berkeley.edu/~bill/papers/112.pdf.
1668                    // sqrt(3)/3*32*32/(128000/2)
1669                    // Also Perron-2011-Journal_of_Geophysical_Research__Earth_Surface.pdf
1670                    let slope_ratio = (mag_slope / max_slope).powi(2);
1671                    let slope_nonlinear_factor =
1672                        slope_ratio * ((3.0 - slope_ratio) / (1.0 - slope_ratio).powi(2));
1673                    max_slopes[posi] += (max_slopes[posi] * slope_nonlinear_factor).min(max_stable);
1674                    if compute_stats && new_h_i > 0.0 {
1675                        sums += mag_slope;
1676                        // Just use the computed rate.
1677                        nland += 1;
1678                        sumh += new_h_i;
1679                        sumsed_land += sed;
1680                    }
1681                }
1682            }
1683            // Set wh to this node's water height (max of receiver's water height and
1684            // this node's height).
1685            wh[posi] = wh_j.max(new_h_i) as Alt;
1686        }
1687
1688        if compute_stats {
1689            sumsed += sed;
1690            let h_i = h[posi];
1691            if h_i > 0.0 {
1692                minh = h_i.min(minh);
1693            }
1694            maxh = h_i.max(maxh);
1695        }
1696    });
1697    debug!(
1698        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (min height: {:?}) \
1699         (avg slope: {:?})\n        (above talus angle, geom. mean slope [actual/critical/ratio]: \
1700         {:?} / {:?} / {:?})\n        (avg sediment thickness [all/land]: {:?} / {:?})\n        \
1701         (num land: {:?}) (num thermal: {:?}) (num corrected: {:?})",
1702        maxh,
1703        avgz(sumh, nland),
1704        minh,
1705        avgz(sums, nland),
1706        geomz(prods_therm, ntherm),
1707        geomz(prodscrit_therm, ntherm),
1708        geomz(prods_therm - prodscrit_therm, ntherm),
1709        avgz(sumsed, newh.len()),
1710        avgz(sumsed_land, nland),
1711        nland,
1712        ntherm,
1713        ncorr,
1714    );
1715
1716    // Apply hillslope diffusion.
1717    diffusion(
1718        nx,
1719        ny,
1720        nx as f64 * dx,
1721        ny as f64 * dy,
1722        dt,
1723        (),
1724        h,
1725        b,
1726        |posi| max_slopes[posi],
1727        -1.0,
1728    );
1729    debug!("Done applying diffusion.");
1730    debug!("Done eroding.");
1731}
1732
1733/// The Planchon-Darboux algorithm for extracting drainage networks.
1734///
1735/// http://horizon.documentation.ird.fr/exl-doc/pleins_textes/pleins_textes_7/sous_copyright/010031925.pdf
1736///
1737/// See https://github.com/mewo2/terrain/blob/master/terrain.js
1738pub(crate) fn fill_sinks<F: Float + Send + Sync>(
1739    map_size_lg: MapSizeLg,
1740    h: impl Fn(usize) -> F + Sync,
1741    is_ocean: impl Fn(usize) -> bool + Sync,
1742) -> Box<[F]> {
1743    // NOTE: We are using the "exact" version of depression-filling, which is slower
1744    // but doesn't change altitudes.
1745    let epsilon = F::zero();
1746    let infinity = F::infinity();
1747    let range = 0..map_size_lg.chunks_len();
1748    let oldh = range
1749        .into_par_iter()
1750        .map(&h)
1751        .collect::<Vec<_>>()
1752        .into_boxed_slice();
1753    let mut newh = oldh
1754        .par_iter()
1755        .enumerate()
1756        .map(|(posi, &h)| {
1757            let is_near_edge = is_ocean(posi);
1758            if is_near_edge {
1759                debug_assert!(h <= F::zero());
1760                h
1761            } else {
1762                infinity
1763            }
1764        })
1765        .collect::<Vec<_>>()
1766        .into_boxed_slice();
1767
1768    loop {
1769        let mut changed = false;
1770        (0..newh.len()).for_each(|posi| {
1771            let nh = newh[posi];
1772            let oh = oldh[posi];
1773            if nh == oh {
1774                return;
1775            }
1776            for nposi in neighbors(map_size_lg, posi) {
1777                let onbh = newh[nposi];
1778                let nbh = onbh + epsilon;
1779                // If there is even one path downhill from this node's original height, fix
1780                // the node's new height to be equal to its original height, and break out of
1781                // the loop.
1782                if oh >= nbh {
1783                    newh[posi] = oh;
1784                    changed = true;
1785                    break;
1786                }
1787                // Otherwise, we know this node's original height is below the new height of all
1788                // of its neighbors.  Then, we try to choose the minimum new
1789                // height among all this node's neighbors that is (plus a
1790                // constant epsilon) below this node's new height.
1791                //
1792                // (If there is no such node, then the node's new height must be (minus a
1793                // constant epsilon) lower than the new height of every
1794                // neighbor, but above its original height.  But this can't be
1795                // true for *all* nodes, because if this is true for any
1796                // node, it is not true of any of its neighbors.  So all neighbors must either
1797                // be their original heights, or we will have another iteration
1798                // of the loop (one of its neighbors was changed to its minimum
1799                // neighbor).  In the second case, in the next round, all
1800                // neighbor heights will be at most nh + epsilon).
1801                if nh > nbh && nbh > oh {
1802                    newh[posi] = nbh;
1803                    changed = true;
1804                }
1805            }
1806        });
1807        if !changed {
1808            return newh;
1809        }
1810    }
1811}
1812
1813/// Algorithm for finding and connecting lakes.  Assumes newh and downhill have
1814/// already been computed.  When a lake's value is negative, it is its own lake
1815/// root, and when it is 0, it is on the boundary of Ω.
1816///
1817/// Returns a 4-tuple containing:
1818/// - The first indirection vector (associating chunk indices with their lake's
1819///   root node).
1820/// - A list of chunks on the boundary (non-lake egress points).
1821/// - The second indirection vector (associating chunk indices with their lake's
1822///   adjacency list).
1823/// - The adjacency list (stored in a single vector), indexed by the second
1824///   indirection vector.
1825pub fn get_lakes<F: FloatCore>(
1826    map_size_lg: MapSizeLg,
1827    h: impl Fn(usize) -> F,
1828    downhill: &mut [isize],
1829) -> (usize, Box<[i32]>, Box<[u32]>, F) {
1830    // Associates each lake index with its root node (the deepest one in the lake),
1831    // and a list of adjacent lakes.  The list of adjacent lakes includes the
1832    // lake index of the adjacent lake, and a node index in the adjacent lake
1833    // which has a neighbor in this lake.  The particular neighbor should be the
1834    // one that generates the minimum "pass height" encountered so far, i.e. the
1835    // chosen pair should minimize the maximum of the heights of the nodes in the
1836    // pair.
1837
1838    // We start by taking steps to allocate an indirection vector to use for storing
1839    // lake indices. Initially, each entry in this vector will contain 0.  We
1840    // iterate in ascending order through the sorted newh.  If the node has
1841    // downhill == -2, it is a boundary node Ω and we store it in the boundary
1842    // vector.  If the node has downhill == -1, it is a fresh lake, and we store 0
1843    // in it.  If the node has non-negative downhill, we use the downhill index
1844    // to find the next node down; if the downhill node has a lake entry < 0,
1845    // then downhill is a lake and its entry can be negated to find an
1846    // (over)estimate of the number of entries it needs.  If the downhill
1847    // node has a non-negative entry, then its entry is the lake index for this
1848    // node, so we should access that entry and increment it, then set our own
1849    // entry to it.
1850    let mut boundary = Vec::with_capacity(downhill.len());
1851    let mut indirection = vec![/*-1i32*/0i32; map_size_lg.chunks_len()].into_boxed_slice();
1852
1853    let mut newh = Vec::with_capacity(downhill.len());
1854
1855    // Now, we know that the sum of all the indirection nodes will be the same as
1856    // the number of nodes.  We can allocate a *single* vector with 8 * nodes
1857    // entries, to be used as storage space, and augment our indirection vector
1858    // with the starting index, resulting in a vector of slices.  As we go, we
1859    // replace each lake entry with its index in the new indirection buffer,
1860    // allowing
1861    let mut lakes = vec![(-1, 0); /*(indirection.len() - boundary.len())*/indirection.len() * 8];
1862    let mut indirection_ = vec![0u32; indirection.len()];
1863    // First, find all the lakes.
1864    let mut lake_roots = Vec::with_capacity(downhill.len()); // Test
1865    (*downhill)
1866        .iter()
1867        .enumerate()
1868        .filter(|(_, &dh_idx)| dh_idx < 0)
1869        .for_each(|(chunk_idx, &dh)| {
1870            if dh == -2 {
1871                // On the boundary, add to the boundary vector.
1872                boundary.push(chunk_idx);
1873            // Still considered a lake root, though.
1874            } else if dh == -1 {
1875                lake_roots.push(chunk_idx);
1876            } else {
1877                panic!("Impossible.");
1878            }
1879            // Find all the nodes uphill from this lake.  Since there is only one outgoing
1880            // edge in the "downhill" graph, this is guaranteed never to visit a
1881            // node more than once.
1882            let start = newh.len();
1883            let indirection_idx = (start * 8) as u32;
1884            // New lake root
1885            newh.push(chunk_idx as u32);
1886            let mut cur = start;
1887            while cur < newh.len() {
1888                let node = newh[cur];
1889                uphill(map_size_lg, downhill, node as usize).for_each(|child| {
1890                    // lake_idx is the index of our lake root.
1891                    indirection[child] = chunk_idx as i32;
1892                    indirection_[child] = indirection_idx;
1893                    newh.push(child as u32);
1894                });
1895                cur += 1;
1896            }
1897            // Find the number of elements pushed.
1898            let length = (cur - start) * 8;
1899            // New lake root (lakes have indirection set to -length).
1900            indirection[chunk_idx] = -(length as i32);
1901            indirection_[chunk_idx] = indirection_idx;
1902        });
1903    assert_eq!(newh.len(), downhill.len());
1904
1905    debug!("Old lake roots: {:?}", lake_roots.len());
1906
1907    let newh = newh.into_boxed_slice();
1908    let mut maxh = -F::infinity();
1909    // Now, we know that the sum of all the indirection nodes will be the same as
1910    // the number of nodes.  We can allocate a *single* vector with 8 * nodes
1911    // entries, to be used as storage space, and augment our indirection vector
1912    // with the starting index, resulting in a vector of slices.  As we go, we
1913    // replace each lake entry with its index in the new indirection buffer,
1914    // allowing
1915    newh.iter().for_each(|&chunk_idx_| {
1916        let chunk_idx = chunk_idx_ as usize;
1917        let lake_idx_ = indirection_[chunk_idx];
1918        let lake_idx = lake_idx_ as usize;
1919        let height = h(chunk_idx_ as usize);
1920        // While we're here, compute the max elevation difference from zero among all
1921        // nodes.
1922        maxh = maxh.max(height.abs());
1923        // For every neighbor, check to see whether it is already set; if the neighbor
1924        // is set, its height is ≤ our height.  We should search through the
1925        // edge list for the neighbor's lake to see if there's an entry; if not,
1926        // we insert, and otherwise we get its height.  We do the same thing in
1927        // our own lake's entry list.  If the maximum of the heights we get out
1928        // from this process is greater than the maximum of this chunk and its
1929        // neighbor chunk, we switch to this new edge.
1930        neighbors(map_size_lg, chunk_idx).for_each(|neighbor_idx| {
1931            let neighbor_height = h(neighbor_idx);
1932            let neighbor_lake_idx_ = indirection_[neighbor_idx];
1933            let neighbor_lake_idx = neighbor_lake_idx_ as usize;
1934            if neighbor_lake_idx_ < lake_idx_ {
1935                // We found an adjacent node that is not on the boundary and has already
1936                // been processed, and also has a non-matching lake.  Therefore we can use
1937                // split_at_mut to get disjoint slices.
1938                let (lake, neighbor_lake) = {
1939                    // println!("Okay, {:?} < {:?}", neighbor_lake_idx, lake_idx);
1940                    let (neighbor_lake, lake) = lakes.split_at_mut(lake_idx);
1941                    (lake, &mut neighbor_lake[neighbor_lake_idx..])
1942                };
1943
1944                // We don't actually need to know the real length here, because we've reserved
1945                // enough spaces that we should always either find a -1 (available slot) or an
1946                // entry for this chunk.
1947                'outer: for pass in lake.iter_mut() {
1948                    if pass.0 == -1 {
1949                        // println!("One time, in my mind, one time... (neighbor lake={:?}
1950                        // lake={:?})", neighbor_lake_idx, lake_idx_);
1951                        *pass = (chunk_idx_ as i32, neighbor_idx as u32);
1952                        // Should never run out of -1s in the neighbor lake if we didn't find
1953                        // the neighbor lake in our lake.
1954                        *neighbor_lake
1955                            .iter_mut()
1956                            .find(|neighbor_pass| neighbor_pass.0 == -1)
1957                            .unwrap() = (neighbor_idx as i32, chunk_idx_);
1958                        // panic!("Should never happen; maybe didn't reserve enough space in
1959                        // lakes?")
1960                        break;
1961                    } else if indirection_[pass.1 as usize] == neighbor_lake_idx_ {
1962                        for neighbor_pass in neighbor_lake.iter_mut() {
1963                            // Should never run into -1 while looping here, since (i, j)
1964                            // and (j, i) should be added together.
1965                            if indirection_[neighbor_pass.1 as usize] == lake_idx_ {
1966                                let pass_height = h(neighbor_pass.1 as usize);
1967                                let neighbor_pass_height = h(pass.1 as usize);
1968                                if height.max(neighbor_height)
1969                                    < pass_height.max(neighbor_pass_height)
1970                                {
1971                                    *pass = (chunk_idx_ as i32, neighbor_idx as u32);
1972                                    *neighbor_pass = (neighbor_idx as i32, chunk_idx_);
1973                                }
1974                                break 'outer;
1975                            }
1976                        }
1977                        // Should always find a corresponding match in the neighbor lake if
1978                        // we found the neighbor lake in our lake.
1979                        let indirection_idx = indirection[chunk_idx];
1980                        let lake_chunk_idx = if indirection_idx >= 0 {
1981                            indirection_idx as usize
1982                        } else {
1983                            chunk_idx
1984                        };
1985                        let indirection_idx = indirection[neighbor_idx];
1986                        let neighbor_lake_chunk_idx = if indirection_idx >= 0 {
1987                            indirection_idx as usize
1988                        } else {
1989                            neighbor_idx
1990                        };
1991                        panic!(
1992                            "For edge {:?} between lakes {:?}, couldn't find partner for pass \
1993                             {:?}. Should never happen; maybe forgot to set both edges?",
1994                            (
1995                                (chunk_idx, uniform_idx_as_vec2(map_size_lg, chunk_idx)),
1996                                (neighbor_idx, uniform_idx_as_vec2(map_size_lg, neighbor_idx))
1997                            ),
1998                            (
1999                                (
2000                                    lake_chunk_idx,
2001                                    uniform_idx_as_vec2(map_size_lg, lake_chunk_idx),
2002                                    lake_idx_
2003                                ),
2004                                (
2005                                    neighbor_lake_chunk_idx,
2006                                    uniform_idx_as_vec2(map_size_lg, neighbor_lake_chunk_idx),
2007                                    neighbor_lake_idx_
2008                                )
2009                            ),
2010                            (
2011                                (pass.0, uniform_idx_as_vec2(map_size_lg, pass.0 as usize)),
2012                                (pass.1, uniform_idx_as_vec2(map_size_lg, pass.1 as usize))
2013                            ),
2014                        );
2015                    }
2016                }
2017            }
2018        });
2019    });
2020
2021    // Now it's time to calculate the lake connections graph T_L covering G_L.
2022    let mut candidates = BinaryHeap::with_capacity(indirection.len());
2023    // let mut pass_flows : Vec<i32> = vec![-1; indirection.len()];
2024
2025    // We start by going through each pass, deleting the ones that point out of
2026    // boundary nodes and adding ones that point into boundary nodes from
2027    // non-boundary nodes.
2028    lakes.iter_mut().for_each(|edge| {
2029        let edge: &mut (i32, u32) = edge;
2030        // Only consider valid elements.
2031        if edge.0 == -1 {
2032            return;
2033        }
2034        // Check to see if this edge points out *from* a boundary node.
2035        // Delete it if so.
2036        let from = edge.0 as usize;
2037        let indirection_idx = indirection[from];
2038        let lake_idx = if indirection_idx < 0 {
2039            from
2040        } else {
2041            indirection_idx as usize
2042        };
2043        if downhill[lake_idx] == -2 {
2044            edge.0 = -1;
2045            return;
2046        }
2047        // This edge is not pointing out from a boundary node.
2048        // Check to see if this edge points *to* a boundary node.
2049        // Add it to the candidate set if so.
2050        let to = edge.1 as usize;
2051        let indirection_idx = indirection[to];
2052        let lake_idx = if indirection_idx < 0 {
2053            to
2054        } else {
2055            indirection_idx as usize
2056        };
2057        if downhill[lake_idx] == -2 {
2058            // Find the pass height
2059            let pass = h(from).max(h(to));
2060            candidates.push(Reverse((
2061                NotNan::new(pass).unwrap(),
2062                (edge.0 as u32, edge.1),
2063            )));
2064        }
2065    });
2066
2067    let mut pass_flows_sorted: Vec<usize> = Vec::with_capacity(indirection.len());
2068
2069    // Now all passes pointing to the boundary are in candidates.
2070    // As long as there are still candidates, we continue...
2071    // NOTE: After a lake is added to the stream tree, the lake bottom's indirection
2072    // entry no longer negates its maximum number of passes, but the lake side
2073    // of the chosen pass.  As such, we should make sure not to rely on using it
2074    // this way afterwards. provides information about the number of candidate
2075    // passes in a lake.
2076    while let Some(Reverse((_, (chunk_idx, neighbor_idx)))) = candidates.pop() {
2077        // We have the smallest candidate.
2078        let lake_idx = indirection_[chunk_idx as usize] as usize;
2079        let indirection_idx = indirection[chunk_idx as usize];
2080        let lake_chunk_idx = if indirection_idx >= 0 {
2081            indirection_idx as usize
2082        } else {
2083            chunk_idx as usize
2084        };
2085        if downhill[lake_chunk_idx] >= 0 {
2086            // Candidate lake has already been connected.
2087            continue;
2088        }
2089        assert_eq!(downhill[lake_chunk_idx], -1);
2090        // Candidate lake has not yet been connected, and is the lowest candidate.
2091        // Delete all other outgoing edges.
2092        let max_len = -if indirection_idx < 0 {
2093            indirection_idx
2094        } else {
2095            indirection[indirection_idx as usize]
2096        } as usize;
2097        // Add this chunk to the tree.
2098        downhill[lake_chunk_idx] = neighbor_idx as isize;
2099        // Also set the indirection of the lake bottom to the negation of the
2100        // lake side of the chosen pass (chunk_idx).
2101        // NOTE: This can't overflow i32 because map_size_lg.chunks_len() should fit
2102        // in an i32.
2103        indirection[lake_chunk_idx] = -(chunk_idx as i32);
2104        // Add this edge to the sorted list.
2105        pass_flows_sorted.push(lake_chunk_idx);
2106        // pass_flows_sorted.push((chunk_idx as u32, neighbor_idx as u32));
2107        for edge in &mut lakes[lake_idx..lake_idx + max_len] {
2108            if *edge == (chunk_idx as i32, neighbor_idx) {
2109                // Skip deleting this edge.
2110                continue;
2111            }
2112            // Delete the old edge, and remember it.
2113            let edge = mem::replace(edge, (-1, 0));
2114            if edge.0 == -1 {
2115                // Don't fall off the end of the list.
2116                break;
2117            }
2118            // Don't add incoming pointers from already-handled lakes or boundary nodes.
2119            let indirection_idx = indirection[edge.1 as usize];
2120            let neighbor_lake_idx = if indirection_idx < 0 {
2121                edge.1 as usize
2122            } else {
2123                indirection_idx as usize
2124            };
2125            if downhill[neighbor_lake_idx] != -1 {
2126                continue;
2127            }
2128            // Find the pass height
2129            let pass = h(edge.0 as usize).max(h(edge.1 as usize));
2130            // Put the reverse edge in candidates, sorted by height, then chunk idx, and
2131            // finally neighbor idx.
2132            candidates.push(Reverse((
2133                NotNan::new(pass).unwrap(),
2134                (edge.1, edge.0 as u32),
2135            )));
2136        }
2137    }
2138    debug!("Total lakes: {:?}", pass_flows_sorted.len());
2139
2140    // Perform the bfs once again.
2141    #[derive(Clone, Copy, PartialEq)]
2142    enum Tag {
2143        UnParsed,
2144        InQueue,
2145        WithRcv,
2146    }
2147    let mut tag = vec![Tag::UnParsed; map_size_lg.chunks_len()];
2148    // TODO: Combine with adding to vector.
2149    let mut filling_queue = Vec::with_capacity(downhill.len());
2150
2151    let mut newh = Vec::with_capacity(downhill.len());
2152    (*boundary)
2153        .iter()
2154        .chain(pass_flows_sorted.iter())
2155        .for_each(|&chunk_idx| {
2156            // Find all the nodes uphill from this lake.  Since there is only one outgoing
2157            // edge in the "downhill" graph, this is guaranteed never to visit a
2158            // node more than once.
2159            let mut start = newh.len();
2160            // First, find the neighbor pass (assuming this is not the ocean).
2161            let neighbor_pass_idx = downhill[chunk_idx];
2162            let first_idx = if neighbor_pass_idx < 0 {
2163                // This is the ocean.
2164                newh.push(chunk_idx as u32);
2165                start += 1;
2166                chunk_idx
2167            } else {
2168                // This is a "real" lake.
2169                let neighbor_pass_idx = neighbor_pass_idx as usize;
2170                // Let's find our side of the pass.
2171                let pass_idx = -indirection[chunk_idx];
2172                // NOTE: Since only lakes are on the boundary, this should be a valid array
2173                // index.
2174                assert!(pass_idx >= 0);
2175                let pass_idx = pass_idx as usize;
2176                // Now, we should recompute flow paths so downhill nodes are contiguous.
2177
2178                /* // Carving strategy: reverse the path from the lake side of the pass to the
2179                // lake bottom, and also set the lake side of the pass's downhill to its
2180                // neighbor pass.
2181                let mut to_idx = neighbor_pass_idx;
2182                let mut from_idx = pass_idx;
2183                // NOTE: Since our side of the lake pass must be in the same basin as chunk_idx,
2184                // and chunk_idx is the basin bottom, we must reach it before we reach an ocean
2185                // node or other node with an invalid index.
2186                while from_idx != chunk_idx {
2187                    // Reverse this (from, to) edge by first replacing to_idx with from_idx,
2188                    // then replacing from_idx's downhill with the old to_idx, and finally
2189                    // replacing from_idx with from_idx's old downhill.
2190                    //
2191                    // println!("Reversing (lake={:?}): to={:?}, from={:?}, dh={:?}", chunk_idx, to_idx, from_idx, downhill[from_idx]);
2192                    from_idx = mem::replace(
2193                        &mut downhill[from_idx],
2194                        mem::replace(
2195                            &mut to_idx,
2196                            // NOTE: This cast should be valid since the node is either a path on the way
2197                            // to a lake bottom, or a lake bottom with an actual pass outwards.
2198                            from_idx
2199                        ) as isize,
2200                    ) as usize;
2201                }
2202                // Remember to set the actual lake's from_idx properly!
2203                downhill[from_idx] = to_idx as isize; */
2204
2205                // TODO: Enqueue onto newh simultaneously with filling (this could help prevent
2206                // needing a special tag just for checking if we are already enqueued).
2207                // Filling strategy: nodes are assigned paths based on cost.
2208                {
2209                    assert!(tag[pass_idx] == Tag::UnParsed);
2210
2211                    filling_queue.push(pass_idx);
2212                    tag[neighbor_pass_idx] = Tag::WithRcv;
2213                    tag[pass_idx] = Tag::InQueue;
2214
2215                    let outflow_coords = uniform_idx_as_vec2(map_size_lg, neighbor_pass_idx);
2216                    let elev = h(neighbor_pass_idx).max(h(pass_idx));
2217
2218                    while let Some(node) = filling_queue.pop() {
2219                        let coords = uniform_idx_as_vec2(map_size_lg, node);
2220
2221                        let mut rcv = -1;
2222                        let mut rcv_cost = -f64::INFINITY; /*f64::EPSILON;*/
2223                        let outflow_distance = (outflow_coords - coords).map(|e| e as f64);
2224
2225                        neighbors(map_size_lg, node).for_each(|ineighbor| {
2226                            if indirection[ineighbor] != chunk_idx as i32
2227                                && ineighbor != chunk_idx
2228                                && ineighbor != neighbor_pass_idx
2229                                || h(ineighbor) > elev
2230                            {
2231                                return;
2232                            }
2233                            let dxy = (uniform_idx_as_vec2(map_size_lg, ineighbor) - coords)
2234                                .map(|e| e as f64);
2235                            let neighbor_distance = /*neighbor_coef * */dxy;
2236                            let tag = &mut tag[ineighbor];
2237                            match *tag {
2238                                Tag::WithRcv => {
2239                                    // TODO: Remove outdated comment.
2240                                    //
2241                                    // vec_to_outflow ⋅ (vec_to_neighbor / |vec_to_neighbor|) =
2242                                    // ||vec_to_outflow||cos Θ
2243                                    //   where θ is the angle between vec_to_outflow and
2244                                    // vec_to_neighbor.
2245                                    //
2246                                    // Which is also the scalar component of vec_to_outflow in the
2247                                    // direction of vec_to_neighbor.
2248                                    let cost = outflow_distance
2249                                        .dot(neighbor_distance / neighbor_distance.magnitude());
2250                                    if cost > rcv_cost {
2251                                        rcv = ineighbor as isize;
2252                                        rcv_cost = cost;
2253                                    }
2254                                },
2255                                Tag::UnParsed => {
2256                                    filling_queue.push(ineighbor);
2257                                    *tag = Tag::InQueue;
2258                                },
2259                                _ => {},
2260                            }
2261                        });
2262                        assert_ne!(rcv, -1);
2263                        downhill[node] = rcv;
2264                        tag[node] = Tag::WithRcv;
2265                    }
2266                }
2267
2268                // Use our side of the pass as the initial node in the stack order.
2269                // TODO: Verify that this stack order will not "double reach" any lake chunks.
2270                neighbor_pass_idx
2271            };
2272            // New lake root
2273            let mut cur = start;
2274            let mut node = first_idx;
2275            loop {
2276                uphill(map_size_lg, downhill, node).for_each(|child| {
2277                    // lake_idx is the index of our lake root.
2278                    // Check to make sure child (flowing into us) is in the same lake.
2279                    if indirection[child] == chunk_idx as i32 || child == chunk_idx {
2280                        newh.push(child as u32);
2281                    }
2282                });
2283
2284                if cur == newh.len() {
2285                    break;
2286                }
2287                node = newh[cur] as usize;
2288                cur += 1;
2289            }
2290        });
2291    assert_eq!(newh.len(), downhill.len());
2292    (boundary.len(), indirection, newh.into_boxed_slice(), maxh)
2293}
2294
2295/// Iterate through set neighbors of multi-receiver flow.
2296pub fn mrec_downhill(
2297    map_size_lg: MapSizeLg,
2298    mrec: &[u8],
2299    posi: usize,
2300) -> impl Clone + Iterator<Item = (usize, usize)> {
2301    let pos = uniform_idx_as_vec2(map_size_lg, posi);
2302    let mrec_i = mrec[posi];
2303    NEIGHBOR_DELTA
2304        .iter()
2305        .enumerate()
2306        .filter(move |&(k, _)| (mrec_i >> k as isize) & 1 == 1)
2307        .map(move |(k, &(x, y))| {
2308            (
2309                k,
2310                vec2_as_uniform_idx(map_size_lg, Vec2::new(pos.x + x, pos.y + y)),
2311            )
2312        })
2313}
2314
2315/// Algorithm for computing multi-receiver flow.
2316///
2317/// * `map_size_lg`: Size of the underlying map.
2318/// * `h`: altitude
2319/// * `downhill`: single receiver
2320/// * `newh`: single receiver stack
2321/// * `wh`: buffer into which water height will be inserted.
2322/// * `nx`, `ny`: resolution in x and y directions.
2323/// * `dx`, `dy`: grid spacing in x- and y-directions
2324/// * `maxh`: maximum |height| among all nodes.
2325///
2326///
2327/// Updates the water height to a nearly planar surface, and returns a 3-tuple
2328/// containing:
2329/// * A bitmask representing which neighbors are downhill.
2330/// * Stack order for multiple receivers (from top to bottom).
2331/// * The weight for each receiver, for each node.
2332pub fn get_multi_rec<F: fmt::Debug + Float + Sync + Into<Compute>>(
2333    map_size_lg: MapSizeLg,
2334    h: impl Fn(usize) -> F + Sync,
2335    downhill: &[isize],
2336    newh: &[u32],
2337    wh: &mut [F],
2338    nx: usize,
2339    ny: usize,
2340    dx: Compute,
2341    dy: Compute,
2342    _maxh: F,
2343    threadpool: &rayon::ThreadPool,
2344) -> (Box<[u8]>, Box<[u32]>, Box<[Computex8]>) {
2345    let nn = nx * ny;
2346    let dxdy = Vec2::new(dx, dy);
2347
2348    /* // set bc
2349    let i1 = 0;
2350    let i2 = nx;
2351    let j1 = 0;
2352    let j2 = ny;
2353    let xcyclic = false;
2354    let ycyclic = false; */
2355    /*
2356      write (cbc,'(i4)') ibc
2357      i1=1
2358      i2=nx
2359      j1=1
2360      j2=ny
2361      if (cbc(4:4).eq.'1') i1=2
2362      if (cbc(2:2).eq.'1') i2=nx-1
2363      if (cbc(1:1).eq.'1') j1=2
2364      if (cbc(3:3).eq.'1') j2=ny-1
2365      xcyclic=.FALSE.
2366      ycyclic=.FALSE.
2367      if (cbc(4:4).ne.'1'.and.cbc(2:2).ne.'1') xcyclic=.TRUE.
2368      if (cbc(1:1).ne.'1'.and.cbc(3:3).ne.'1') ycyclic=.TRUE.
2369    */
2370    assert_eq!(nn, wh.len());
2371
2372    // fill the local minima with a nearly planar surface
2373    // See https://matthew-brett.github.io/teaching/floating_error.html;
2374    // our absolute error is bounded by β^(e-(p-1)), where e is the exponent of the
2375    // largest value we care about.  In this case, since we are summing up to nn
2376    // numbers, we are bounded from above by nn * |maxh|; however, we only need
2377    // to invoke this when we actually encounter a number, so we compute it
2378    // dynamically. for nn + |maxh|
2379    // TODO: Consider that it's probably not possible to have a downhill path the
2380    // size of the whole grid... either measure explicitly (maybe in get_lakes)
2381    // or work out a more precise upper bound (since using nn * 2 * (maxh +
2382    // epsilon) makes f32 not work very well).
2383    let deltah = F::epsilon() + F::epsilon();
2384    newh.iter().for_each(|&ijk| {
2385        let ijk = ijk as usize;
2386        let h_i = h(ijk);
2387        let ijr = downhill[ijk];
2388        wh[ijk] = if ijr >= 0 {
2389            let ijr = ijr as usize;
2390            let wh_j = wh[ijr];
2391            if wh_j >= h_i {
2392                let deltah = deltah * wh_j.abs();
2393                wh_j + deltah
2394            } else {
2395                h_i
2396            }
2397        } else {
2398            h_i
2399        };
2400    });
2401
2402    let mut wrec = Vec::<Computex8>::with_capacity(nn);
2403    let mut mrec = Vec::with_capacity(nn);
2404    let mut don_vis = Vec::with_capacity(nn);
2405
2406    // loop on all nodes
2407    (0..nn)
2408        .into_par_iter()
2409        .map(|ij| {
2410            // TODO: SIMDify?  Seems extremely amenable to that.
2411            let wh_ij = wh[ij];
2412            let mut mrec_ij = 0u8;
2413            let mut ndon_ij = 0u8;
2414            let neighbor_iter = |posi| {
2415                let pos = uniform_idx_as_vec2(map_size_lg, posi);
2416                NEIGHBOR_DELTA
2417                    .iter()
2418                    .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y))
2419                    .enumerate()
2420                    .filter(move |&(_, pos)| {
2421                        pos.x >= 0 && pos.y >= 0 && pos.x < nx as i32 && pos.y < ny as i32
2422                    })
2423                    .map(move |(k, pos)| (k, vec2_as_uniform_idx(map_size_lg, pos)))
2424            };
2425
2426            neighbor_iter(ij).for_each(|(k, ijk)| {
2427                let wh_ijk = wh[ijk];
2428                if wh_ij > wh_ijk {
2429                    // Set neighboring edge lower than this one as being downhill.
2430                    // NOTE: relying on at most 8 neighbors.
2431                    mrec_ij |= 1 << k;
2432                } else if wh_ijk > wh_ij {
2433                    // Avoiding ambiguous cases.
2434                    ndon_ij += 1;
2435                }
2436            });
2437            (mrec_ij, (ndon_ij, ndon_ij))
2438        })
2439        .unzip_into_vecs(&mut mrec, &mut don_vis);
2440
2441    let czero = <Compute as Zero>::zero();
2442    let (wrec, stack) = threadpool.join(
2443        || {
2444            (0..nn)
2445                .into_par_iter()
2446                .map(|ij| {
2447                    let mut sumweight = czero;
2448                    let mut wrec = [czero; 8];
2449                    let mut nrec = 0;
2450                    mrec_downhill(map_size_lg, &mrec, ij).for_each(|(k, ijk)| {
2451                        let lrec_ijk = ((uniform_idx_as_vec2(map_size_lg, ijk)
2452                            - uniform_idx_as_vec2(map_size_lg, ij))
2453                        .map(|e| e as Compute)
2454                            * dxdy)
2455                            .magnitude();
2456                        let wrec_ijk = (wh[ij] - wh[ijk]).into() / lrec_ijk;
2457                        // NOTE: To emulate single-direction flow, uncomment this line.
2458                        // let wrec_ijk = if ijk as isize == downhill[ij] { <Compute as One>::one()
2459                        // } else { <Compute as Zero>::zero() };
2460                        wrec[k] = wrec_ijk;
2461                        sumweight += wrec_ijk;
2462                        nrec += 1;
2463                    });
2464                    let slope = sumweight / (nrec as Compute).max(1.0);
2465                    let p_mfd_exp = 0.5 + 0.6 * slope;
2466                    sumweight = czero;
2467                    wrec.iter_mut().for_each(|wrec_k| {
2468                        let wrec_ijk = wrec_k.powf(p_mfd_exp);
2469                        sumweight += wrec_ijk;
2470                        *wrec_k = wrec_ijk;
2471                    });
2472                    if sumweight > czero {
2473                        wrec.iter_mut().for_each(|wrec_k| {
2474                            *wrec_k /= sumweight;
2475                        });
2476                    }
2477                    wrec
2478                })
2479                .collect_into_vec(&mut wrec);
2480            wrec
2481        },
2482        || {
2483            let mut stack = Vec::with_capacity(nn);
2484            let mut parse = Vec::with_capacity(nn);
2485
2486            // we go through the nodes
2487            (0..nn).for_each(|ij| {
2488                let (ndon_ij, _) = don_vis[ij];
2489                // when we find a "summit" (ie a node that has no donors)
2490                // we parse it (put it in a stack called parse)
2491                if ndon_ij == 0 {
2492                    parse.push(ij);
2493                }
2494                // we go through the parsing stack
2495                while let Some(ijn) = parse.pop() {
2496                    // we add the node to the stack
2497                    stack.push(ijn as u32);
2498                    mrec_downhill(map_size_lg, &mrec, ijn).for_each(|(_, ijr)| {
2499                        let (_, ref mut vis_ijr) = don_vis[ijr];
2500                        if *vis_ijr >= 1 {
2501                            *vis_ijr -= 1;
2502                            if *vis_ijr == 0 {
2503                                parse.push(ijr);
2504                            }
2505                        }
2506                    });
2507                }
2508            });
2509
2510            assert_eq!(stack.len(), nn);
2511            stack
2512        },
2513    );
2514
2515    (
2516        mrec.into_boxed_slice(),
2517        stack.into_boxed_slice(),
2518        wrec.into_boxed_slice(),
2519    )
2520}
2521
2522/// Perform erosion n times.
2523#[expect(clippy::too_many_arguments)]
2524pub fn do_erosion(
2525    map_size_lg: MapSizeLg,
2526    _max_uplift: f32,
2527    n_steps: usize,
2528    seed: &RandomField,
2529    rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
2530    oldh: impl Fn(usize) -> f32 + Sync,
2531    oldb: impl Fn(usize) -> f32 + Sync,
2532    is_ocean: impl Fn(usize) -> bool + Sync,
2533    uplift: impl Fn(usize) -> f64 + Sync,
2534    n: impl Fn(usize) -> f32 + Sync,
2535    theta: impl Fn(usize) -> f32 + Sync,
2536    kf: impl Fn(usize) -> f64 + Sync,
2537    kd: impl Fn(usize) -> f64 + Sync,
2538    g: impl Fn(usize) -> f32 + Sync,
2539    epsilon_0: impl Fn(usize) -> f32 + Sync,
2540    alpha: impl Fn(usize) -> f32 + Sync,
2541    // scaling factors
2542    height_scale: impl Fn(f32) -> Alt + Sync,
2543    k_d_scale: f64,
2544    k_da_scale: impl Fn(f64) -> f64,
2545    threadpool: &rayon::ThreadPool,
2546    report_progress: &dyn Fn(f64),
2547) -> (Box<[Alt]>, Box<[Alt]> /* , Box<[Alt]> */) {
2548    debug!("Initializing erosion arrays...");
2549    let oldh_ = (0..map_size_lg.chunks_len())
2550        .into_par_iter()
2551        .map(|posi| oldh(posi) as Alt)
2552        .collect::<Vec<_>>()
2553        .into_boxed_slice();
2554    // Topographic basement (The height of bedrock, not including sediment).
2555    let mut b = (0..map_size_lg.chunks_len())
2556        .into_par_iter()
2557        .map(|posi| oldb(posi) as Alt)
2558        .collect::<Vec<_>>()
2559        .into_boxed_slice();
2560    // Stream power law slope exponent--link between channel slope and erosion rate.
2561    let n = (0..map_size_lg.chunks_len())
2562        .into_par_iter()
2563        .map(&n)
2564        .collect::<Vec<_>>()
2565        .into_boxed_slice();
2566    // Stream power law concavity index (θ = m/n), turned into an exponent on
2567    // drainage (which is a proxy for discharge according to Hack's Law).
2568    let m = (0..map_size_lg.chunks_len())
2569        .into_par_iter()
2570        .map(|posi| theta(posi) * n[posi])
2571        .collect::<Vec<_>>()
2572        .into_boxed_slice();
2573    // Stream power law erodability constant for fluvial erosion (bedrock)
2574    let kf = (0..map_size_lg.chunks_len())
2575        .into_par_iter()
2576        .map(&kf)
2577        .collect::<Vec<_>>()
2578        .into_boxed_slice();
2579    // Stream power law erodability constant for hillslope diffusion (bedrock)
2580    let kd = (0..map_size_lg.chunks_len())
2581        .into_par_iter()
2582        .map(&kd)
2583        .collect::<Vec<_>>()
2584        .into_boxed_slice();
2585    // Deposition coefficient
2586    let g = (0..map_size_lg.chunks_len())
2587        .into_par_iter()
2588        .map(&g)
2589        .collect::<Vec<_>>()
2590        .into_boxed_slice();
2591    let epsilon_0 = (0..map_size_lg.chunks_len())
2592        .into_par_iter()
2593        .map(&epsilon_0)
2594        .collect::<Vec<_>>()
2595        .into_boxed_slice();
2596    let alpha = (0..map_size_lg.chunks_len())
2597        .into_par_iter()
2598        .map(&alpha)
2599        .collect::<Vec<_>>()
2600        .into_boxed_slice();
2601    let mut wh = vec![0.0; map_size_lg.chunks_len()].into_boxed_slice();
2602    // TODO: Don't do this, maybe?
2603    // (To elaborate, maybe we should have varying uplift or compute it some other
2604    // way).
2605    let uplift = (0..oldh_.len())
2606        .into_par_iter()
2607        .map(|posi| uplift(posi) as f32)
2608        .collect::<Vec<_>>()
2609        .into_boxed_slice();
2610    let sum_uplift = uplift
2611        .into_par_iter()
2612        .cloned()
2613        .map(|e| e as f64)
2614        .sum::<f64>();
2615    debug!("Sum uplifts: {:?}", sum_uplift);
2616
2617    let max_uplift = uplift
2618        .into_par_iter()
2619        .cloned()
2620        .max_by(|a, b| a.partial_cmp(b).unwrap())
2621        .unwrap();
2622    let max_g = g
2623        .into_par_iter()
2624        .cloned()
2625        .max_by(|a, b| a.partial_cmp(b).unwrap())
2626        .unwrap();
2627    debug!("Max uplift: {:?}", max_uplift);
2628    debug!("Max g: {:?}", max_g);
2629    // Height of terrain, including sediment.
2630    let mut h = oldh_;
2631    // Bedrock transport coefficients (diffusivity) in m^2 / year, for sediment.
2632    // For now, we set them all to be equal
2633    // on land, but in theory we probably want to at least differentiate between
2634    // soil, bedrock, and sediment.
2635    let kdsed = 1.5e-2 / 4.0;
2636    let kdsed = kdsed * k_d_scale;
2637    let n = |posi: usize| n[posi];
2638    let m = |posi: usize| m[posi];
2639    let kd = |posi: usize| kd[posi];
2640    let kf = |posi: usize| kf[posi];
2641    let g = |posi: usize| g[posi];
2642    let epsilon_0 = |posi: usize| epsilon_0[posi];
2643    let alpha = |posi: usize| alpha[posi];
2644
2645    (0..n_steps).for_each(|i| {
2646        debug!("Erosion iteration #{:?}", i);
2647
2648        // Print out the percentage complete. Do this at most 20 times.
2649        if i % std::cmp::max(n_steps / 20, 1) == 0 {
2650            let pct = (i as f64 / n_steps as f64) * 100.0;
2651            report_progress(pct);
2652            info!("{:.2}% complete", pct);
2653        }
2654
2655        erode(
2656            map_size_lg,
2657            &mut h,
2658            &mut b,
2659            &mut wh,
2660            max_uplift,
2661            max_g,
2662            kdsed,
2663            seed,
2664            rock_strength_nz,
2665            |posi| uplift[posi],
2666            n,
2667            m,
2668            kf,
2669            kd,
2670            g,
2671            epsilon_0,
2672            alpha,
2673            &is_ocean,
2674            &height_scale,
2675            &k_da_scale,
2676            threadpool,
2677        );
2678    });
2679    (h, b)
2680}