From c92ff34e15ca45adc612fd3fb4d91cbd656bf301 Mon Sep 17 00:00:00 2001 From: Joshua Yanovski Date: Tue, 3 Dec 2019 19:14:29 +0100 Subject: [PATCH] Fix sediment transport, add hack for sediment. --- world/src/sim/diffusion.rs | 5 +- world/src/sim/erosion.rs | 251 ++++++++++++++++++++++--------------- world/src/sim/mod.rs | 100 ++++++++++----- world/src/sim/util.rs | 13 +- 4 files changed, 232 insertions(+), 137 deletions(-) diff --git a/world/src/sim/diffusion.rs b/world/src/sim/diffusion.rs index 30a9503f0f..3393f423f2 100644 --- a/world/src/sim/diffusion.rs +++ b/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; } } diff --git a/world/src/sim/erosion.rs b/world/src/sim/erosion.rs index db9997d45c..91bddc4906 100644 --- a/world/src/sim/erosion.rs +++ b/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> + Sync)) -> Box<[f64]> { +fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn> + 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> + Sync /// /// fn erode( - h: &mut [f32], - b: &mut [f32], - wh: &mut [f32], + h: &mut [Alt], + b: &mut [Alt], + 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::>().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 elev = vec![/*-1i32*/0.0f64; 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 deltah = 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.0 as Compute; 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.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 old_b_i = b[posi] as f64; - let uplift_i = uplift(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,17 +882,33 @@ 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 { - 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; + 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; - sumh += new_h_i; - sums += mag_slope; + nland += 1; + sumsed_land += sed; + sumh += new_h_i; + sums += mag_slope; + } } } else { let lposj = lake_sill[posj]; @@ -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; - maxh = h[posi].max(maxh); + 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); + sum_err += (h[posi] as Compute - hp[posi]).powi(2); } - err = (sum_err / newh.len() as f64).sqrt(); + 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 { - f64::INFINITY - } else { - sumh / nland as f64 - }, - if nland == 0 { - f64::INFINITY - } else { - sums / nland as f64 - }, + avgz(sumh, nland), + avgz(sums, nland), + avgz(sumsed, newh.len()), + avgz(sumsed_land, nland), 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)*/; - max_slopes[posi] = if (old_h_i - old_b_i as f64) > sediment_thickness && kdsed > 0.0 { + (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/; + 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,29 +1131,35 @@ 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 { - 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; - // max_slopes[posi] = /*(mag_slope - max_slope) * */max_slopes[posi].max(kdsed); - /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi); - sums += mag_slope; */ - /* if kd_factor < 1.0 { - max_slopes[posi] /= kd_factor; - } else { - max_slopes[posi] *= kd_factor; - } */ - // max_slopes[posi] *= kd_factor; - nland += 1; - sumh += new_h_i; + 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; + // max_slopes[posi] = /*(mag_slope - max_slope) * */max_slopes[posi].max(kdsed); + /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi); + sums += mag_slope; */ + /* if kd_factor < 1.0 { + max_slopes[posi] /= kd_factor; + } else { + max_slopes[posi] *= kd_factor; + } */ + // 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; } - ntherm += 1; + 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; - sums += mag_slope; - // Just use the computed rate. - nland += 1; - sumh += new_h_i; + 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; - maxh = h[posi].max(maxh); + 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 { - f64::INFINITY - } else { - sumh / nland as f64 - }, - if nland == 0 { - f64::INFINITY - } else { - sums / nland as f64 - }, + avgz(sumh, nland), + avgz(sums, nland), + avgz(sumsed, newh.len()), + avgz(sumsed_land, nland), 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(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 neighbor_pass_height = h[pass.1 as usize]; + let pass_height = h(neighbor_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::>() .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::>() .into_boxed_slice(); // Stream power law erodability constant for fluvial erosion (bedrock) diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs index fe85b86312..5c2a670315 100644 --- a/world/src/sim/mod.rs +++ b/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]>, - basement: Box<[f32]>, + alt: Box<[Alt]>, + 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_small_steps = 8;//50;//8;//8;//8;//8;//8; // 8 + 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 = 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) + - /*exp_inverse_cdf*/(oheight.max(0.0).min(max_epsilon).abs())/*.powf(2.0)*/ * 0.2 * udist_2.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)*/; */ + (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 = 1e-6: normal-low (dyke [exposed]) + // kf = 2e-5: normal (dike [unexposed]) + // 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) { old_height(posi) } else { 0.0 }, + |posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.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| basement[posi], + |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32, + |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 (boundary_len, indirection, water_alt_pos, _) = get_lakes(&/*water_alt*/alt, &mut dh); + 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*/|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 basement_pre = gen_cdf.basement[posi]; + let alt_pre = gen_cdf.alt[posi] as f32; + 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]; diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs index 046cb05581..ced1fc93fa 100644 --- a/world/src/sim/util.rs +++ b/world/src/sim/util.rs @@ -247,12 +247,13 @@ pub fn uphill<'a>(dh: &'a [isize], posi: usize) -> impl Clone + Iterator bool + Sync) -> Box<[isize]> { +pub fn downhill(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() - .enumerate() - .map(|(posi, &nh)| { + (0..WORLD_SIZE.x * WORLD_SIZE.y).into_par_iter() + // .enumerate() + .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 {