Changes to worldgen, adding more sedmient etc.

This commit is contained in:
Joshua Yanovski 2020-01-16 22:42:51 +01:00
parent ebe0d14eab
commit 1358f1dffa
11 changed files with 336 additions and 138 deletions

27
Cargo.lock generated
View File

@ -671,7 +671,7 @@ dependencies = [
"rand_core 0.5.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rand_os 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)",
"rand_xoshiro 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"serde 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
"serde_derive 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
"serde_json 1.0.42 (registry+https://github.com/rust-lang/crates.io-index)",
@ -1572,7 +1572,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"ahash 0.2.18 (registry+https://github.com/rust-lang/crates.io-index)",
"autocfg 0.1.7 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"serde 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
]
@ -1590,7 +1590,7 @@ version = "0.6.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"atom 0.3.5 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
@ -1703,7 +1703,7 @@ version = "0.1.16"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"byteorder 1.3.2 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
@ -2749,17 +2749,17 @@ dependencies = [
[[package]]
name = "rayon"
version = "1.2.1"
version = "1.3.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"crossbeam-deque 0.7.2 (registry+https://github.com/rust-lang/crates.io-index)",
"either 1.5.3 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon-core 1.6.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon-core 1.7.0 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "rayon-core"
version = "1.6.1"
version = "1.7.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"crossbeam-deque 0.7.2 (registry+https://github.com/rust-lang/crates.io-index)",
@ -3103,7 +3103,7 @@ dependencies = [
"arrayvec 0.4.12 (registry+https://github.com/rust-lang/crates.io-index)",
"hashbrown 0.6.3 (registry+https://github.com/rust-lang/crates.io-index)",
"mopa 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"shred-derive 0.6.1 (registry+https://github.com/rust-lang/crates.io-index)",
"smallvec 0.6.13 (registry+https://github.com/rust-lang/crates.io-index)",
]
@ -3206,7 +3206,7 @@ dependencies = [
"hashbrown 0.6.3 (registry+https://github.com/rust-lang/crates.io-index)",
"hibitset 0.6.2 (registry+https://github.com/rust-lang/crates.io-index)",
"log 0.4.8 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"serde 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
"shred 0.9.4 (registry+https://github.com/rust-lang/crates.io-index)",
"shrev 1.1.1 (registry+https://github.com/rust-lang/crates.io-index)",
@ -3599,7 +3599,7 @@ dependencies = [
"notify 5.0.0-pre.1 (registry+https://github.com/rust-lang/crates.io-index)",
"parking_lot 0.9.0 (registry+https://github.com/rust-lang/crates.io-index)",
"rand 0.7.2 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"ron 0.5.1 (registry+https://github.com/rust-lang/crates.io-index)",
"serde 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
"serde_derive 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
@ -3700,6 +3700,7 @@ dependencies = [
"faster 0.5.0 (git+https://github.com/AdamNiederer/faster.git?rev=6f99e0396e9992222bb33e8fd1e84347b410d9c0)",
"glsl-to-spirv 0.1.7 (registry+https://github.com/rust-lang/crates.io-index)",
"hashbrown 0.6.3 (registry+https://github.com/rust-lang/crates.io-index)",
"itertools 0.8.2 (registry+https://github.com/rust-lang/crates.io-index)",
"lazy_static 1.4.0 (registry+https://github.com/rust-lang/crates.io-index)",
"log 0.4.8 (registry+https://github.com/rust-lang/crates.io-index)",
"minifb 0.13.0 (git+https://github.com/emoon/rust_minifb.git)",
@ -3710,7 +3711,7 @@ dependencies = [
"pretty_env_logger 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rand 0.7.2 (registry+https://github.com/rust-lang/crates.io-index)",
"rand_chacha 0.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
"rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)",
"ron 0.5.1 (registry+https://github.com/rust-lang/crates.io-index)",
"roots 0.0.5 (registry+https://github.com/rust-lang/crates.io-index)",
"serde 1.0.102 (registry+https://github.com/rust-lang/crates.io-index)",
@ -4432,8 +4433,8 @@ dependencies = [
"checksum rand_xoshiro 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)" = "0e18c91676f670f6f0312764c759405f13afb98d5d73819840cf72a518487bff"
"checksum range-alloc 0.1.0 (registry+https://github.com/rust-lang/crates.io-index)" = "dd5927936723a9e8b715d37d7e4b390455087c4bdf25b9f702309460577b14f9"
"checksum raw-window-handle 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)" = "9db80d08d3ed847ce4fb3def46de0af4bfb6155bd09bd6eaf28b5ac72541c1f1"
"checksum rayon 1.2.1 (registry+https://github.com/rust-lang/crates.io-index)" = "43739f8831493b276363637423d3622d4bd6394ab6f0a9c4a552e208aeb7fddd"
"checksum rayon-core 1.6.1 (registry+https://github.com/rust-lang/crates.io-index)" = "f8bf17de6f23b05473c437eb958b9c850bfc8af0961fe17b4cc92d5a627b4791"
"checksum rayon 1.3.0 (registry+https://github.com/rust-lang/crates.io-index)" = "db6ce3297f9c85e16621bb8cca38a06779ffc31bb8184e1be4bed2be4678a098"
"checksum rayon-core 1.7.0 (registry+https://github.com/rust-lang/crates.io-index)" = "08a89b46efaf957e52b18062fb2f4660f8b8a4dde1807ca002690868ef2c85a9"
"checksum rdrand 0.4.0 (registry+https://github.com/rust-lang/crates.io-index)" = "678054eb77286b51581ba43620cc911abf02758c91f93f479767aed0f90458b2"
"checksum redox_syscall 0.1.56 (registry+https://github.com/rust-lang/crates.io-index)" = "2439c63f3f6139d1b57529d16bc3b8bb855230c8efcc5d3a896c8bea7c3b1e84"
"checksum redox_users 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)" = "4ecedbca3bf205f8d8f5c2b44d83cd0690e39ee84b951ed649e9f1841132b66d"

View File

@ -14,4 +14,4 @@ num_cpus = "1.10.1"
log = "0.4.8"
specs = "0.15.1"
vek = { version = "0.9.9", features = ["serde"] }
hashbrown = { version = "0.6.2", features = ["serde", "nightly"] }
hashbrown = { version = "0.6.2", features = ["rayon", "serde", "nightly"] }

View File

@ -20,10 +20,10 @@ ron = "0.5.1"
bincode = "1.2.0"
log = "0.4.8"
rand = "0.7.2"
rayon = "1.2.0"
rayon = "^1.3.0"
lazy_static = "1.4.0"
lz4-compress = "0.1.1"
hashbrown = { version = "0.6.2", features = ["serde", "nightly"] }
hashbrown = { version = "0.6.2", features = ["rayon", "serde", "nightly"] }
find_folder = "0.3.0"
parking_lot = "0.9.0"
crossbeam = "=0.7.2"

View File

@ -21,7 +21,7 @@ serde = "1.0.102"
serde_derive = "1.0.102"
rand = "0.7.2"
chrono = "0.4.9"
hashbrown = { version = "0.6.2", features = ["serde", "nightly"] }
hashbrown = { version = "0.6.2", features = ["rayon", "serde", "nightly"] }
crossbeam = "=0.7.2"
prometheus = "0.7"
prometheus-static-metric = "0.2"

View File

@ -59,7 +59,7 @@ treeculler = { git = "https://gitlab.com/yusdacra/treeculler.git" }
rodio = { git = "https://github.com/RustAudio/rodio", rev = "e5474a2"}
cpal = "0.10"
crossbeam = "=0.7.2"
hashbrown = { version = "0.6.2", features = ["serde", "nightly"] }
hashbrown = { version = "0.6.2", features = ["rayon", "serde", "nightly"] }
chrono = "0.4.9"
rust-argon2 = "0.5"
bincode = "1.2"

View File

@ -9,18 +9,19 @@ bincode = "1.2.0"
common = { package = "veloren-common", path = "../common" }
bitvec = "0.15.2"
faster = { git = "https://github.com/AdamNiederer/faster.git", rev = "6f99e0396e9992222bb33e8fd1e84347b410d9c0" }
itertools = "0.8.2"
vek = "0.9.9"
noise = { version = "0.6.0", default-features = false }
num = "0.2.0"
ordered-float = "1.0"
hashbrown = { version = "0.6.2", features = ["rayon", "serde"] }
hashbrown = { version = "0.6.2", features = ["rayon", "serde", "nightly"] }
lazy_static = "1.4.0"
log = "0.4.8"
rand = "0.7.2"
rand_chacha = "0.2.1"
arr_macro = "0.1.2"
packed_simd = "0.3.3"
rayon = "1.2.0"
rayon = "^1.3.0"
roots = "0.0.5"
serde = "1.0.102"
serde_derive = "1.0.102"

View File

@ -37,7 +37,7 @@ fn main() {
_map_file.push(map_file);
let world = World::generate(
1337,
5284,
WorldOpts {
seed_elements: false,
// world_file: sim::FileOpts::Load(_map_file),

View File

@ -443,8 +443,6 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
// sim.get_interpolated(wpos, |chunk| chunk.alt)?
};
let basement = alt
+ sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?;
// Find the average distance to each neighboring body of water.
let mut river_count = 0.0f64;
@ -831,6 +829,8 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
/* .mul((1.0 - chaos).min(1.0).max(0.3))
.mul(1.0 - humidity) */
/* .mul(32.0) */;
let basement = alt_
+ sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?;
let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor);
let alt = alt_ + riverless_alt_delta;

View File

@ -7,6 +7,7 @@ use arr_macro::arr;
use bitvec::prelude::{bitbox, bitvec, BitBox};
use common::{terrain::TerrainChunkSize, vol::RectVolSize};
use faster::*;
use itertools::izip;
use noise::{NoiseFn, Point3};
use num::{Float, FromPrimitive, One, Zero};
use ordered_float::NotNan;
@ -27,6 +28,24 @@ pub type Alt = f64;
pub type Compute = f64;
pub type Computex8 = [Compute; 8];
/// This is a fast approximation of powf. This should only be used when minor accuracy loss is acceptable.
#[inline(always)]
#[allow(unsafe_code)]
fn approx_powf(b: f32, e: f32) -> f32 {
unsafe {
let b = b as f64;
let e = e as f64;
union Swagger {
f: f64,
a: [i32; 2],
}
let mut b = Swagger { f: b };
b.a[1] = (e * (b.a[1] as f64 - 1072632447.0) + 1072632447.0) as i32;
b.a[0] = 0;
b.f as f32
}
}
/// 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]> {
@ -742,8 +761,8 @@ fn erode(
let k_df = 1.0e-4/*0.0*/;
// Debris flow area coefficient (m^(-2q)).
let q = 0.2;
let q_ = 1.5/*1.0*/;
let k_da = /*5.0*/2.5 * 16.0.powf(q);
let q_ = /*1.0*/1.5/*1.0*/;
let k_da = /*5.0*//*2.5*//*5.0*//*2.5*//*5.0*/2.5 * 4.0.powf(2.0 * q);
let nx = WORLD_SIZE.x;
let ny = WORLD_SIZE.y;
let dx = TerrainChunkSize::RECT_SIZE.x as f64/* * height_scale*//* / CONFIG.mountain_scale as f64*/;
@ -797,8 +816,9 @@ fn erode(
/* let epsilon_0 = 2.68e-4;
let alpha = 3e-2;
let epsilon_0_tot = epsilon_0 * dt; */
// Net precipitation rate (m / year)
let p = 1.0 * height_scale;
// Spatio-temporal variation in net precipitation rate ((m / year) / (m / year)) (ratio of
// precipitation rate at chunk relative to mean precipitation rate p₀).
let p = 1.0/*0.25*//* * height_scale*/;
/* let n = 2.4;// 1.0;//1.5;//2.4;//1.0;
let m = n * 0.5;// n * 0.4;// 0.96;// 0.4;//0.6;//0.96;//0.4; */
// Stream power erosion constant (bedrock), in m^(1-2m) / year (times dt).
@ -833,20 +853,22 @@ fn erode(
) = rayon::join(
|| {
let mut dh = downhill(
|posi| h[posi], /* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/
|posi| h[posi], /* + uplift(posi) as Alt*/
/* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/
|posi| {
is_ocean(posi)
&& h[posi]/* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/ <= 0.0
&& h[posi]/* + uplift(posi) as Alt*//* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/ <= 0.0
},
);
log::debug!("Computed downhill...");
let (boundary_len, indirection, newh, maxh) = get_lakes(
|posi| h[posi], /* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/
|posi| h[posi], /* + uplift(posi) as Alt*/
/* + a[posi].max(0.0)*//* + uplift(posi) as Alt*/
&mut dh,
);
log::debug!("Got lakes...");
let (mrec, mstack, mwrec) = get_multi_rec(
|posi| h[posi],
|posi| h[posi], /* + uplift(posi) as Alt*/
&dh,
&newh,
wh,
@ -870,7 +892,8 @@ fn erode(
max_slope
},
|| {
h.to_vec().into_boxed_slice()
h.to_vec() /*.map(|posi| h[mstack[posi]])*/
.into_boxed_slice()
/* rayon::join(
|| {
// Store the elevation at t
@ -892,9 +915,9 @@ fn erode(
// TODO: Make more principled.
let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan(); //1.0;
type SimdType = f64;
type MaskType = m64;
let simd_func = f64s/*f32s*/;
type SimdType = f32;
type MaskType = m32;
let simd_func = /*f64s*/f32s;
// Precompute factors for Stream Power Law.
let czero = </*Compute*/SimdType as Zero>::zero();
@ -971,8 +994,8 @@ fn erode(
// Higher rock strength tends to lead to higher curvature?
let kd_factor =
// 1.0;
(1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/;
let k_da = k_da / /*max_slope*/kd_factor;
(/*1.0 / */(max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0/* * q*/))/*.min(kdsed)*/;
let k_da = k_da /*/ /*max_slope*/*/ * kd_factor;// .powf(q/* / q_*/);
// let k_df = /*uplift_i*/0.05e-3 / (1.0 + k_da * /*chunk_area_pow*/(10_000.0).powf(q)) / max_slope.powf(q_);
// let k = (uplift_i - kf(posi) * dt * (p * chunk_area).powf(m_f(posi) as f64) * max_slope.powf(n_f(posi) as f64)).max(0.0) / (1.0 + k_da * chunk_area_pow) / max_slope.powf(q_);
@ -1026,8 +1049,8 @@ fn erode(
//
// (eliminating EΔt maintains the sign, but it's somewhat imprecise;
// we can address this later, e.g. by assigning a debris flow / fluvial erosion ratio).
let chunk_neutral_area = /*10.0e6*/1.0e6/*0.1e6*//*100.0 * 100.0*/; // 1 km^2 * (1000 m / km)^2 = 1e6 m^2
let k = (mwrec_kk * (uplift_i + max_uplift as f64 * g_i/* / p*/)) / (1.0 + k_da * (mwrec_kk * chunk_neutral_area).powf(q)) / max_slope.powf(q_);
let chunk_neutral_area = /*10.0e6*//*1.0e6*/0.1e6/*0.01e6*//*100.0 * 100.0*/; // 1 km^2 * (1000 m / km)^2 = 1e6 m^2
let k = (mwrec_kk * (uplift_i + max_uplift as f64 * g_i / p)) / (1.0 + k_da * (mwrec_kk * chunk_neutral_area).powf(q)) / max_slope.powf(q_);
// ∆p = ||chunk_i - rec_i,kk||
// k = k_df * Δt / (Δp)^(q_)
@ -1117,6 +1140,30 @@ fn erode(
let mut n_gs_stream_power_law = 0;
let max_n_gs_stream_power_law = /*199*/99;
let mut mstack_inv = vec![0usize; dh.len()];
let mut h_t_stack = vec![Zero::zero(); dh.len()];
let mut dh_stack = vec![0isize; dh.len()];
let mut h_stack = vec![Zero::zero(); dh.len()];
let mut b_stack = vec![Zero::zero(); dh.len()];
let mut area_stack = vec![Zero::zero(); dh.len()];
assert!(mstack.len() == dh.len());
assert!(b.len() == dh.len());
assert!(h_t.len() == dh.len());
let mut mstack_inv = &mut *mstack_inv;
mstack.iter().enumerate().for_each(|(stacki, &posi)| {
let posi = posi as usize;
mstack_inv[posi] = stacki;
dh_stack[stacki] = /*if dh_i < 0 {
dh_i
} else {
mstack[dh_i]
};*/dh[posi];
h_t_stack[stacki] = h_t[posi]/* + uplift(posi) as Alt*/;
h_stack[stacki] = h[posi];
b_stack[stacki] = b[posi];
area_stack[stacki] = area[posi];
});
while err > tol && n_gs_stream_power_law < max_n_gs_stream_power_law {
log::debug!("Stream power iteration #{:?}", n_gs_stream_power_law);
@ -1135,6 +1182,7 @@ fn erode(
prods_therm = 0.0;
prodscrit_therm = 0.0;
let start_time = Instant::now();
// Keep track of how many iterations we've gone to test for convergence.
n_gs_stream_power_law += 1;
@ -1144,10 +1192,11 @@ fn erode(
|| */
{
// guess/update the elevation at t+Δt (k)
h_p.par_iter_mut()
.zip_eq(/*h_*/ h.par_iter()) /*.zip(wh.par_iter())*/
.for_each(|((mut h_p, h_)/*, wh_*/)| {
*h_p = (*h_)/*.max(*wh_)*/ as Compute;
(&mut *h_p, &*h_stack)
.into_par_iter()
// .enumerate()
.for_each(|(/*stacki, */((mut h_p, h_)/*, wh_*/))| {
*h_p = (*h_)/*.max(*wh_)*/ as Compute/* + (0.5 - (stacki & 1) as Compute) * err as Compute*/;
});
/* hp.par_iter_mut().zip(h.par_iter()).for_each(|(mut hp, h)| {
*hp = *h as Compute;
@ -1166,7 +1215,11 @@ fn erode(
|| */
{
// calculate erosion/deposition of sediment at each node
deltah.par_iter_mut().enumerate().for_each(|(posi, mut deltah)| {
(&*mstack, &mut *deltah, &*h_t_stack, &*h_stack).into_par_iter()/*.enumerate()*/
.for_each(|(/*stacki, */(&posi, mut deltah, &h_t_i, &h_i))| {
let posi = posi as usize;
/* assert_eq!(mstack_inv[posi], stacki);
assert_eq!(h_t[posi], h_t_i); */
let uplift_i = uplift(posi) as Alt;
// h_j(t, FINAL) + U_j * Δt - h_j(t + Δt, k)
/* debug_assert!(((h[posi] + a[posi].max(0.0)) - h_[posi]).abs() <= 1.0e-6);
@ -1174,7 +1227,7 @@ fn erode(
let foo = (ht[posi] + uplift_i - h[posi] - (at[posi] - a[posi].max(0.0)/*.max(0.0)*/)) as Compute;
// let foo = (ht[posi] + at[posi] + uplift_i - (h[posi] + a[posi].max(0.0))) as Compute;
let bar = (at[posi]/* + uplift_i*/ - a[posi].max(0.0)/* - (h[posi] - ht[posi]*/)) as Compute; */
let delta = (/*ht[posi] + at[posi].max(0.0)*/h_t[posi] + uplift_i - /*h_*/h[posi]/*.max(wh[posi])*/) as Compute;
let delta = (/*ht[posi] + at[posi].max(0.0)*/h_t_i + uplift_i - /*h_*/h_i/*.max(wh[posi])*/) as Compute;
/* if (delta - (foo + bar)).abs() > 1.0e-6 {
println!("dh: {:?}, foo: {:?}, bar: {:?}, delta: {:?}, ht: {:?}, at: {:?}, h: {:?}, a: {:?}, h_: {:?}",
dh[posi], foo, bar, delta,
@ -1200,20 +1253,29 @@ fn erode(
)*/
},
);
log::debug!("(Done precomputation).");
log::debug!(
"(Done precomputation, time={:?}ms).",
start_time.elapsed().as_millis()
);
// sum the erosion in stack order
//
// After:
// deltah_i = Σ{j ∈ {i} upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_i(t + Δt, k))
/*newh.iter().rev()*/
mstack.into_iter().for_each(|&posi| {
let start_time = Instant::now();
izip!(&*mstack, &*dh_stack, /*&mut *deltah, */&h_t_stack, &*h_p).enumerate().for_each(|(stacki, (&posi, &posj, /*deltah_i, */&h_t_i, &h_p_i))| {
let posi = posi as usize;
let posj = dh[posi];
/* assert_eq!(mstack_inv[posi], stacki);
assert_eq!(posj, dh[posi]);
assert_eq!(h_t_i, h_t[posi]);
/* let posj = dh[posi]; */ */
let deltah_i = deltah[/*posi*/stacki];
if posj < 0 {
lake_silt[posi] = deltah[posi];
lake_silt[stacki] = deltah_i;// deltah[posi];
/* lake_sediment[posi] = deltah_sediment[posi];
lake_alluvium[posi] = deltah_alluvium[posi]; */
} else {
let posj = posj as usize;
let uplift_i = uplift(posi) as Alt;
/* if (deltah[posi] - (deltah_sediment[posi] + deltah_alluvium[posi])).abs() > 1.0e-2 {
println!("deltah_sediment: {:?}, deltah_alluvium: {:?}, deltah: {:?}, hp: {:?}, ap: {:?}, h_p: {:?}, ht: {:?}, at: {:?}, h: {:?}, a: {:?}, h_: {:?}",
@ -1223,15 +1285,15 @@ fn erode(
h[posi], a[posi], h_[posi]);
debug_assert_eq!(deltah[posi], deltah_sediment[posi] + deltah_alluvium[posi]);
} */
deltah[posi] -= ((/*ht[posi] + at[posi].max(0.0)*/h_t[posi] + uplift_i) as Compute - h_p[posi]);
let uphill_deltah_i = /*deltah[posi]*/deltah_i - ((/*ht[posi] + at[posi].max(0.0)*//*h_t[posi]*/h_t_i + uplift_i) as Compute - /*h_p[posi]*/h_p_i);
/* deltah_sediment[posi] -= ((ht[posi] + uplift_i - (at[posi] - ap[posi].max(0.0))) as Compute - hp[posi]);
deltah_alluvium[posi] -= ((at[posi]/* + uplift_i - (hp[posi] - ht[posi])*/) as Compute - ap[posi].max(0.0)); */
let lposi = lake_sill[posi];
if lposi == posi as isize {
if deltah[posi] <= 0.0 {
lake_silt[posi] = 0.0;
} else {
lake_silt[posi] = deltah[posi];
let lposi = lake_sill[stacki];
if lposi == stacki as isize {
if /*deltah[posi]*/uphill_deltah_i <= 0.0 {
lake_silt[stacki] = 0.0;
} else {
lake_silt[stacki] = /*deltah[posi]*/uphill_deltah_i;
}
/* if deltah_sediment[posi] <= 0.0 {
lake_sediment[posi] = 0.0;
@ -1244,11 +1306,12 @@ fn erode(
lake_alluvium[posi] = deltah_alluvium[posi];
} */
}
deltah[posi] += (/* ht[posi] + at[posi].max(0.0) */h_t[posi] + uplift_i) as Compute - h_p[posi];
// /*deltah[posi]*/deltah_i += (/* ht[posi] + at[posi].max(0.0) *//*h_t[posi]*/h_t_i + uplift_i) as Compute - /*h_p[posi]*/h_p_i;
let mwrec_i = &mwrec[posi];
mrec_downhill(&mrec, posi).for_each(|(k, posj)| {
deltah[posj] += deltah[posi] * mwrec_i[k]/*mwrec_i.extract(k)*/;
});
let stack_posj = mstack_inv[posj];
deltah[stack_posj] += deltah_i * mwrec_i[k];/*mwrec_i.extract(k)*/
})/*.sum::<Compute>()*/;
/* let posj = posj as usize;
deltah[posj] += deltah[posi]; */
@ -1265,7 +1328,10 @@ fn erode(
debug_assert_eq!(deltah[posi], deltah_sediment[posi] + deltah_alluvium[posi]);
} */
});
log::debug!("(Done sediment transport computation).");
log::debug!(
"(Done sediment transport computation, time={:?}ms).",
start_time.elapsed().as_millis()
);
// do ij=nn,1,-1
// ijk=stack(ij)
// ijr=rec(ijk)
@ -1285,18 +1351,26 @@ fn erode(
// endif
// enddo
elev.par_iter_mut().enumerate().for_each(|(posi, mut elev)| {
let start_time = Instant::now();
(&*mstack, &mut *elev, &*dh_stack, &*h_t_stack, &*area_stack, &*deltah, &*h_p, &*b_stack).into_par_iter()./*enumerate().*/for_each(|(/*stacki, */(&posi, mut elev, &dh_i, &h_t_i, &area_i, &deltah_i, &h_p_i, &b_i))| {
let posi = posi as usize;
/* assert_eq!(mstack_inv[posi], stacki);
assert_eq!(dh_i, dh[posi]);
assert_eq!(h_t_i, h_t[posi]);
assert_eq!(area_i, area[posi]);
assert_eq!(b_i, b[posi]); */
let uplift_i = uplift(posi) as Alt;
// let delta_a = a[posi] - at[posi];
if dh[posi] < 0 {
if dh_i < 0 {
// *elev = (ht[posi] + uplift_i - delta_a) as Compute;
*elev = (/*ht[posi] + at[posi].max(0.0)*/h_t[posi] + uplift_i) as Compute;
*elev = (/*ht[posi] + at[posi].max(0.0)*/h_t_i + uplift_i) as Compute;
} else {
// let old_h_after_uplift_i = (ht[posi] + uplift_i - delta_a) as Compute;
let old_h_after_uplift_i = (/*ht[posi] + at[posi].max(0.0)*/h_t[posi] + uplift_i) as Compute;
let area_i = area[posi] as Compute;
let old_h_after_uplift_i = (/*ht[posi] + at[posi].max(0.0)*/h_t_i + uplift_i) as Compute;
let area_i = area_i as Compute;
// debug_assert!(area_i > 0.0);
let uphill_silt_i = deltah[posi] - (old_h_after_uplift_i - h_p[posi]);
let uphill_silt_i = deltah_i - (old_h_after_uplift_i - h_p_i);
// let uphill_sediment_i = deltah_sediment[posi] - (old_h_after_uplift_i - hp[posi]);
/* let uphill_sediment_alluvium_i =
deltah_sediment[posi] - (ht[posi] + uplift_i - (at[posi] - ap[posi].max(0.0)) - hp[posi]) +
@ -1305,8 +1379,8 @@ fn erode(
deltah_sediment[posi] + deltah_alluvium[posi] -
2.0 * (old_h_after_uplift_i - (hp[posi] + ap[posi])); */
// let g_i = g(posi) as Compute;
let old_b_i = b[posi];
let sed = (h_t[posi] - old_b_i) as f64;
let old_b_i = b_i;
let sed = (h_t_i - old_b_i) as f64;
let g_i = if sed > sediment_thickness {
g_fs_mult_sed * g(posi) as Compute
} else {
@ -1316,7 +1390,7 @@ fn erode(
// actually was material to deposit. The current assumption is that as long as we
// are storing at most as much sediment as there actually was along the river, we
// are in the clear.
let g_i_ratio = (g_i / area_i)/*.min(1.0)*/;
let g_i_ratio = (g_i / (p * area_i))/*.min(1.0)*/;
// One side of nonlinear equation (23):
//
// h_i(t) + U_i * Δt + G / (p̃ * Ã_i) * Σ{j ∈ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_j(t + Δt, k))
@ -1334,7 +1408,10 @@ fn erode(
} */
}
});
log::debug!("(Done elevation estimation).");
log::debug!(
"(Done elevation estimation, time={:?}ms).",
start_time.elapsed().as_millis()
);
let start_time = Instant::now();
// TODO: Consider taking advantage of multi-receiver flow here.
@ -1350,29 +1427,37 @@ fn erode(
let mut simd_buf2 = [0.0; 8];
let mut simd_buf3 = [0.0; 8];
/*&*newh*/
mstack.into_iter().rev().for_each(|&posi| {
itertools::izip!(&*mstack, &*elev, /*&mut *wh, */&*b_stack/*, &mut *h_stack*/, &*h_t_stack, &*dh_stack, &*h_p)
.enumerate()
.rev()
.for_each(|(stacki, (&posi, &elev_i, /*wh_i, */&b_i, /*h_i, */&h_t_i, &dh_i, &h_p_i))| {
let mut iteration_error = 0.0;
let posi = posi as usize;
let old_elev_i = /*h*/elev[posi] as f64;
let old_wh_i = wh[posi];
let old_b_i = b[posi];
let old_ht_i = /*ht*/h_t[posi];
/* assert_eq!(mstack_inv[posi], stacki);
assert_eq!(dh_i, dh[posi]);
assert_eq!(h_t_i, h_t[posi]);
assert_eq!(b_i, b[posi]); */
let old_elev_i = /*h*//*elev[posi]*/elev_i as f64;
let old_wh_i = wh[posi]/*wh_i*/;
let old_b_i = /*b[posi]*/b_i;
let old_ht_i = /*ht*//*h_t[posi]*/h_t_i;
let sed = (old_ht_i - old_b_i) as f64;
let posj = dh[posi];
let posj = /*dh[posi]*/dh_i;
if posj < 0 {
if posj == -1 {
panic!("Disconnected lake!");
}
if /*ht*/h_t[posi] > 0.0 {
if /*ht*//*h_t[posi]*/h_t_i > 0.0 {
log::warn!("Ocean above zero?");
}
// wh for oceans is always at least min_erosion_height.
let uplift_i = uplift(posi) as Alt;
// wh[posi] = min_erosion_height.max(ht[posi] + uplift_i);
wh[posi] = min_erosion_height.max(/*ht[posi] + at[posi].max(0.0)*/h_t[posi] + uplift_i);
wh[posi] = min_erosion_height.max(/*ht[posi] + at[posi].max(0.0)*//*h_t[posi]*/h_t_i + uplift_i);
// debug_assert!(wh[posi].is_normal() || wh[posi] == 0.0);
lake_sill[posi] = posi as isize;
lake_water_volume[posi] = 0.0;
lake_sill[stacki] = posi as isize;
lake_water_volume[stacki] = 0.0;
// max_slopes[posi] = kd(posi);
// Egress with no outgoing flows.
} else {
@ -1397,7 +1482,7 @@ fn erode(
// let a_j = a[posj] as f64;
let wh_j = wh[posj] as f64;
// let old_a_i = a[posi] as f64;
let old_h_i = h[posi] as f64;
let old_h_i = /*h*/h_stack[/*posi*/stacki] as f64;
let mut new_h_i = /*old_elev_i*//*old_h_i + old_a_i.max(0.0)*/old_h_i/*h[posi] as f64*//* + uplift_i*/;
/* let mut df_part;
@ -1473,10 +1558,15 @@ fn erode(
let mut f = h0;
let mut df = 1.0;
mrec_downhill(&mrec, posi).for_each(|(kk, posj)| {
let posj_stack = mstack_inv[posj];
/* if posj_stack <= stacki {
println!("Huh? {:?} {:?}", posj_stack, stacki);
}
assert!(posj_stack > stacki); */
// This can happen in cases where receiver kk is neither uphill of
// nor downhill from posi's direct receiver.
if /*new_ht_i/*old_ht_i*/ >= (h_t[posj]/* + uplift(posj) as f64*/)*/old_elev_i /*>=*/> wh[posj] as f64/*h[posj]*//*h_j*/ {
let h_j = h[posj] as f64;
let h_j = h_stack[posj_stack] as f64;
let elev_j = h_j/* + a_j.max(0.0)*/;
/* let flux = /*k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance;*/k_fs_fact[kk] + k_df_fact[kk];
assert!(flux.is_normal() && flux.is_positive() || flux == 0.0);
@ -1494,6 +1584,7 @@ fn erode(
let omega2 = 0.875f64 / q_;/*if q_ < 0.5 { 0.875f64/* * q_*/ } else { 0.875f64 / q_ }*/;
let omega = omega1.max(omega2);
let tolp = 1.0e-3/*1.0e-4*/;
// let tolp = tol;
let mut errp = 2.0 * tolp;
// let h0 = old_elev_i + (new_h_i - old_h_i);
// let mut count = 0;
@ -1509,6 +1600,7 @@ fn erode(
let mut rec_heights = [0.0; 8];//f64s(0.0);//f64x8::splat(0.0);
let mut mask = [MaskType::new(false); 8];
mrec_downhill(&mrec, posi).for_each(|(kk, posj)| {
let posj_stack = mstack_inv[posj];
if old_elev_i > wh[posj] as f64 {
// k_fs_weights[kk] = k_fs_fact[kk] as SimdType;
// k_fs_weights[max] = k_fs_fact[kk] as SimdType;
@ -1517,7 +1609,7 @@ fn erode(
// k_df_weights[max] = k_df_fact[kk] as SimdType;
// /*k_df_weights = */k_df_weights.replace(max, k_df_fact[kk]/*.extract(kk)*/ as f64);
mask[kk] = MaskType::new(true);
rec_heights[kk] = h[posj] as SimdType;
rec_heights[kk] = h_stack[posj_stack] as SimdType;
// rec_heights[max] = h[posj] as SimdType;
// /*rec_heights = */rec_heights.replace(max, h[posj] as f64);
// max += 1;
@ -1546,7 +1638,8 @@ fn erode(
/* let czero = &czero[..max];
let n_weights = &n_weights[..max];
let q__weights = &q__weights[..max]; */
while errp > tolp {
// let mut count = 0;
while errp > tolp /*&& count < 4 */{
// count += 1;
/* if count > -1 {
println!("posi={:?} errp={:?} tolp={:?} h0={:?} new_h_i={:?} elev_j={:?} area={:?}", posi, errp, tolp, h0, new_h_i, elev_j, area[posi]);
@ -1666,14 +1759,16 @@ fn erode(
df += (n as SimdType * dh_fs_sample + q_ as SimdType * dh_df_sample) as f64;
}
} */
(0..8).for_each(|kk| {
izip!(&mask, &rec_heights, k_fs_fact, k_df_fact).for_each(|(&mask_kk, &rec_heights_kk, &k_fs_fact_kk, &k_df_fact_kk)| {
//if /*new_ht_i/*old_ht_i*/ > (h_t[posj]/* + uplift(posj) as f64*/)*/old_elev_i /*>=*/> wh[posj] as f64/*h[posj]*//*h_j*/ {
if mask[kk].test() {
let h_j = rec_heights[kk];
if /*mask[kk]*/mask_kk.test() {
let h_j = rec_heights_kk;//rec_heights[kk];
let elev_j = h_j/* + a_j.max(0.0)*/;
let dh = 0.0.max((new_h_i as SimdType - elev_j)/*.abs()*/);
let dh_fs_sample = k_fs_fact[kk] as /*f64*/SimdType * dh.powf(n as SimdType - 1.0);
let dh_df_sample = k_df_fact[kk] as /*f64*/SimdType * dh.powf(q_ as SimdType - 1.0);
let powf = |a: SimdType, b| a.powf(b);
// let powf = |a, b| approx_powf(a as f32, b as f32) as SimdType;
let dh_fs_sample = /*k_fs_fact[kk]*/k_fs_fact_kk as /*f64*/SimdType * powf(dh, n as SimdType - 1.0);
let dh_df_sample = /*k_df_fact[kk]*/k_df_fact_kk as /*f64*/SimdType * powf(dh, q_ as SimdType - 1.0);
// Want: h_i(t+Δt) = h0 - fact * (h_i(t+Δt) - h_j(t+Δt))^n
// Goal: h_i(t+Δt) - h0 + fact * (h_i(t+Δt) - h_j(t+Δt))^n = 0
f += ((dh_fs_sample + dh_df_sample) * dh) as f64;
@ -1703,12 +1798,40 @@ fn erode(
// h_i(t+∆t, k+1) = ...
new_h_i = new_h_i * (1.0 - omega) + hn * omega;
}
/* let mut f = new_h_i - h0;
let mut df = 1.0;
(0..8).for_each(|kk| {
//if /*new_ht_i/*old_ht_i*/ > (h_t[posj]/* + uplift(posj) as f64*/)*/old_elev_i /*>=*/> wh[posj] as f64/*h[posj]*//*h_j*/ {
if mask[kk].test() {
let h_j = rec_heights[kk];
let elev_j = h_j/* + a_j.max(0.0)*/;
let dh = 0.0.max((new_h_i as SimdType - elev_j)/*.abs()*/);
let powf = |a: SimdType, b| a.powf(b);
// let powf = |a, b| approx_powf(a as f32, b as f32) as SimdType;
let dh_fs_sample = k_fs_fact[kk] as /*f64*/SimdType * powf(dh, n as SimdType - 1.0);
let dh_df_sample = k_df_fact[kk] as /*f64*/SimdType * powf(dh, q_ as SimdType - 1.0);
// Want: h_i(t+Δt) = h0 - fact * (h_i(t+Δt) - h_j(t+Δt))^n
// Goal: h_i(t+Δt) - h0 + fact * (h_i(t+Δt) - h_j(t+Δt))^n = 0
f += ((dh_fs_sample + dh_df_sample) * dh) as f64;
// ∂h_i(t+Δt)/∂n = 1 + fact * n * (h_i(t+Δt) - h_j(t+Δt))^(n - 1)
df += (n as SimdType * dh_fs_sample + q_ as SimdType * dh_df_sample) as f64;
}
//}
});
// hn = h_i(t+Δt, k) - (h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) - h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)
let hn = new_h_i - f / df;
// errp = |(h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) - h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)|
errp = (hn - new_h_i).abs(); */
// iteration_error += errp;//.abs();
// assert!((h0 - new_h_i).abs() <= tolp);
// Correct new_h_i to keep it at or under h0.
/* if h0 < new_h_i {
println!("Huh? h0={:?}, new_h_i={:?},", h0, new_h_i);
// debug_assert!(new_h_i <= h0);
} */
new_h_i = h0.min(new_h_i);
// new_h_i = h0.min(new_h_i);
/* let omega = 0.875f64 / n;
let tolp = 1.0e-3;
let mut errp = 2.0 * tolp;
@ -1800,8 +1923,8 @@ fn erode(
}
}
} */
lake_sill[posi] = posi as isize;
lake_water_volume[posi] = 0.0;
lake_sill[stacki] = posi as isize;
lake_water_volume[stacki] = 0.0;
// If we dipped below the receiver's water level, set our height to the receiver's
// water level.
@ -1876,8 +1999,9 @@ fn erode(
// h_[posi] = old_elev_i;
// new_h_i = old_elev_i - df_part.max(0.0);
new_h_i = old_elev_i;
let lposj = lake_sill[posj];
lake_sill[posi] = lposj;
let posj_stack = mstack_inv[posj];
let lposj = lake_sill[posj_stack];
lake_sill[stacki] = lposj;
if lposj >= 0 {
let lposj = lposj as usize;
lake_water_volume[lposj] += (wh_j - old_elev_i) as Compute;
@ -1952,7 +2076,8 @@ fn erode(
sums += mag_slope;
// Just use the computed rate.
} */
h[posi] = new_h_i as Alt;
h_stack[/*(posi*/stacki] = new_h_i as Alt;
// println!("Delta: {:?} - {:?} = {:?}", new_h_i, h_p_i, new_h_i - h_p_i);
// a[posi] = df_part as Alt;
// Make sure to update the basement as well!
// b[posi] = (old_b_i + uplift_i).min(new_h_i) as f32;
@ -1962,7 +2087,7 @@ fn erode(
if compute_stats {
sumsed += sed;
// suma += a[posi];
let h_i = h[posi];
let h_i = h_stack[/*posi*/stacki];
if h_i > 0.0 {
minh = h_i.min(minh);
}
@ -1980,7 +2105,7 @@ fn erode(
debug_assert_eq!(h[posi] + a[posi].max(0.0), h_[posi]);
} */
// Add sum of squares of errors.
sum_err += (h[posi]/*.max(wh[posi])*/ as Compute/* + a[posi].max(0.0) as Compute*/ - /*hp*/h_p[posi] as Compute/* - ap[posi].max(0.0) as Compute*/).powi(2);
sum_err += (iteration_error + h_stack[/*posi*/stacki]/*.max(wh[posi])*/ as Compute/* + a[posi].max(0.0) as Compute*/ - /*hp*//*h_p[/*posi*/stacki]*/h_p_i as Compute/* - ap[posi].max(0.0) as Compute*/).powi(2);
});
log::debug!(
"(Done erosion computation, time={:?}ms)",
@ -1988,6 +2113,7 @@ fn erode(
);
err = (sum_err / /*newh*/mstack.len() as Compute).sqrt();
log::debug!("(RMSE: {:?})", err);
/* if max_g == 0.0 {
err = 0.0;
} */
@ -2000,6 +2126,14 @@ fn erode(
}
}
(&*mstack_inv, &mut *h)
.into_par_iter()
.enumerate()
.for_each(|(posi, (&stacki, h))| {
assert!(posi == mstack[stacki] as usize);
*h = h_stack[stacki];
});
//b=min(h,b)
// update the basement
@ -2098,21 +2232,25 @@ fn erode(
log::debug!("Done updating basement and applying soil production...");
// update the height to reflect sediment flux.
h.par_iter_mut().enumerate().for_each(|(posi, mut h)| {
let lposi = lake_sill[posi];
if lposi >= 0 {
let lposi = lposi as usize;
if lake_water_volume[lposi] > 0.0 {
// +max(0.d0,min(lake_sediment(lake_sill(ij)),lake_water_volume(lake_sill(ij))))/
// lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
*h += (0.0.max(
/*lake_sediment[lposi]*/ lake_silt[posi].min(lake_water_volume[lposi]),
) / lake_water_volume[lposi]
* (wh[posi]/* - a[posi].max(0.0)*/ - *h) as Compute)
as Alt;
(&mut *h, &*mstack_inv)
.into_par_iter()
.enumerate()
.for_each(|(posi, (mut h, &stacki))| {
let lposi = lake_sill[stacki];
if lposi >= 0 {
let lposi = lposi as usize;
if lake_water_volume[lposi] > 0.0 {
// +max(0.d0,min(lake_sediment(lake_sill(ij)),lake_water_volume(lake_sill(ij))))/
// lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
*h += (0.0.max(
/*lake_sediment[lposi]*/
lake_silt[stacki].min(lake_water_volume[lposi]),
) / lake_water_volume[lposi]
* (wh[posi]/* - a[posi].max(0.0)*/ - *h) as Compute)
as Alt;
}
}
}
});
});
// do ij=1,nn
// if (lake_sill(ij).ne.0) then
// if (lake_water_volume(lake_sill(ij)).gt.0.d0) h(ij)=h(ij) &
@ -3341,7 +3479,8 @@ pub fn do_erosion(
// Bedrock transport coefficients (diffusivity) in m^2 / year. For now, we set them all to be equal
// on land, but in theory we probably want to at least differentiate between soil, bedrock, and
// sediment.
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
let height_scale = 1.0 / 4.0; // 1.0 / CONFIG.mountain_scale as f64;
let time_scale = 1.0; //1.0 / 4.0; // 4.0.powf(-n)
let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
let k_fb = /*(erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*//*5.0*//*2.5*//*1.5*/4.0/*1.0*//*3.75*/ * max_uplift as f64) / dt;*/
@ -3350,7 +3489,7 @@ pub fn do_erosion(
/*1e-2*//*0.25e-2*/1e-2 / 1.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/* * k_fb / 2e-5 */;
let kdsed =
/*1.5e-2*//*1e-4*//*1.25e-2*//*1.5e-2 */1.5e-2/ 1.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/*1.5e-2*//*1e-4*//*1.25e-2*//*1.5e-2 */1.5e-2 / 1.0 * height_scale * height_scale / time_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/* * k_fb / 2e-5 */;
// let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
let n = |posi: usize| n[posi];

View File

@ -107,7 +107,7 @@ pub(crate) struct GenCtx {
pub fast_turb_y_nz: FastNoise,
pub town_gen: StructureGen2d,
// pub loc_gen: StructureGen2d,
pub river_seed: RandomField,
pub rock_strength_nz: Fbm,
pub uplift_nz: Worley,
@ -175,6 +175,7 @@ impl WorldSim {
let uplift_scale = /*512.0*//*256.0*/128.0;
let uplift_turb_scale = uplift_scale / 4.0/*32.0*//*64.0*/;
// NOTE: Changing order will significantly change WorldGen, so try not to!
let gen_ctx = GenCtx {
turb_x_nz: SuperSimplex::new().set_seed(rng.gen()),
turb_y_nz: SuperSimplex::new().set_seed(rng.gen()),
@ -271,9 +272,10 @@ impl WorldSim {
.set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * uplift_scale))
// .set_displacement(/*0.5*/0.0)
.set_displacement(/*0.5*/1.0)
.set_range_function(RangeFunction::Euclidean)
.set_range_function(RangeFunction::Euclidean),
// .enable_range(true),
// g_nz: RidgedMulti::new()
// loc_gen: StructureGen2d::new(rng.gen(), 2048, 1024),
};
let river_seed = &gen_ctx.river_seed;
@ -290,7 +292,7 @@ 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 = /*8.0*//*32.0*//*1.0*//*32.0*//*32.0*//*64.0*/16.0/*128.0*//*1.0*//*0.2 * /*100.0*/250.0*//*128.0*//*16.0*//*128.0*//*32.0*/ * height_scale;
let max_erosion_per_delta_t = /*8.0*//*32.0*//*1.0*//*32.0*//*32.0*//*16.0*//*64.0*//*32.0*/16.0/*128.0*//*1.0*//*0.2 * /*100.0*/250.0*//*128.0*//*16.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;
@ -692,6 +694,7 @@ impl WorldSim {
let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
let rock_strength_div_factor = /*8.0*/(2.0 * TerrainChunkSize::RECT_SIZE.x as f64) / 8.0;
let time_scale = 1.0; //4.0/*4.0*/;
let n_func = |posi| {
if is_ocean_fn(posi) {
return 1.0;
@ -734,12 +737,15 @@ impl WorldSim {
let theta_func = |posi| 0.4;
let kf_func = {
|posi| {
let m_i = (theta_func(posi) * n_func(posi)) as f64;
// let precip_mul = (0.25).powf(m);
if is_ocean_fn(posi) {
return 1.0e-4;
// return 2.0e-5;
// return 2.0e-6;
// return 2.0e-10;
// return 0.0;
// multiplied by height_scale^(2m) to account for change in area.
return 1.0e-4 * 4.0.powf(2.0 * m_i)/* / time_scale*/; // .powf(-(1.0 - 2.0 * m_i))/* * 4.0*/;
// return 2.0e-5;
// return 2.0e-6;
// return 2.0e-10;
// return 0.0;
}
let wposf = (uniform_idx_as_vec2(posi)
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@ -801,7 +807,8 @@ impl WorldSim {
// ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight_2) * (1.5e-4 - 2.0e-6) + 2.0e-6)
// ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight) * (1.5e-4 - 2.0e-6) + 2.0e-6)
// 2e-5
2.5e-6 * 16.0.powf(0.4) /* / 4.0 * 0.25 *//* * 4.0*/
// multiplied by height_scale^(2m) to account for change in area.
2.5e-6 * 4.0.powf(/*-(1.0 - 2.0 * m_i)*/ 2.0 * m_i) /* / time_scale*//* / 4.0 * 0.25 *//* * 4.0*/
// 2.9e-10
// ((1.0 - uheight) * (5e-5 - 2.9e-10) + 2.9e-10)
// ((1.0 - uheight) * (5e-5 - 2.9e-14) + 2.9e-14)
@ -809,8 +816,9 @@ impl WorldSim {
};
let kd_func = {
|posi| {
// let height_scale = 1.0 / 4.0.powf(-n_i);
if is_ocean_fn(posi) {
return 1.0e-2;
return 1.0e-2 * 4.0.powf(-2.0) * time_scale;
}
let wposf = (uniform_idx_as_vec2(posi)
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@ -841,8 +849,10 @@ impl WorldSim {
// kd = 1e-1: high (mountain, dike)
// kd = 1.5e-2: normal-high (plateau [fan sediment])
// kd = 1e-2: normal (plateau)
1.0e-2 * (1.0 / 16.0) // m_old^2 / y * (1 m_new / 4 m_old)^2
// (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2)
// multiplied by height_scale² to account for change in area, then divided by
// time_scale to account for lower dt.
1.0e-2 * 4.0.powf(-2.0) * time_scale // m_old^2 / y * (1 m_new / 4 m_old)^2
// (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2)
}
};
let g_func = |posi| {
@ -907,12 +917,42 @@ impl WorldSim {
// 2.0
// 0.0
1.0
// 4.0
// 1.0
// 1.5
};
let epsilon_0_func = |posi| {
let n_i = n_func(posi);
// height_scale is roughly [using Hack's Law with b = 2 and SPL without debris flow or
// hillslopes] equal to the ratio of the old to new area, to the power of -n_i.
let height_scale = 4.0.powf(-n_i);
if is_ocean_fn(posi) {
// marine: ε₀ = 2.078e-3
return 2.078e-3;
// divide by height scale, multiplied by time_scale, cancels out to 1; idea is that
// we are finishing in time τ = τ' * height_scale. We have production
// rate
//
// ∆P = ε₀ e^(-αH) Δt
// = ε₀ e^(-α' / height_scale * (H' * height_scale)) (Δt' * height_scale)
// = ε₀ e^(-α' H') (Δt' * height_scale)
//
// while the old production rate was
//
// ∆P' = ε₀' e^(-α'H') Δt'.
//
// BUT, we don't actually want the same production rate, but rather the same
// *relative* production rate, which means we actually want to multiply by
// height_scale again... this entails multiplying the right hand side by the
// production rate, which gets us
//
// ΔP = ΔP' * height_scale
// ΔP / height_scale = ΔP'
//
// so to equate them we need
//
// ε₀ e^(-α' H') (Δt' * height_scale) / height_scale = ε₀' e^(-α' H') Δt'
// ε₀ = ε₀'
return 2.078e-3 * height_scale/* * time_scale*/;
// return 5.0;
}
let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@ -954,9 +994,9 @@ impl WorldSim {
.max(-1.0)
.mul(0.5)
.add(0.5);
let center = /*0.25*/0.4;
let dmin = center - /*0.15;//0.05*/0.05;
let dmax = center + /*0.05*//*0.10*/0.05; //0.05;
let center = /*0.25*/0.4 / 4.0;
let dmin = center - /*0.15;//0.05*/0.05 / 4.0;
let dmax = center + /*0.05*//*0.10*/0.05 / 4.0; //0.05;
let log_odds = |x: f64| logit(x) - logit(center);
let ustrength = logistic_cdf(
1.0 * logit(rock_strength.min(1.0f64 - 1e-7).max(1e-7))
@ -969,13 +1009,26 @@ impl WorldSim {
// Point Reyes: ε₀ = 8.1e-5
// Nunnock River (fractured granite, least weathered?): ε₀ = 5.3e-5
// The stronger the rock, the lower the production rate of exposed bedrock.
// ((1.0 - ustrength) * (/*3.18e-4*/2.078e-3 - 5.3e-5) + 5.3e-5) as f32
0.0
// divide by height scale, then multiplied by time_scale, cancels out.
((1.0 - ustrength) * (/*3.18e-4*/2.078e-3 - 5.3e-5) + 5.3e-5) as f32 * height_scale
/* * time_scale*/
// 0.0
};
let alpha_func = |posi| {
let n_i = n_func(posi);
// height_scale is roughly [using Hack's Law with b = 2 and SPL without debris flow or
// hillslopes] equal to the ratio of the old to new area, to the power of -n_i.
let height_scale = 4.0.powf(-n_i);
if is_ocean_fn(posi) {
// marine: α = 3.7e-2
return 3.7e-2;
// divided by height_scale; idea is that since the final height itself will be
// the old height * height scale, and we take the rate as ε₀ * e^(-αH), to keep
// the rate of rate of change in soil production consistent we must divide H by
// height_scale.
//
// αH = α(H' * height_scale) = α'H'
// α = α' / height_scale
return 3.7e-2 / height_scale;
}
let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
.map(|e| e as f64);
@ -1016,9 +1069,9 @@ impl WorldSim {
.max(-1.0)
.mul(0.5)
.add(0.5);
let center = /*0.25*/0.4;
let dmin = center - /*0.15;//0.05*/0.05;
let dmax = center + /*0.05*//*0.10*/0.05; //0.05;
let center = /*0.25*/0.4 / 4.0;
let dmin = center - /*0.15;//0.05*/0.05 / 4.0;
let dmax = center + /*0.05*//*0.10*/0.05 / 4.0; //0.05;
let log_odds = |x: f64| logit(x) - logit(center);
let ustrength = logistic_cdf(
1.0 * logit(rock_strength.min(1.0f64 - 1e-7).max(1e-7))
@ -1031,7 +1084,8 @@ impl WorldSim {
// Nunnock river (fractured granite, least weathered?): α = 2e-3
// Point Reyes: α = 1.6e-2
// The stronger the rock, the faster the decline in soil production.
(ustrength * (4.2e-2 - 1.6e-2) + 1.6e-2) as f32
// divided by height_scale.
(ustrength * (4.2e-2 - 1.6e-2) + 1.6e-2) as f32 / height_scale
};
let uplift_fn = |posi| {
if is_ocean_fn(posi) {
@ -1664,7 +1718,6 @@ impl WorldSim {
/// Prepare the world for simulation
pub fn seed_elements(&mut self) {
let mut rng = self.rng.clone();
let random_loc = RandomField::new(rng.gen());
let cell_size = 16;
let grid_size = WORLD_SIZE / cell_size;

View File

@ -20,11 +20,15 @@ impl<V: Default> Default for SmallCache<V> {
fn default() -> Self {
Self {
index: [None; CACHE_LEN + 9],
data: arr![V::default(); /*521*//*137*/41], // TODO: Use CACHE_LEN
data: arr![V::default(); /*521*//*137*/41], // TODO: Use CACHE_LEN + 9
random: 1,
}
}
}
/* pub struct SmalllCacheKeys<'id, V> {
} */
impl<V: Default> SmallCache<V> {
pub fn get<F: FnOnce(Vec2<i32>) -> V>(&mut self, key: Vec2<i32>, f: F) -> &V {
let idx = calc_idx(key) % CACHE_LEN;