Fix sediment transport, add hack for sediment.

world/src/sim/diffusion.rs

@@ -1,4 +1,5 @@
 use rayon::prelude::*;
+use super::Alt;
 
 /// From https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
 ///
@@ -36,7 +37,7 @@ use rayon::prelude::*;
   implicit none
 */
 pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
-                 h: &mut [/*f64*/f32], b: &mut [/*f64*/f32],
+                 h: &mut [Alt], b: &mut [Alt],
                  kd: impl Fn(usize) -> f64, kdsed: f64,
                  ) {
     let aij = |i: usize, j: usize| j * nx + i;
@@ -370,7 +371,7 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
         for i in 0..nx {
             ij = aij(i, j);
             // FIXME: Track total erosion and erosion rate.
-            h[ij] = zintp[ij] as f32;
+            h[ij] = zintp[ij] as Alt;
         }
     }
 

world/src/sim/erosion.rs

@@ -13,6 +13,9 @@ use std::{
 };
 use vek::*;
 
+pub type Alt = f32;
+pub type Compute = f64;
+
 /// Compute the water flux at all chunks, given a list of chunk indices sorted by increasing
 /// height.
 pub fn get_drainage(newh: &[u32], downhill: &[isize], _boundary_len: usize) -> Box<[f32]> {
@@ -288,7 +291,12 @@ pub fn get_rivers(
 
         // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
         // pushed towards the point identified by pass_idx).
-        let pass_idx = (-indirection[lake_idx]) as usize;
+        let pass_idx = if downhill[lake_idx] < 0 {
+            // Flows into nothing, so this lake is its own pass.
+            lake_idx
+        } else {
+            (-indirection[lake_idx]) as usize
+        };
 
         // Add our spline derivative to the downhill river (weighted by the chunk's drainage).
         // NOTE: Don't add the spline derivative to the lake side of the pass for our own lake,
@@ -486,7 +494,7 @@ pub fn get_rivers(
 /// Precompute the maximum slope at all points.
 ///
 /// TODO: See if allocating in advance is worthwhile.
-fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
+fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
     const MIN_MAX_ANGLE: f64 = 15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
     const MAX_MAX_ANGLE: f64 = 60.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
     const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE;
@@ -606,9 +614,9 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
 ///
 ///
 fn erode(
-    h: &mut [f32],
+    h: &mut [Alt],
-    b: &mut [f32],
+    b: &mut [Alt],
-    wh: &mut [f32],
+    wh: &mut [Alt],
     is_done: &mut BitBox,
     done_val: bool,
     erosion_base: f32,
@@ -623,6 +631,7 @@ fn erode(
     g: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
 ) {
+    let compute_stats = true;
     log::debug!("Done draining...");
     let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
     let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
@@ -671,15 +680,15 @@ fn erode(
         // (2.244*(5.010e-4)/512)/(2.244*(5.010e-4)/2500) = 4.88...
         // 2.444 * 5
     // Stream power erosion constant (sediment), in m^(1-2m) / year (times dt).
-    let k_fs_mult = 2.0;/*1.5*/;
+    let k_fs_mult = 2.0;//2.0;/*1.5*/;
     // let k_fs = k_fb * 1.0/*1.5*//*2.0*//*2.0*//*4.0*/;
     // u = k * h_max / 2.244
     // let uplift_scale = erosion_base as f64 + (k_fb * mmaxh / 2.244 / 5.010e-4 as f64 * mmaxh as f64) * dt;
     let ((dh, indirection, newh, maxh, area), (mut max_slopes, ht)) = rayon::join(
         || {
-            let mut dh = downhill(h, |posi| is_ocean(posi) && h[posi] <= 0.0);
+            let mut dh = downhill(|posi| h[posi], |posi| is_ocean(posi) && h[posi] <= 0.0);
             log::debug!("Computed downhill...");
-            let (boundary_len, indirection, newh, maxh) = get_lakes(&h, &mut dh);
+            let (boundary_len, indirection, newh, maxh) = get_lakes(|posi| h[posi], &mut dh);
             log::debug!("Got lakes...");
             let area = get_drainage(&newh, &dh, boundary_len);
             log::debug!("Got flux...");
@@ -695,6 +704,7 @@ fn erode(
                 || {
                     // Store the elevation at t
                     h.to_vec().into_boxed_slice()
+                    // h.into_par_iter().map(|e| e as f64).collect::<Vec<_>>().into_boxed_slice()
                 },
             )
         },
@@ -707,13 +717,13 @@ fn erode(
     // TODO: Make more principled.
     let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();//1.0;
 
-    let mut lake_water_volume = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut lake_water_volume = vec![/*-1i32*/0.0 as Compute; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
-    let mut elev = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut elev = vec![/*-1i32*/0.0 as Compute; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
-    let mut hp = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut hp = vec![/*-1i32*/0.0 as Compute; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
-    let mut deltah = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut deltah = vec![/*-1i32*/0.0 as Compute; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
 
     // calculate the elevation / SPL, including sediment flux
-    let tol = 1.0e-4f64 * (maxh as f64 + 1.0);
+    let tol = 1.0e-4 as Compute * (maxh as Compute + 1.0);
     let mut err = 2.0 * tol;
 
     // Some variables for tracking statistics, currently only for debugging purposes.
@@ -721,11 +731,14 @@ fn erode(
     let mut nland = 0usize;
     let mut sums = 0.0;
     let mut sumh = 0.0;
+    let mut sumsed = 0.0;
+    let mut sumsed_land = 0.0;
     let mut ntherm = 0usize;
+    let avgz = |x, y: usize| if y == 0 { f64::INFINITY } else { x / y as f64 };
 
     // Gauss-Seidel iteration
 
-    let mut lake_sediment = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut lake_sediment = vec![/*-1i32*/0.0 as Compute; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
     let mut lake_sill = vec![/*-1i32*/-1isize; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
 
     let mut n_gs_stream_power_law = 0;
@@ -737,6 +750,8 @@ fn erode(
         nland = 0usize;
         sums = 0.0;
         sumh = 0.0;
+        sumsed = 0.0;
+        sumsed_land = 0.0;
         ntherm = 0usize;
 
         // Keep track of how many iterations we've gone to test for convergence.
@@ -746,13 +761,14 @@ fn erode(
             || {
                 // guess/update the elevation at t+Δt (k)
                 hp.par_iter_mut().zip(h.par_iter()).for_each(|(mut hp, h)| {
-                    *hp = *h as f64;
+                    *hp = *h as Compute;
                 });
             },
             || {
                 // calculate erosion/deposition at each node
                 deltah.par_iter_mut().enumerate().for_each(|(posi, mut deltah)| {
-                    *deltah = (ht[posi] - h[posi]) as f64;
+                    let uplift_i = uplift(posi) as Alt;
+                    *deltah = (ht[posi] + uplift_i - h[posi]) as Compute;
                 });
             },
         );
@@ -764,8 +780,9 @@ fn erode(
             if posj < 0 {
                 lake_sediment[posi] = deltah[posi];
             } else {
+                let uplift_i = uplift(posi) as Alt;
                 let posj = posj as usize;
-                deltah[posi] -= (ht[posi] as f64 - hp[posi]);
+                deltah[posi] -= ((ht[posi] + uplift_i) as Compute - hp[posi]);
                 let lposi = lake_sill[posi];
                 if lposi == posi as isize {
                     if deltah[posi] <= 0.0 {
@@ -774,7 +791,7 @@ fn erode(
                         lake_sediment[posi] = deltah[posi];
                     }
                 }
-                deltah[posi] += ht[posi] as f64 - hp[posi];
+                deltah[posi] += (ht[posi] + uplift_i) as Compute - hp[posi];
                 deltah[posj] += deltah[posi];
             }
         }
@@ -799,9 +816,11 @@ fn erode(
 
         elev.par_iter_mut().enumerate().for_each(|(posi, mut elev)| {
             if dh[posi] < 0 {
-                *elev = ht[posi] as f64;
+                *elev = ht[posi] as Compute;
             } else {
-                *elev = ht[posi] as f64 + (deltah[posi] - (ht[posi] as f64 - hp[posi])) * g(posi) as f64 / area[posi] as f64;
+                let uplift_i = uplift(posi) as Alt;
+                assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.0);
+                *elev = (ht[posi] + uplift_i) as Compute + (deltah[posi] - ((ht[posi] + uplift_i) as Compute - hp[posi])) * g(posi) as Compute / area[posi] as Compute;
             }
         });
 
@@ -809,6 +828,10 @@ fn erode(
         let mut sum_err = 0.0f64;
         for &posi in &*newh {
             let posi = posi as usize;
+            let old_h_i = /*h*/elev[posi] as f64;
+            let old_b_i = b[posi];
+            let sed = (ht[posi] - old_b_i) as f64;
+
             let posj = dh[posi];
             if posj < 0 {
                 if posj == -1 {
@@ -833,20 +856,18 @@ fn erode(
                 // want to multiply neighbor_distance by height_scale and area[posi] by
                 // height_scale^2, to make sure we were working in the correct units for dz
                 // (which has height_scale height unit = 1.0 meters).
-                let old_h_i = /*h*/elev[posi] as f64;
+                /* let uplift_i = uplift(posi) as f64;
-                let old_b_i = b[posi] as f64;
+                assert!(uplift_i.is_normal() && uplift_i == 0.0 || uplift_i.is_positive()); */
-                let uplift_i = uplift(posi) as f64;
-                assert!(uplift_i.is_normal() && uplift_i == 0.0 || uplift_i.is_positive());
                 // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
                 // NOTE: posj has already been computed since it's downhill from us.
                 // Therefore, we can rely on wh being set to the water height for that node.
                 let h_j = h[posj] as f64;
                 let wh_j = wh[posj] as f64;
-                let mut new_h_i = old_h_i + uplift_i;
+                let mut new_h_i = old_h_i/* + uplift_i*/;
                 // Only perform erosion if we are above the water level of the previous node.
                 if new_h_i > wh_j {
                     // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
-                    let k = if (old_h_i - old_b_i as f64) > sediment_thickness {
+                    let k = if sed > sediment_thickness {
                         // Sediment
                         // k_fs
                         k_fs_mult * kf(posi) as f64
@@ -861,18 +882,34 @@ fn erode(
                     new_h_i = (new_h_i + flux * h_j) / (1.0 + flux);
                     lake_sill[posi] = posi as isize;
                     lake_water_volume[posi] = 0.0;
+
+                    /* // Thermal erosion (landslide)
+                    let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                    let mag_slope = dz/*.abs()*/ / neighbor_distance;
+                    let max_slope = max_slopes[posi] as f64;
+                    if mag_slope > max_slope {
+                        let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
+                        new_h_i = (ht[posi] as f64 + l * (mag_slope - max_slope));
+                        if compute_stats && new_h_i > wh_j {
+                            ntherm += 1;
+                        }
+                    } */
+
                     // If we dipped below the receiver's water level, set our height to the receiver's
                     // water level.
                     if new_h_i <= wh_j {
                         new_h_i = wh_j;
                     } else {
+                        if compute_stats {
                             let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                             let mag_slope = dz/*.abs()*/ / neighbor_distance;
 
                             nland += 1;
+                            sumsed_land += sed;
                             sumh += new_h_i;
                             sums += mag_slope;
                         }
+                    }
                 } else {
                     let lposj = lake_sill[posj];
                     lake_sill[posi] = lposj;
@@ -883,7 +920,7 @@ fn erode(
                 }
                 // Set max_slope to this node's water height (max of receiver's water height and
                 // this node's height).
-                wh[posi] = wh_j.max(new_h_i) as f32;
+                wh[posi] = wh_j.max(new_h_i) as Alt;
                 // Prevent erosion from dropping us below our receiver, unless we were already below it.
                 // new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
                 // Find out if this is a lake bottom.
@@ -940,19 +977,22 @@ fn erode(
                         sums += mag_slope;
                         // Just use the computed rate.
                     } */
-                    h[posi] = new_h_i as f32;
+                    h[posi] = new_h_i as Alt;
                     // Make sure to update the basement as well!
                     // b[posi] = (old_b_i + uplift_i).min(new_h_i) as f32;
                 }
             }
             // *is_done.at(posi) = done_val;
+            if compute_stats {
                 maxh = h[posi].max(maxh);
-
+                sumsed += sed;
-            // Add sum of squares of errors.
-            sum_err += (h[posi] as f64 - hp[posi]).powi(2);
             }
 
-        err = (sum_err / newh.len() as f64).sqrt();
+            // Add sum of squares of errors.
+            sum_err += (h[posi] as Compute - hp[posi]).powi(2);
+        }
+
+        err = (sum_err / newh.len() as Compute).sqrt();
         if max_g == 0.0 {
             err = 0.0;
         }
@@ -966,7 +1006,7 @@ fn erode(
     // update the basement
     b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| {
         let old_b_i = *b;
-        let uplift_i = uplift(posi);
+        let uplift_i = uplift(posi) as Alt;
 
         *b = (old_b_i + uplift_i).min(*h);
     });
@@ -982,7 +1022,7 @@ fn erode(
                 *h +=
                     (0.0.max(lake_sediment[lposi].min(lake_water_volume[lposi])) /
                     lake_water_volume[lposi] *
-                    (wh[posi] - *h) as f64) as f32;
+                    (wh[posi] - *h) as Compute) as Alt;
             }
         }
     });
@@ -995,18 +1035,14 @@ fn erode(
     // enddo
 
     log::info!(
-        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?}) (num land: {:?}) (num thermal: {:?})",
+        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n        \
+        (old avg sediment thickness [all/land]: {:?} / {:?})\n        \
+        (num land: {:?}) (num thermal: {:?})",
         maxh,
-        if nland == 0 {
+        avgz(sumh, nland),
-            f64::INFINITY
+        avgz(sums, nland),
-        } else {
+        avgz(sumsed, newh.len()),
-            sumh / nland as f64
+        avgz(sumsed_land, nland),
-        },
-        if nland == 0 {
-            f64::INFINITY
-        } else {
-            sums / nland as f64
-        },
         nland,
         ntherm,
     );
@@ -1015,25 +1051,30 @@ fn erode(
     maxh = 0.0;
     sumh = 0.0;
     sums = 0.0;
+    sumsed = 0.0;
+    sumsed_land = 0.0;
     nland = 0usize;
     ntherm = 0usize;
     for &posi in &*newh {
         let posi = posi as usize;
-        let posj = dh[posi];
+        let old_h_i = h/*b*/[posi] as f64;
-        let old_h_i = /*h*/b[posi] as f64;
         let old_b_i = b[posi] as f64;
-        let max_slope = max_slopes[posi] as f64;
+        let sed = (old_h_i - old_b_i) as f64;
+
+        let max_slope = max_slopes[posi];
         // Remember k_d for this chunk in max_slopes.
         let kd_factor =
             // 1.0;
-            (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/1.0))/*.min(kdsed)*/;
+            (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/;
-        max_slopes[posi] = if (old_h_i - old_b_i as f64) > sediment_thickness && kdsed > 0.0 {
+        max_slopes[posi] = if sed > sediment_thickness && kdsed > 0.0 {
             // Sediment
             kdsed/* * kd_factor*/
         } else {
             // Bedrock
             kd(posi) as f64 / kd_factor
         };
+
+        let posj = dh[posi];
         if posj < 0 {
             if posj == -1 {
                 panic!("Disconnected lake!");
@@ -1047,7 +1088,7 @@ fn erode(
             let wh_j = wh[posj] as f64;
             // If you're on the lake bottom and not right next to your neighbor, don't compute a
             // slope.
-            let mut new_h_i = old_h_i;//old_b_i;
+            let mut new_h_i = old_h_i;/*old_b_i;*/
             if
             /* !is_lake_bottom */ /* !fake_neighbor */
             wh_j < old_h_i {
@@ -1090,10 +1131,12 @@ fn erode(
                     // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
                     // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
                     let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
-                    new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
+                    new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*/h_k + dh + l * (mag_slope - max_slope));
+                    // new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
                     if new_h_i <= wh_j {
                         new_h_i = wh_j;
                     } else {
+                        if compute_stats {
                             let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                             let slope = dz/*.abs()*/ / neighbor_distance;
                             sums += slope;
@@ -1108,11 +1151,15 @@ fn erode(
                             // max_slopes[posi] *= kd_factor;
                             nland += 1;
                             sumh += new_h_i;
+                            sumsed_land += sed;
+                        }
                         // let slope = dz.signum() * max_slope;
                         // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
                         // sums += max_slope;
                     }
+                    if compute_stats {
                         ntherm += 1;
+                    }
                 } else {
                     /*if kd_factor < 1.0 {
                         max_slopes[posi] *= kd_factor;
@@ -1121,37 +1168,39 @@ fn erode(
                         max_slopes[posi] *= kd_factor;
                     } */
                     // max_slopes[posi] *= kd_factor;
+                    if compute_stats {
                         sums += mag_slope;
                         // Just use the computed rate.
                         nland += 1;
                         sumh += new_h_i;
+                        sumsed_land += sed;
                     }
-                // h/*b*/[posi] = old_h_i.min(new_h_i) as f32;
+                }
+                h/*b*/[posi] = old_h_i.min(new_h_i) as Alt;
                 // Make sure to update the basement as well!
                 // b[posi] = old_b_i.min(new_h_i) as f32;
-                b[posi] = old_b_i.min(old_b_i + (new_h_i - old_h_i)) as f32;
+                b[posi] = old_b_i.min(old_b_i + (new_h_i - old_h_i)) as Alt;
                 // sumh += new_h_i;
             }
             // Set wh to this node's water height (max of receiver's water height and
             // this node's height).
-            wh[posi] = wh_j.max(new_h_i) as f32;
+            wh[posi] = wh_j.max(new_h_i) as Alt;
         }
         // *is_done.at(posi) = done_val;
+        if compute_stats {
+            sumsed += sed;
             maxh = h[posi].max(maxh);
         }
+    }
     log::debug!(
-        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?}) (num land: {:?}) (num thermal: {:?})",
+        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n        \
+        (avg sediment thickness [all/land]: {:?} / {:?})\n        \
+        (num land: {:?}) (num thermal: {:?})",
         maxh,
-        if nland == 0 {
+        avgz(sumh, nland),
-            f64::INFINITY
+        avgz(sums, nland),
-        } else {
+        avgz(sumsed, newh.len()),
-            sumh / nland as f64
+        avgz(sumsed_land, nland),
-        },
-        if nland == 0 {
-            f64::INFINITY
-        } else {
-            sums / nland as f64
-        },
         nland,
         ntherm,
     );
@@ -1257,7 +1306,7 @@ pub fn fill_sinks(
 /// - A list of chunks on the boundary (non-lake egress points).
 /// - The second indirection vector (associating chunk indices with their lake's adjacency list).
 /// - The adjacency list (stored in a single vector), indexed by the second indirection vector.
-pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[u32]>, f32) {
+pub fn get_lakes<F: Float>(h: impl Fn(usize) -> F, downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[u32]>, F) {
     // Associates each lake index with its root node (the deepest one in the lake), and a list of
     // adjacent lakes.  The list of adjacent lakes includes the lake index of the adjacent lake,
     // and a node index in the adjacent lake which has a neighbor in this lake.  The particular
@@ -1330,7 +1379,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
     log::debug!("Old lake roots: {:?}", lake_roots.len());
 
     let newh = newh.into_boxed_slice();
-    let mut maxh = -f32::INFINITY;
+    let mut maxh = -F::infinity();
     // Now, we know that the sum of all the indirection nodes will be the same as the number of
     // nodes.  We can allocate a *single* vector with 8 * nodes entries, to be used as storage
     // space, and augment our indirection vector with the starting index, resulting in a vector of
@@ -1340,7 +1389,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
         let chunk_idx = chunk_idx_ as usize;
         let lake_idx_ = indirection_[chunk_idx];
         let lake_idx = lake_idx_ as usize;
-        let height = h[chunk_idx_ as usize];
+        let height = h(chunk_idx_ as usize);
         // While we're here, compute the max elevation difference from zero among all nodes.
         maxh = maxh.max(height.abs());
         // For every neighbor, check to see whether it is already set; if the neighbor is set,
@@ -1350,7 +1399,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
         // of the heights we get out from this process is greater than the maximum of this
         // chunk and its neighbor chunk, we switch to this new edge.
         for neighbor_idx in neighbors(chunk_idx) {
-            let neighbor_height = h[neighbor_idx];
+            let neighbor_height = h(neighbor_idx);
             let neighbor_lake_idx_ = indirection_[neighbor_idx];
             let neighbor_lake_idx = neighbor_lake_idx_ as usize;
             if neighbor_lake_idx_ < lake_idx_ {
@@ -1384,8 +1433,8 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
                             // Should never run into -1 while looping here, since (i, j)
                             // and (j, i) should be added together.
                             if indirection_[neighbor_pass.1 as usize] == lake_idx_ {
-                                let pass_height = h[neighbor_pass.1 as usize];
+                                let pass_height = h(neighbor_pass.1 as usize);
-                                let neighbor_pass_height = h[pass.1 as usize];
+                                let neighbor_pass_height = h(pass.1 as usize);
                                 if height.max(neighbor_height)
                                     < pass_height.max(neighbor_pass_height)
                                 {
@@ -1477,7 +1526,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
         };
         if downhill[lake_idx] == -2 {
             // Find the pass height
-            let pass = h[from].max(h[to]);
+            let pass = h(from).max(h(to));
             candidates.push(Reverse((
                 NotNan::new(pass).unwrap(),
                 (edge.0 as u32, edge.1),
@@ -1546,7 +1595,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
                 continue;
             }
             // Find the pass height
-            let pass = h[edge.0 as usize].max(h[edge.1 as usize]);
+            let pass = h(edge.0 as usize).max(h(edge.1 as usize));
             // Put the reverse edge in candidates, sorted by height, then chunk idx, and finally
             // neighbor idx.
             candidates.push(Reverse((
@@ -1655,16 +1704,16 @@ pub fn do_erosion(
     kf: impl Fn(usize) -> f32 + Sync,
     kd: impl Fn(usize) -> f32 + Sync,
     g: impl Fn(usize) -> f32 + Sync,
-) -> (Box<[f32]>, Box<[f32]>) {
+) -> (Box<[Alt]>, Box<[Alt]>) {
     let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
-        .map(|posi| oldh(posi))
+        .map(|posi| oldh(posi) as Alt)
         .collect::<Vec<_>>()
         .into_boxed_slice();
     // Topographic basement (The height of bedrock, not including sediment).
     let mut b = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
-        .map(|posi| oldb(posi))
+        .map(|posi| oldb(posi) as Alt)
         .collect::<Vec<_>>()
         .into_boxed_slice();
     // Stream power law erodability constant for fluvial erosion (bedrock)

world/src/sim/mod.rs

@@ -7,7 +7,7 @@ mod util;
 // Reexports
 pub use self::diffusion::diffusion;
 pub use self::erosion::{
-    do_erosion, fill_sinks, get_drainage, get_lakes, get_rivers, RiverData, RiverKind,
+    Alt, do_erosion, fill_sinks, get_drainage, get_lakes, get_rivers, RiverData, RiverKind,
 };
 pub use self::location::Location;
 pub use self::settlement::Settlement;
@@ -62,8 +62,8 @@ struct GenCdf {
     humid_base: InverseCdf,
     temp_base: InverseCdf,
     chaos: InverseCdf,
-    alt: Box<[f32]>,
+    alt: Box<[Alt]>,
-    basement: Box<[f32]>,
+    basement: Box<[Alt]>,
     water_alt: Box<[f32]>,
     dh: Box<[isize]>,
     /// NOTE: Until we hit 4096 × 4096, this should suffice since integers with an absolute value
@@ -231,6 +231,7 @@ impl WorldSim {
                 .set_displacement(/*0.5*/1.0)
                 .set_range_function(RangeFunction::Euclidean)
                 // .enable_range(true),
+            // g_nz: RidgedMulti::new()
         };
 
         let river_seed = &gen_ctx.river_seed;
@@ -247,12 +248,12 @@ impl WorldSim {
             .set_scale(1.0 / rock_strength_scale); */
 
         let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
-        let max_erosion_per_delta_t = 128.0/*128.0*//*32.0*/ * height_scale;
+        let max_erosion_per_delta_t = /*8.0*/64.0/*128.0*//*32.0*/ * height_scale;
         let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
         let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
         let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
-        let n_steps = /*200*/50; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
+        let n_steps = 50; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
-        let n_small_steps = 8;//50;//8;//8;//8;//8;//8; // 8
+        let n_small_steps = 25;//50;//8;//8;//8;//8;//8; // 8
 
 
         // Logistic regression.  Make sure x ∈ (0, 1).
@@ -548,8 +549,9 @@ impl WorldSim {
                             // uheight_1;
                             // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5));
                             // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin()
-                            uheight /* * 0.8*/ * udist_1.powf(2.0) +
+                            /* (uheight /* * 0.8*//* * udist_1.powf(2.0)*/ +
-                            /*exp_inverse_cdf*/(oheight.max(0.0).min(max_epsilon).abs())/*.powf(2.0)*/ * 0.2 * udist_2.powf(2.0);
+                            /*exp_inverse_cdf*/(oheight./*max(0.0).min(max_epsilon).abs()*/)/*.powf(2.0)*/ * 0.2/* * udist_2.powf(2.0)*/; */
+                            (uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
                             // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
                         // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
 
@@ -645,8 +647,8 @@ impl WorldSim {
                     .add(0.5);
                 // kf = 1.5e-4: high-high (plateau [fan sediment])
                 // kf = 1e-4: high (plateau)
-                // kf = 2e-5: normal (dyke [unexposed])
+                // kf = 2e-5: normal (dike [unexposed])
-                // kf = 1e-6: normal-low (dyke [exposed])
+                // kf = 1e-6: normal-low (dike [exposed])
                 // kf = 2e-6: low (mountain)
                 ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
             }
@@ -678,7 +680,7 @@ impl WorldSim {
                     .max(-1.0)
                     .mul(0.5)
                     .add(0.5);
-                // kd = 1e-1: high (mountain, dyke)
+                // kd = 1e-1: high (mountain, dike)
                 // kd = 1.5e-2: normal-high (plateau [fan sediment])
                 // kd = 1e-2: normal (plateau)
                 1.0e-2
@@ -690,7 +692,41 @@ impl WorldSim {
                 if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 {
                     return 0.0;
                 }
-                1.0
+                let wposf = (uniform_idx_as_vec2(posi)
+                    * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+                .map(|e| e as f64);
+                /* let turb_wposf =
+                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                let turb = Vec2::new(
+                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                // let turb = Vec2::zero();
+                let turb_wposf = wposf + turb; */
+                let turb_wposf = wposf;
+                let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
+                    .min(1.0)
+                    .max(-1.0)
+                    .mul(0.5)
+                    .add(0.5); */
+                    chaos[posi].1;
+
+                assert!(uchaos <= 1.24);
+
+                // G = d* v_s / p_0, where
+                //  v_s is the settling velocity of sediment grains
+                //  p_0 is the mean precipitation rate
+                //  d* is the sediment concentration ratio (between concentration near riverbed
+                //  interface, and average concentration over the water column).
+                //  d* varies with Rouse number which defines relative contribution of bed, suspended,
+                //  and washed loads.
+                //
+                // G is typically on the order of 1 or greater.  However, we are only guaranteed to
+                // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.24].
+                ((1.24 - uchaos) / 1.24).powf(0.75) * 1.24
+                // 1.0
+                // 1.0
+                // 1.5
             };
         let uplift_fn =
             |posi| {
@@ -729,7 +765,7 @@ impl WorldSim {
                             .mul(0.045)*/)
                 };
                 let height =
-                    ((/*old_height_uniform*/uplift_uniform[posi].0/*1*/ - alt_old_min_uniform) as f64
+                    ((/*old_height_uniform*/uplift_uniform[posi]./*0*/1 - alt_old_min_uniform) as f64
                     / (alt_old_max_uniform - alt_old_min_uniform) as f64)
                     /*((old_height(posi) - alt_old_min) as f64
                     / (alt_old_max - alt_old_min) as f64)*/
@@ -787,20 +823,26 @@ impl WorldSim {
                 // tan(pi/6)*32 ~ 18
                 // tan(54/360*2*pi)*32
                 // let height = 1.0f64;
-                /* let turb_wposf =
+                let turb_wposf =
                     wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
                 let turb = Vec2::new(
                     gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
                     gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                ) * 64.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
                 let turb_wposf = wposf + turb;
-                let height = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
                     /* .min(0.5)
                     .max(-0.5)*/
                     .min(1.0)
                     .max(-1.0)
                     .mul(0.5)
-                    .add(0.5); */
+                    .add(0.5);
+                // u = 1e-3: normal-high (dike, mountain)
+                // u = 0.5: normal (mid example in Yuan, average mountain uplift)
+                // u = 0.2: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05).
+                // u = 0: low (plateau [fan, altitude = 0.0])
+                // let height = uheight;
+                // let height = 1.0f64;
 
                 // let height = 1.0 / 7.0f64;
                 // let height = 0.0 / 31.0f64;
@@ -866,7 +908,7 @@ impl WorldSim {
                 };
 
                 /* 0.0 */
-                old_height_uniform(posi)/*.powf(2.0)*/ - 0.5
+                (old_height_uniform(posi)/*.powf(2.0)*/ - 0.5)/* * CONFIG.mountain_scale as f32*/
                 // uplift_fn(posi) * (CONFIG.mountain_scale / max_erosion_per_delta_t as f32)
                 // 0.0
                 /* // 0.0
@@ -887,8 +929,8 @@ impl WorldSim {
             n_steps,
             &river_seed,
             &rock_strength_nz,
-            |posi| alt_func(posi),
+            |posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.0 },
-            |posi| alt_func(posi),// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
+            |posi| alt_func(posi) - if is_ocean_fn(posi) { 0.0 } else { 1400.0 },// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
             is_ocean_fn,
             uplift_fn,
             |posi| kf_func(posi),
@@ -902,8 +944,8 @@ impl WorldSim {
             n_small_steps,
             &river_seed,
             &rock_strength_nz,
-            |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi],
+            |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32,
-            |posi| basement[posi],
+            |posi| basement[posi] as f32,
             is_ocean_fn,
             |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
             |posi| kf_func(posi),
@@ -911,10 +953,10 @@ impl WorldSim {
             |posi| g_func(posi),
         );
 
-        let is_ocean = get_oceans(|posi| alt[posi]);
+        let is_ocean = get_oceans(|posi| alt[posi] as f32);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
-        let mut dh = downhill(&alt, /*old_height*/ is_ocean_fn);
+        let mut dh = downhill(|posi| alt[posi] as f32/*&alt*/, /*old_height*/ is_ocean_fn);
-        let (boundary_len, indirection, water_alt_pos, _) = get_lakes(&/*water_alt*/alt, &mut dh);
+        let (boundary_len, indirection, water_alt_pos, _) = get_lakes(/*&/*water_alt*/alt*/|posi| alt[posi] as f32, &mut dh);
         let flux_old = get_drainage(&water_alt_pos, &dh, boundary_len);
 
         let water_height_initial = |chunk_idx| {
@@ -958,7 +1000,7 @@ impl WorldSim {
             };
             // Use the maximum of the pass height and chunk height as the parameter to fill_sinks.
             let chunk_alt = alt[chunk_idx];
-            chunk_alt.max(chunk_water_alt)
+            chunk_alt.max(chunk_water_alt) as f32
         };
 
         let water_alt = fill_sinks(water_height_initial, is_ocean_fn);
@@ -1045,7 +1087,7 @@ impl WorldSim {
                             } else {
                                 // A version of alt that is uniform over *non-water* (or land-adjacent water)
                                 // chunks.
-                                Some(alt[posi])
+                                Some(alt[posi] as f32)
                             }
                         })
                     },
@@ -1579,8 +1621,8 @@ impl SimChunk {
 
         let _map_edge_factor = map_edge_factor(posi);
         let (_, chaos) = gen_cdf.chaos[posi];
-        let alt_pre = gen_cdf.alt[posi];
+        let alt_pre = gen_cdf.alt[posi] as f32;
-        let basement_pre = gen_cdf.basement[posi];
+        let basement_pre = gen_cdf.basement[posi] as f32;
         let water_alt_pre = gen_cdf.water_alt[posi];
         let downhill_pre = gen_cdf.dh[posi];
         let flux = gen_cdf.flux[posi];

world/src/sim/util.rs

@@ -247,12 +247,13 @@ pub fn uphill<'a>(dh: &'a [isize], posi: usize) -> impl Clone + Iterator<Item =
 /// Compute the neighbor "most downhill" from all chunks.
 ///
 /// TODO: See if allocating in advance is worthwhile.
-pub fn downhill(h: &[f32], is_ocean: impl Fn(usize) -> bool + Sync) -> Box<[isize]> {
+pub fn downhill<F: Float>(h: impl Fn(usize) -> F + Sync, is_ocean: impl Fn(usize) -> bool + Sync) -> Box<[isize]> {
     // Constructs not only the list of downhill nodes, but also computes an ordering (visiting
     // nodes in order from roots to leaves).
-    h.par_iter()
+    (0..WORLD_SIZE.x * WORLD_SIZE.y).into_par_iter()
-        .enumerate()
+        // .enumerate()
-        .map(|(posi, &nh)| {
+        .map(|(posi/*, &nh*/)| {
+            let nh = h(posi);
             let _pos = uniform_idx_as_vec2(posi);
             if is_ocean(posi) {
                 -2
@@ -260,7 +261,7 @@ pub fn downhill(h: &[f32], is_ocean: impl Fn(usize) -> bool + Sync) -> Box<[isiz
                 let mut best = -1;
                 let mut besth = nh;
                 for nposi in neighbors(posi) {
-                    let nbh = h[nposi];
+                    let nbh = h(nposi);
                     if nbh < besth {
                         besth = nbh;
                         best = nposi as isize;
@@ -285,6 +286,8 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
         let posi = vec2_as_uniform_idx(pos);
         if oldh(posi) <= 0.0 {
             stack.push(posi);
+        } else {
+            panic!("Hopefully no border tiles are above sea level.");
         }
     };
     for x in 0..WORLD_SIZE.x as i32 {