diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs index dd231eac49..fcf7eacffc 100644 --- a/world/src/column/mod.rs +++ b/world/src/column/mod.rs @@ -436,13 +436,14 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { (temp < CONFIG.snow_temp || downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/ is_rocky { - sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? + // sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)? - // sim.get_interpolated(wpos, |chunk| chunk.alt)? + sim.get_interpolated(wpos, |chunk| chunk.alt)? } else { - sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? - // sim.get_interpolated(wpos, |chunk| chunk.alt)? + // sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? + sim.get_interpolated(wpos, |chunk| chunk.alt)? }; + let basement = sim.get_interpolated/*get_interpolated_monotone*/(wpos, |chunk| chunk.basement)?; // Find the average distance to each neighboring body of water. let mut river_count = 0.0f64; @@ -828,6 +829,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor); let alt = alt_ + riverless_alt_delta; + let basement = basement.min(alt); let rock = (sim.gen_ctx.small_nz.get( Vec3::new(wposf.x, wposf.y, alt as f64) @@ -917,6 +919,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { // For below desert humidity, we are always sand or rock, depending on altitude and // temperature. let ground = Rgb::lerp( + cliff, Rgb::lerp( dead_tundra, sand, @@ -924,10 +927,11 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { .div(CONFIG.desert_temp.sub(CONFIG.snow_temp)) .mul(/*4.5*/ 0.5), ), - cliff, - alt.sub(CONFIG.sea_level) + alt.sub(basement) + .mul(0.25) + /* alt.sub(CONFIG.sea_level) .sub(CONFIG.mountain_scale * 0.25) - .div(CONFIG.mountain_scale * 0.125), + .div(CONFIG.mountain_scale * 0.125), */ ); // From desert to forest humidity, we go from tundra to dirt to grass to moss to sand, // depending on temperature. @@ -1091,14 +1095,15 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { water_level, warp_factor, surface_color: Rgb::lerp( - sand, + Rgb::lerp(cliff, sand, alt.sub(basement).mul(0.25)), // Land ground, // Beach ((ocean_level - 1.0) / 2.0).max(0.0), ), - sub_surface_color: /*warm_grass*/dirt, - tree_density, + sub_surface_color: /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)), + // No growing directly on bedrock. + tree_density: if alt - basement < 2.0 { 0.0 } else { tree_density }, forest_kind: sim_chunk.forest_kind, close_structures: self.gen_close_structures(wpos), cave_xy, diff --git a/world/src/sim/erosion.rs b/world/src/sim/erosion.rs index f7ffc57119..be34953985 100644 --- a/world/src/sim/erosion.rs +++ b/world/src/sim/erosion.rs @@ -521,7 +521,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn> + Sync // sure anyone would be able to usefully tweak them on a per-map basis? Plus it's just a // hacky heuristic anyway. let center = 0.25; - let dmin = center - 0.15;//0.05; + let dmin = center - /*0.15;//0.05*/0.05; let dmax = center + 0.05;//0.05; let log_odds = |x: f64| logit(x) - logit(center); let rock_strength = logistic_cdf( @@ -592,6 +592,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn> + Sync fn erode( h: &mut [f32], b: &mut [f32], + wh: &mut [f32], is_done: &mut BitBox, done_val: bool, erosion_base: f32, @@ -625,7 +626,7 @@ fn erode( let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4; println!("dt={:?}", dt); // Landslide constant: ideally scaled to 10e-2 m / y^-1 - let l = /*200.0 * max_uplift as f64;*/10.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt; + let l = /*200.0 * max_uplift as f64;*/1.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt; // Net precipitation rate (m / year) let p = 1.0 * height_scale; // Stream power erosion constant (bedrock), in m^(1-2m) / year (times dt). @@ -638,12 +639,13 @@ fn erode( // 8e-6 * dt // 2e-5 * dt; // 2.244/2048*5*32/(250000/4)*10^6 - erosion_base as f64 + 2.244 / mmaxh as f64 * 5.0 * max_uplift as f64; + // ln(tan(30/360*2*pi))-ln(tan(6/360*2*pi))*1500 = 3378 + erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*/5.0/*3.75*/ * max_uplift as f64; // 1e-5 * dt; // (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 = k_fb * 1.0/* 2.0*//*4.0*/; + let k_fs = k_fb * /*1.0*//*2.0*/1.0/*4.0*/; let ((dh, indirection, newh, area), mut max_slopes) = rayon::join( || { let mut dh = downhill(h, |posi| is_ocean(posi)); @@ -677,11 +679,12 @@ fn erode( if posj == -1 { panic!("Disconnected lake!"); } + // wh for oceans is always 0. + wh[posi] = 0.0; // max_slopes[posi] = kd(posi); // Egress with no outgoing flows. } else { // *is_done.at(posi) = done_val; - nland += 1; let posj = posj as usize; let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64); @@ -694,8 +697,9 @@ fn erode( // 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[posi] as f64; + let old_b_i = b[posi] as f64; let uplift_i = uplift(posi) as f64; - let k = if (old_h_i - b[posi] as f64) > 1.0e-6 { + let k = if (old_h_i - old_b_i as f64) > 1.0e-6 { // Sediment k_fs } else { @@ -706,13 +710,30 @@ fn erode( let flux = k * (p * chunk_area * area[posi] as f64).sqrt() / neighbor_distance; // 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; // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||) - let mut new_h_i = (old_h_i + (uplift_i + flux * h_j)) / (1.0 + flux); + 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 { + new_h_i = (new_h_i + flux * h_j) / (1.0 + flux); + // 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 { + nland += 1; + sumh += new_h_i; + } + } + // 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; // 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. - let indirection_idx = indirection[posi]; + /* let indirection_idx = indirection[posi]; let is_lake_bottom = indirection_idx < 0; let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0; // Test the slope. @@ -736,12 +757,109 @@ fn erode( let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64); let neighbor_distance = (neighbor_coef * dxy).magnitude(); let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/; - let mag_slope = dz/*.abs()*/ / neighbor_distance; + let mag_slope = dz/*.abs()*/ / neighbor_distance; */ // If you're on the lake bottom and not right next to your neighbor, don't compute a // slope. if /* !is_lake_bottom */ /* !fake_neighbor */ true { + /* if + /* !is_lake_bottom && */ + mag_slope > max_slope { + // println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i); + // Thermal erosion says this can't happen, so we reduce dh_i to make the slope + // exactly max_slope. + // 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)); + let dz = (new_h_i - /*h_j*/h_k).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) * */kd(posi); + sums += mag_slope; */ + // let slope = dz.signum() * max_slope; + // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j; + // sums += max_slope; + } else { + // max_slopes[posi] = 0.0; + sums += mag_slope; + // Just use the computed rate. + } */ + h[posi] = new_h_i as f32; + // 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); + } + log::info!( + "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})", + maxh, + if nland == 0 { + f64::INFINITY + } else { + sumh / nland as f64 + }, + if nland == 0 { + f64::INFINITY + } else { + sums / nland as f64 + }, + ); + + // Apply thermal erosion. + maxh = 0.0; + sumh = 0.0; + for &posi in &*newh { + let posi = posi as usize; + let posj = dh[posi]; + if posj < 0 { + if posj == -1 { + panic!("Disconnected lake!"); + } + // max_slopes[posi] = kd(posi); + // Egress with no outgoing flows. + } else { + let posj = posj as usize; + let old_h_i = h[posi] as f64; + let old_b_i = b[posi] as f64; + // Find the water height for this chunk's receiver; we only apply thermal erosion + // for chunks above water. + 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. + if + /* !is_lake_bottom */ /* !fake_neighbor */ + wh_j < old_h_i { + let mut new_h_i = old_h_i; + let h_j = h[posj] as f64; + /* let indirection_idx = indirection[posi]; + let is_lake_bottom = indirection_idx < 0; + let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0; */ + // Test the slope. + let max_slope = max_slopes[posi] as f64; + // Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute + // uplift to other neighbors. + let (posk, h_k) = /* neighbors(posi) + .filter(|&posk| *is_done.at(posk) == done_val) + // .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk)) + .map(|posk| (posk, h[posk] as f64)) + // .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0) + .min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap()) + .unwrap_or((posj, h_j)); */ + (posj, h_j); + // .max(h_j); + let (posk, h_k) = if h_k < h_j { + (posk, h_k) + } else { + (posj, h_j) + }; + let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64); + let neighbor_distance = (neighbor_coef * dxy).magnitude(); + let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/; + let mag_slope = dz/*.abs()*/ / neighbor_distance; if /* !is_lake_bottom && */ mag_slope > max_slope { @@ -767,7 +885,7 @@ fn erode( } h[posi] = new_h_i as f32; // Make sure to update the basement as well! - b[posi] = (old_h_i + uplift_i).min(new_h_i) as f32; + b[posi] = old_b_i.min(new_h_i) as f32; sumh += new_h_i; } } @@ -775,7 +893,7 @@ fn erode( maxh = h[posi].max(maxh); } log::debug!( - "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})", + "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})", maxh, if nland == 0 { f64::INFINITY @@ -1255,6 +1373,7 @@ pub fn do_erosion( .map(|posi| oldb(posi)) .collect::>() .into_boxed_slice(); + let mut wh = vec![0.0; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice(); // TODO: Don't do this, maybe? // (To elaborate, maybe we should have varying uplift or compute it some other way). let uplift = (0..oldh_.len()) @@ -1282,8 +1401,15 @@ pub fn do_erosion( // 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 kdsed = /*1.5e-2*//*1e-4*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */; - let kd_bedrock = /*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */; + 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/*3.75*/ * max_uplift as f64) / dt; + let kd_bedrock = + /*1e-2*//*0.25e-2*/1e-2 * 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 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */ + /* * k_fb / 2e-5 */; let kd = |posi: usize| if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock }; // Hillslope diffusion coefficient for sediment. let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y]; @@ -1292,6 +1418,7 @@ pub fn do_erosion( erode( &mut h, &mut b, + &mut wh, &mut is_done, // The value to use to indicate that erosion is complete on a chunk. Should toggle // once per iteration, to avoid having to reset the bits, and start at true, since diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs index bddec66218..4a4fcdd709 100644 --- a/world/src/sim/mod.rs +++ b/world/src/sim/mod.rs @@ -190,7 +190,7 @@ impl WorldSim { let river_seed = RandomField::new(rng.gen()); let rock_lacunarity = 0.5/*HybridMulti::DEFAULT_LACUNARITY*/; let rock_strength_nz = /*Fbm*/HybridMulti_/*BasicMulti*//*Fbm*/::new() - .set_octaves(/*6*//*5*//*4*//*5*/4) + .set_octaves(/*6*//*5*//*4*//*5*//*4*/6) // persistence = lacunarity^(-(1.0 - fractal increment)) .set_persistence(/*0.9*/ /*2.0*//*1.5*//*HybridMulti::DEFAULT_LACUNARITY*/rock_lacunarity.powf(-(1.0 - 0.9))) // 256*32/2^4 @@ -199,7 +199,7 @@ impl WorldSim { // (0.5^(-(1.0-0.9)))^1/256/32*2^4*256*4 // (2^(-(1.0-0.9)))^4 // 16.0 - .set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (/*64.0*/256.0/*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */)) + .set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (64.0/*256.0*//*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */)) // .set_persistence(/*0.9*/ /*2.0*/0.67) // .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (2.0 * 32.0)) // .set_lacunarity(0.5) @@ -216,12 +216,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 = 32.0 * height_scale; + let max_erosion_per_delta_t = /*32.0*//*128.0*/128.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 = 150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200; - let n_small_steps = 0;//8;//8; // 8 + let n_steps = /*100*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200; + let n_small_steps = 8;//8;//8; // 8 // fractal dimension should be between 0 and 0.9999... // (lacunarity^octaves)^(-H) = persistence^(octaves) @@ -525,7 +525,7 @@ impl WorldSim { /* let erosion_factor = |x: f64| logistic_cdf(logistic_base * if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * log_odds(x))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs( ).powf(erosion_pow).mul(0.5))*/; */ 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)*/; +).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*/.powf(2.0); let uplift_fn = |posi| { if is_ocean_fn(posi) { @@ -606,7 +606,7 @@ impl WorldSim { // assert!(alt_main >= 0.0); let (bump_factor, bump_max) = if /*height < f32::EPSILON as f64 * 0.5/*false*/*/ - false { + true { ( /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/ (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64, @@ -615,11 +615,13 @@ impl WorldSim { } else { (0.0, 0.0) }; - let height = 1.0f64; + // let height = 1.0f64; // let height = 1.0 / 7.0f64; let height = height /* .mul(15.0 / 16.0) .add(1.0 / 16.0) */ + /* .mul(5.0 / 8.0) + .add(3.0 / 8.0) */ .mul(7.0 / 8.0) .add(1.0 / 8.0) .mul(max_erosion_per_delta_t) @@ -670,7 +672,7 @@ impl WorldSim { /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0) .mul(0.045)*/) }; - old_height_uniform(posi) + old_height_uniform(posi)/*.powf(2.0)*/ /* // 0.0 // -/*CONFIG.sea_level / CONFIG.mountain_scale*//* 0.75 */1.0 // ((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) as f32 @@ -1485,8 +1487,8 @@ impl SimChunk { _ => {} } - // No trees in the ocean or with zero humidity (currently) - let tree_density = if is_underwater { + // No trees in the ocean, with zero humidity (currently), or directly on bedrock. + let tree_density = if is_underwater/* || alt - basement.min(alt) < 2.0 */ { 0.0 } else { let tree_density = (gen_ctx.tree_nz.get((wposf.div(1024.0)).into_array()))