Mirrored_Repos/veloren
1
0
Fork 0
mirror of https://gitlab.com/veloren/veloren.git synced 2024-08-30 18:12:32 +00:00
Code Issues Packages Projects Releases Wiki Activity

Map saving, soil production, speedup attempts.

Browse Source
This commit is contained in:
Joshua Yanovski 2019-12-11 10:14:50 +01:00
parent 33b0063e2c
commit ee5d383c46
10 changed files with 682 additions and 171 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -33,6 +33,7 @@
*.log *.log
settings.ron settings.ron
run.sh run.sh
maps
screenshots screenshots
todo.txt todo.txt

7
server/src/lib.rs
View File

@ -46,7 +46,7 @@ use std::{
}; };
use uvth::{ThreadPool, ThreadPoolBuilder}; use uvth::{ThreadPool, ThreadPoolBuilder};
use vek::*; use vek::*;
use world::{sim::{WORLD_SIZE, WorldOpts}, World}; use world::{sim::{FileOpts, WORLD_SIZE, WorldOpts}, World};
const CLIENT_TIMEOUT: f64 = 20.0; // Seconds const CLIENT_TIMEOUT: f64 = 20.0; // Seconds
pub enum Event { pub enum Event {
@ -106,6 +106,11 @@ impl Server {
let world = World::generate(settings.world_seed, WorldOpts { let world = World::generate(settings.world_seed, WorldOpts {
seed_elements: true, seed_elements: true,
world_file: if let Some(ref file) = settings.map_file {
FileOpts::Load(file.clone())
} else {
FileOpts::Generate
},
..WorldOpts::default() ..WorldOpts::default()
}); });

5
server/src/settings.rs
View File

@ -15,6 +15,7 @@ pub struct ServerSettings {
//pub login_server: whatever //pub login_server: whatever
pub start_time: f64, pub start_time: f64,
pub admins: Vec<String>, pub admins: Vec<String>,
pub map_file: Option<PathBuf>,
} }
impl Default for ServerSettings { impl Default for ServerSettings {
@ -27,6 +28,7 @@ impl Default for ServerSettings {
server_description: "This is the best Veloren server.".to_owned(), server_description: "This is the best Veloren server.".to_owned(),
max_players: 100, max_players: 100,
start_time: 9.0 * 3600.0, start_time: 9.0 * 3600.0,
map_file: None,
admins: [ admins: [
"Pfau", "Pfau",
"zesterer", "zesterer",
@ -97,12 +99,13 @@ impl ServerSettings {
[127, 0, 0, 1], [127, 0, 0, 1],
pick_unused_port().expect("Failed to find unused port!"), pick_unused_port().expect("Failed to find unused port!"),
)), )),
world_seed: 1337,
server_name: "Singleplayer".to_owned(), server_name: "Singleplayer".to_owned(),
server_description: "Who needs friends anyway?".to_owned(), server_description: "Who needs friends anyway?".to_owned(),
max_players: 100, max_players: 100,
start_time: 9.0 * 3600.0, start_time: 9.0 * 3600.0,
admins: vec!["singleplayer".to_string()], // TODO: Let the player choose if they want to use admin commands or not admins: vec!["singleplayer".to_string()], // TODO: Let the player choose if they want to use admin commands or not
..
Self::load() // Fill in remaining fields from settings.ron.
} }
} }

2
world/Cargo.toml
View File

@ -5,6 +5,7 @@ authors = ["Joshua Barretto <joshua.s.barretto@gmail.com>"]
edition = "2018" edition = "2018"
[dependencies] [dependencies]
bincode = "1.2.0"
common = { package = "veloren-common", path = "../common" } common = { package = "veloren-common", path = "../common" }
bitvec = "0.15.2" bitvec = "0.15.2"
vek = "0.9.9" vek = "0.9.9"
@ -20,6 +21,7 @@ arr_macro = "0.1.2"
rayon = "1.2.0" rayon = "1.2.0"
roots = "0.0.5" roots = "0.0.5"
serde = "1.0.102" serde = "1.0.102"
serde_derive = "1.0.102"
ron = "0.5.1" ron = "0.5.1"
pretty_env_logger = "0.3.0" pretty_env_logger = "0.3.0"

15
world/examples/water.rs
View File

@ -1,9 +1,9 @@
use common::{terrain::TerrainChunkSize, vol::RectVolSize}; use common::{terrain::TerrainChunkSize, vol::RectVolSize};
// use self::Mode::*; // use self::Mode::*;
use std::{f32, f64}; use std::{f32, f64, path::PathBuf};
use vek::*; use vek::*;
use veloren_world::{ use veloren_world::{
sim::{RiverKind, WorldOpts, WORLD_SIZE}, sim::{self, RiverKind, WorldOpts, WORLD_SIZE},
util::Sampler, util::Sampler,
World, CONFIG, World, CONFIG,
}; };
@ -29,8 +29,17 @@ const H: usize = 1024;
fn main() { fn main() {
pretty_env_logger::init(); pretty_env_logger::init();
let map_file =
// "map_1575990726223.bin";
// "map_1575987666972.bin";
"map_1576046079066.bin";
let mut _map_file = PathBuf::from("./maps");
_map_file.push(map_file);
let world = World::generate(1337, WorldOpts { let world = World::generate(1337, WorldOpts {
seed_elements: false, seed_elements: false,
// world_file: sim::FileOpts::Load(_map_file),
world_file: sim::FileOpts::Save,
..WorldOpts::default() ..WorldOpts::default()
}); });
@ -102,7 +111,7 @@ fn main() {
let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32); let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32);
let downhill_pos = (downhill let downhill_pos = (downhill
.map(|downhill_pos| downhill_pos/*.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| e / sz as i32)*/) .map(|downhill_pos| downhill_pos/*.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| e / sz as i32)*/)
.unwrap_or(pos) .unwrap_or(pos + Vec2::new(1, 0))
- pos)/* * scale*/ - pos)/* * scale*/
+ pos; + pos;
let downhill_alt = sampler let downhill_alt = sampler

9
world/src/block/mod.rs
View File

@ -61,9 +61,17 @@ impl<'a> BlockGen<'a> {
{ {
let cliff_pos3d = Vec3::from(*cliff_pos); let cliff_pos3d = Vec3::from(*cliff_pos);
// Conservative range of height: [15.70, 49.33]
let height = (RandomField::new(seed + 1).get(cliff_pos3d) % 64) as f32 let height = (RandomField::new(seed + 1).get(cliff_pos3d) % 64) as f32
// [0, 63] / (1 + 3 * [0.12, 1.32]) + 3 =
// [0, 63] / (1 + [0.36, 3.96]) + 3 =
// [0, 63] / [1.36, 4.96] + 3 =
// [0, 63] / [1.36, 4.96] + 3 =
// (height min) [0, 0] + 3 = [3, 3]
// (height max) [12.70, 46.33] + 3 = [15.70, 49.33]
/ (1.0 + 3.0 * cliff_sample.chaos) / (1.0 + 3.0 * cliff_sample.chaos)
+ 3.0; + 3.0;
// COnservative range of radius: [8, 47]
let radius = RandomField::new(seed + 2).get(cliff_pos3d) % 48 + 8; let radius = RandomField::new(seed + 2).get(cliff_pos3d) % 48 + 8;
max_height.max( max_height.max(
@ -203,6 +211,7 @@ impl<'a> BlockGen<'a> {
world.sim().gen_ctx.fast_turb_x_nz.get(wposf.div(25.0)) as f32, world.sim().gen_ctx.fast_turb_x_nz.get(wposf.div(25.0)) as f32,
world.sim().gen_ctx.fast_turb_y_nz.get(wposf.div(25.0)) as f32, world.sim().gen_ctx.fast_turb_y_nz.get(wposf.div(25.0)) as f32,
) * 8.0; ) * 8.0;
// let turb = Lerp::lerp(0.0, turb, warp_factor);
let wpos_turb = Vec2::from(wpos).map(|e: i32| e as f32) + turb; let wpos_turb = Vec2::from(wpos).map(|e: i32| e as f32) + turb;
let cliff_height = Self::get_cliff_height( let cliff_height = Self::get_cliff_height(

6
world/src/column/mod.rs
View File

@ -583,7 +583,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
.gen_ctx .gen_ctx
.small_nz .small_nz
.get((wposf_turb.div(128.0)).into_array()) as f32) .get((wposf_turb.div(128.0)).into_array()) as f32)
.mul(24.0); .mul(4.0);
let alt_for_river = alt let alt_for_river = alt
+ if overlap_count == 0.0 { + if overlap_count == 0.0 {
@ -806,7 +806,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
} else { } else {
(false, alt_for_river, downhill_water_alt, 1.0) (false, alt_for_river, downhill_water_alt, 1.0)
}; };
let warp_factor = 0.0; // let warp_factor = 0.0;
let riverless_alt_delta = (sim let riverless_alt_delta = (sim
.gen_ctx .gen_ctx
@ -1117,7 +1117,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
.powf(2.0) .powf(2.0)
.neg() .neg()
.add(1.0) .add(1.0)
.mul((1.15 - chaos).min(1.0)) .mul((1.32 - chaos).min(1.0))
}; };
let cave_xy = cave_at(wposf); let cave_xy = cave_at(wposf);
let cave_alt = alt - 24.0 let cave_alt = alt - 24.0

356
world/src/sim/erosion.rs
View File

@ -9,7 +9,9 @@ use rayon::prelude::*;
use std::{ use std::{
cmp::{Ordering, Reverse}, cmp::{Ordering, Reverse},
collections::BinaryHeap, collections::BinaryHeap,
f32, f64, mem, u32, f32, f64, mem,
path::PathBuf,
u32,
}; };
use vek::*; use vek::*;
@ -495,9 +497,9 @@ pub fn get_rivers(
/// ///
/// TODO: See if allocating in advance is worthwhile. /// TODO: See if allocating in advance is worthwhile.
fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> { fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
const MIN_MAX_ANGLE: f64 = 15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI; let min_max_angle = (15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI).tan();
const MAX_MAX_ANGLE: f64 = 60.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI; let max_max_angle = (45.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI).tan();
const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE; let max_angle_range = max_max_angle - min_max_angle;
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64; let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
h.par_iter() h.par_iter()
.enumerate() .enumerate()
@ -505,7 +507,7 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64; let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64;
// Normalized to be between 6 and and 54 degrees. // Normalized to be between 6 and and 54 degrees.
let div_factor = /*32.0*//*16.0*//*64.0*//*256.0*/8.0/*8.0*//*1.0*//*4.0*/ * /*1.0*/16.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */; let div_factor = /*32.0*//*16.0*//*64.0*//*256.0*/8.0/*8.0*//*1.0*//*4.0*//* * /*1.0*/16.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */*/;
let rock_strength = rock_strength_nz let rock_strength = rock_strength_nz
.get([ .get([
wposf.x, /* / div_factor*/ wposf.x, /* / div_factor*/
@ -548,7 +550,7 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
+ 1.0 * log_odds((wposz / CONFIG.mountain_scale as f64).abs().min(dmax).max(dmin)), + 1.0 * log_odds((wposz / CONFIG.mountain_scale as f64).abs().min(dmax).max(dmin)),
); );
// let rock_strength = 0.5; // let rock_strength = 0.5;
let max_slope = (rock_strength * MAX_ANGLE_RANGE + MIN_MAX_ANGLE).tan(); let max_slope = rock_strength * max_angle_range + min_max_angle;
// let max_slope = /*30.0.to_radians().tan();*/3.0.sqrt() / 3.0; // let max_slope = /*30.0.to_radians().tan();*/3.0.sqrt() / 3.0;
max_slope max_slope
}) })
@ -607,12 +609,22 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
/// ///
/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90 /// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90
/// ///
/// The approximate equation for soil production function (predictng the rate at which bedrock
/// turns into soil, increasing the distance between the basement and altitude) is taken from
/// equation (11) from [2]. This (among numerous other sources) also includes at least one
/// prediction that hillslope diffusion should be nonlinear, which we sort of attempt to
/// approximate.
///
/// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al.. /// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al..
/// Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion. /// Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
/// Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175. /// Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175.
/// ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩ /// ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
/// ///
/// /// [2] William E. Dietrich, Dino G. Bellugi, Leonard S. Sklar, Jonathan D. Stock
/// Geomorphic Transport Laws for Predicting Landscape Form and Dynamics.
/// Prediction in Geomorphology, Geophysical Monograph 135.
/// Copyright 2003 by the American Geophysical Union
/// 10.1029/135GM09
fn erode( fn erode(
h: &mut [Alt], h: &mut [Alt],
b: &mut [Alt], b: &mut [Alt],
@ -626,14 +638,17 @@ fn erode(
_seed: &RandomField, _seed: &RandomField,
rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync), rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
uplift: impl Fn(usize) -> f32 + Sync, uplift: impl Fn(usize) -> f32 + Sync,
kf: impl Fn(usize) -> f32, n_f: impl Fn(usize) -> f32 + Sync,
kd: impl Fn(usize) -> f32, m_f: impl Fn(usize) -> f32 + Sync,
kf: impl Fn(usize) -> f64 + Sync,
kd: impl Fn(usize) -> f64,
g: impl Fn(usize) -> f32 + Sync, g: impl Fn(usize) -> f32 + Sync,
is_ocean: impl Fn(usize) -> bool + Sync, is_ocean: impl Fn(usize) -> bool + Sync,
) { ) {
let compute_stats = true; let compute_stats = true;
log::debug!("Done draining..."); log::debug!("Done draining...");
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64; let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
let min_erosion_height = 0.0;//-f64::INFINITY as Alt;
let mmaxh = CONFIG.mountain_scale as f64 * height_scale; let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
// Minimum sediment thickness before we treat erosion as sediment based. // Minimum sediment thickness before we treat erosion as sediment based.
let sediment_thickness = 1.0; let sediment_thickness = 1.0;
@ -656,11 +671,58 @@ fn erode(
// 5e-7 // 5e-7
let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / /*5.010e-4*/1e-3; let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / /*5.010e-4*/1e-3;
println!("dt={:?}", dt); println!("dt={:?}", dt);
let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
let chunk_area = neighbor_coef.x * neighbor_coef.y;
let min_length = neighbor_coef.reduce_partial_min();
let max_stable = /*max_slope * */min_length * min_length / (dt/* / 2.0*/);//1.0/* + /*max_uplift as f64 / dt*/sed / dt*/;
// Landslide constant: ideally scaled to 10e-2 m / y^-1 // Landslide constant: ideally scaled to 10e-2 m / y^-1
let l = /*200.0 * max_uplift as f64;*/1.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);
let l_tot = l * dt;
// ε₀ = 0.000268 m/y, α = 0.03 (1/m). This is part of the soil production approximate
// equation:
//
// -∂z_b / ∂t = ε₀ * e^(-αH)
//
// where
// z_b is the elevation of the soil-bedrock interface (i.e. the basement),
// ε₀ is the production rate of exposed bedrock (H = 0),
// H is the soil thickness normal to the ground surface,
// and α is a parameter (units of 1 / length).
//
// Note that normal depth at i, for us, will be interpreted as the soil depth vector,
// sed_i = (0, h_i - b_i),
// projected onto the bedrock slope vector,
// bedrock_surface_i = (||p_i - p_j||, b_i - b_j),
// yielding the soil depth vector
// H_i = sed_i - sed_i ⋅ bedrock_surface_i / (bedrock_surface_i ⋅ bedrock_surface_i) * bedrock_surface_i
//
// = (0, h_i - b_i) -
// (0 * ||p_i - p_j|| + (h_i - b_i) * (b_i - b_j)) / (||p_i - p_j||^2 + (b_i - b_j)^2)
// * (||p_i - p_j||, b_i - b_j)
// = (0, h_i - b_i) -
// ((h_i - b_i) * (b_i - b_j)) / (||p_i - p_j||^2 + (b_i - b_j)^2)
// * (||p_i - p_j||, b_i - b_j)
// = (h_i - b_i) *
// ((0, 1) - (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2) * (||p_i - p_j||, b_i - b_j))
// H_i_fact = (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2)
// H_i = (h_i - b_i) * (((0, 1) - H_i_fact * (||p_i - p_j||, b_i - b_j)))
// = (h_i - b_i) * (-H_i_fact * ||p_i - p_j||, 1 - H_i_fact * (b_i - b_j))
// ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (b_i - b_j))^2))
//
// where j is i's receiver and ||p_i - p_j|| is the horizontal displacement between them. The
// idea here is that we first compute the hypotenuse between the horizontal and vertical
// displacement of bedrock (getting the horizontal component of the triangle), and then this is
// taken as one of the non-hypotenuse sides of the triangle whose other non-hypotenuse side is
// the normal height H_i, while their square adds up to the vertical displacement (h_i - b_i).
// If h and b have different slopes, this may not work completely correctly, but this is
// probably fine as an approximation.
let epsilon_0 = 2.68e-4;
let alpha = 3e-2;
let epsilon_0_tot = epsilon_0 * dt;
// Net precipitation rate (m / year) // Net precipitation rate (m / year)
let p = 1.0 * height_scale; let p = 1.0 * height_scale;
let m = 0.4; /* 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). // Stream power erosion constant (bedrock), in m^(1-2m) / year (times dt).
let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64; let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
// 2.244*(5.010e-4)/512*5- (1.097e-5) // 2.244*(5.010e-4)/512*5- (1.097e-5)
@ -711,9 +773,38 @@ fn erode(
); );
assert!(h.len() == dh.len() && dh.len() == area.len()); assert!(h.len() == dh.len() && dh.len() == area.len());
// Precompute factors for Stream Power Law.
let k_fact = dh.par_iter().enumerate()
.map(|(posi, &posj)| {
if posj < 0 {
// Egress with no outgoing flows, no stream power.
0.0
} else {
let posj = posj as usize;
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
let neighbor_distance = (neighbor_coef * dxy).magnitude();
let old_b_i = b[posi];
let sed = (ht[posi] - old_b_i) as f64;
let k = if sed > sediment_thickness {
// Sediment
// k_fs
k_fs_mult * kf(posi)
} else {
// Bedrock
// k_fb
kf(posi)
} * dt;
let n = n_f(posi) as f64;
let m = m_f(posi) as f64;
k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance.powf(n)
}
})
.collect::<Vec<f64>>();
log::info!("Computed stream power factors...");
// max angle of slope depends on rock strength, which is computed from noise function. // max angle of slope depends on rock strength, which is computed from noise function.
let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
let chunk_area = neighbor_coef.x * neighbor_coef.y;
// TODO: Make more principled. // TODO: Make more principled.
let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();//1.0; let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();//1.0;
@ -727,6 +818,7 @@ fn erode(
let mut err = 2.0 * tol; let mut err = 2.0 * tol;
// Some variables for tracking statistics, currently only for debugging purposes. // Some variables for tracking statistics, currently only for debugging purposes.
let mut minh = f64::INFINITY as Alt;
let mut maxh = 0.0; let mut maxh = 0.0;
let mut nland = 0usize; let mut nland = 0usize;
let mut sums = 0.0; let mut sums = 0.0;
@ -747,6 +839,7 @@ fn erode(
// Reset statistics in each loop. // Reset statistics in each loop.
maxh = 0.0; maxh = 0.0;
minh = f64::INFINITY as Alt;
nland = 0usize; nland = 0usize;
sums = 0.0; sums = 0.0;
sumh = 0.0; sumh = 0.0;
@ -774,6 +867,9 @@ fn erode(
); );
// sum the erosion in stack order // sum the erosion in stack order
//
// After:
// deltah_i = Σ{j ∈ {i} ∪ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_j(t + Δt, k))
for &posi in newh.iter().rev() { for &posi in newh.iter().rev() {
let posi = posi as usize; let posi = posi as usize;
let posj = dh[posi]; let posj = dh[posi];
@ -820,12 +916,19 @@ fn erode(
} else { } else {
let uplift_i = uplift(posi) as Alt; let uplift_i = uplift(posi) as Alt;
assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.0); assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.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))
//
// where
//
// Ã_i = A_i / (∆x∆y) = N_i, number of cells upstream of cell i.
*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; *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;
} }
}); });
// Iterate in ascending height order. // Iterate in ascending height order.
let mut sum_err = 0.0f64; let mut sum_err = 0.0 as Compute;
for &posi in &*newh { for &posi in &*newh {
let posi = posi as usize; let posi = posi as usize;
let old_h_i = /*h*/elev[posi] as f64; let old_h_i = /*h*/elev[posi] as f64;
@ -837,8 +940,11 @@ fn erode(
if posj == -1 { if posj == -1 {
panic!("Disconnected lake!"); panic!("Disconnected lake!");
} }
// wh for oceans is always 0. if ht[posi] > 0.0 {
wh[posi] = 0.0; log::warn!("Ocean above zero?");
}
// wh for oceans is always at least min_erosion_height.
wh[posi] = min_erosion_height.max(ht[posi]);
lake_sill[posi] = posi as isize; lake_sill[posi] = posi as isize;
lake_water_volume[posi] = 0.0; lake_water_volume[posi] = 0.0;
// max_slopes[posi] = kd(posi); // max_slopes[posi] = kd(posi);
@ -846,11 +952,11 @@ fn erode(
} else { } else {
// *is_done.at(posi) = done_val; // *is_done.at(posi) = done_val;
let posj = posj as usize; let posj = posj as usize;
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64); // let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
// Has an outgoing flow edge (posi, posj). // Has an outgoing flow edge (posi, posj).
// flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1 // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
let neighbor_distance = (neighbor_coef * dxy).magnitude(); // let neighbor_distance = (neighbor_coef * dxy).magnitude();
// Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5, // Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
// we have meters^(1) / meters^(1), so they should cancel out. Otherwise, we would // we have meters^(1) / meters^(1), so they should cancel out. Otherwise, we would
// want to multiply neighbor_distance by height_scale and area[posi] by // want to multiply neighbor_distance by height_scale and area[posi] by
@ -863,34 +969,82 @@ fn erode(
// Therefore, we can rely on wh being set to the water height for that node. // Therefore, we can rely on wh being set to the water height for that node.
let h_j = h[posj] as f64; let h_j = h[posj] as f64;
let wh_j = wh[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*/h[posi] as f64/* + uplift_i*/;
// Only perform erosion if we are above the water level of the previous node. // Only perform erosion if we are above the water level of the previous node.
if new_h_i > wh_j { if old_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||) // 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 sed > sediment_thickness { /* let k = if sed > sediment_thickness {
// Sediment // Sediment
// k_fs // k_fs
k_fs_mult * kf(posi) as f64 k_fs_mult * kf(posi)
} else { } else {
// Bedrock // Bedrock
// k_fb // k_fb
kf(posi) as f64 kf(posi)
} * dt; } * dt;
// let k = k * uplift_i / max_uplift as f64; // let k = k * uplift_i / max_uplift as f64;
let flux = k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance; let n = n_f(posi) as f64;
assert!(flux.is_normal() && flux.is_positive()); let m = m_f(posi) as f64; */
new_h_i = (new_h_i + flux * h_j) / (1.0 + flux); let n = n_f(posi) as f64;
if /*n == 1.0*/(n - 1.0).abs() <= 1.0e-3/*f64::EPSILON*/ {
let flux = /*k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance;*/k_fact[posi];
assert!(flux.is_normal() && flux.is_positive() || flux == 0.0);
new_h_i = (/*new_h_i*/old_h_i + flux * h_j) / (1.0 + flux);
} else {
// Local Newton-Raphson
let omega = 0.875f64 / n;
let tolp = 1.0e-3;
let mut errp = 2.0 * tolp;
let h0 = old_h_i;
let fact = k_fact[posi];// k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance.powf(n);
while errp > tolp {
let mut f = new_h_i - h0;
let mut df = 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 += fact * 0.0.max(new_h_i - h_j).powf(n);
// ∂h_i(t+Δt)/∂n = 1 + fact * n * (h_i(t+Δt) - h_j(t+Δt))^(n - 1)
df += fact * n * 0.0.max(new_h_i - h_j).powf(n - 1.0);
// 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();
// h_i(t+∆t, k+1) = ...
new_h_i = new_h_i * (1.0 - omega) + hn * omega;
}
/* omega=0.875d0/n
tolp=1.d-3
errp=2.d0*tolp
h0=elev(ijk)
do while (errp.gt.tolp)
f=h(ijk)-h0
df=1.d0
if (ht(ijk).gt.ht(ijr)) then
fact = kfint(ijk)*dt*a(ijk)**m/length(ijk)**n
f=f+fact*max(0.d0,h(ijk)-h(ijr))**n
df=df+fact*n*max(0.d0,h(ijk)-h(ijr))**(n-1.d0)
endif
hn=h(ijk)-f/df
errp=abs(hn-h(ijk))
h(ijk)=h(ijk)*(1.d0-omega)+hn*omega
enddo */
}
lake_sill[posi] = posi as isize; lake_sill[posi] = posi as isize;
lake_water_volume[posi] = 0.0; lake_water_volume[posi] = 0.0;
/* // Thermal erosion (landslide) /* // 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 dz = (new_h_i - /*h_j*//*h_k*//*wh_j*/h_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
let mag_slope = dz/*.abs()*/ / neighbor_distance; let mag_slope = dz/*.abs()*/ / neighbor_distance;
let max_slope = max_slopes[posi] as f64; let max_slope = max_slopes[posi] as f64;
if mag_slope > max_slope { if mag_slope > max_slope {
let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/; 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)); // new_h_i = (ht[posi] as f64 + l_tot * (mag_slope - max_slope));
if compute_stats && new_h_i > wh_j { // new_h_i = new_h_i - l_tot * (mag_slope - max_slope);
// new_h_i = new_h_i - l_tot * (mag_slope - max_slope);
// new_h_i = hp[posj] + dh;
new_h_i = /*old_h_i.max*/(/*wh_j*//*ht[posi] as Compute*//*h_j*/hp[posj]/*ht[posj] as Compute*/ + dh).max(new_h_i - l_tot * (mag_slope - max_slope));
if compute_stats/* && new_h_i > wh_j*/ {
ntherm += 1; ntherm += 1;
} }
} */ } */
@ -900,7 +1054,9 @@ fn erode(
if new_h_i <= wh_j { if new_h_i <= wh_j {
new_h_i = wh_j; new_h_i = wh_j;
} else { } else {
if compute_stats { if compute_stats && new_h_i > 0.0 {
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
let neighbor_distance = (neighbor_coef * dxy).magnitude();
let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/; 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 mag_slope = dz/*.abs()*/ / neighbor_distance;
@ -911,6 +1067,7 @@ fn erode(
} }
} }
} else { } else {
new_h_i = old_h_i;
let lposj = lake_sill[posj]; let lposj = lake_sill[posj];
lake_sill[posi] = lposj; lake_sill[posi] = lposj;
if lposj >= 0 { if lposj >= 0 {
@ -963,7 +1120,7 @@ fn erode(
// max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE // 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. // 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*/; 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 - l_tot * (mag_slope - max_slope));
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/; 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; let slope = dz/*.abs()*/ / neighbor_distance;
sums += slope; sums += slope;
@ -984,8 +1141,12 @@ fn erode(
} }
// *is_done.at(posi) = done_val; // *is_done.at(posi) = done_val;
if compute_stats { if compute_stats {
maxh = h[posi].max(maxh);
sumsed += sed; sumsed += sed;
let h_i = h[posi];
if h_i > 0.0 {
minh = h_i.min(minh);
}
maxh = h_i.max(maxh);
} }
// Add sum of squares of errors. // Add sum of squares of errors.
@ -1004,12 +1165,58 @@ fn erode(
//b=min(h,b) //b=min(h,b)
// update the basement // update the basement
b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| { //
// NOTE: Despite this not quite applying since basement order and height order differ, we once
// again borrow the implicit FastScape stack order. If this becomes a problem we can easily
// compute a separate stack order just for basement.
for &posi in &*newh {
let posi = posi as usize;
let old_b_i = b[posi];
let h_i = h[posi];
let uplift_i = uplift(posi) as Alt;
// First, add uplift...
let mut new_b_i = (old_b_i + uplift_i).min(h_i);
let posj = dh[posi];
// Sediment height normal to bedrock. NOTE: Currently we can actually have sedment and
// bedrock slope at different heights, meaning there's no uniform slope. There are
// probably more correct ways to account for this, such as averaging, integrating, or doing
// things by mass / volume instead of height, but for now we use the time-honored
// technique of ignoring the problem.
let h_normal = if posj < 0 {
// Egress with no outgoing flows; for now, we assume this means normal and vertical
// coincide.
(h_i - new_b_i) as f64
} else {
let posj = posj as usize;
let b_j = b[posj];
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
let neighbor_distance_squared = (neighbor_coef * dxy).magnitude_squared();
let vertical_sed = (h_i - new_b_i) as f64;
let db = (new_b_i - b_j) as f64;
// H_i_fact = (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2)
let h_i_fact = db / (neighbor_distance_squared + db * db);
let h_i_vertical = 1.0 - h_i_fact * db;
// ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (b_i - b_j))^2))
vertical_sed * (h_i_fact * h_i_fact * neighbor_distance_squared + h_i_vertical * h_i_vertical).sqrt()
};
// Rate of sediment production: -∂z_b / ∂t = ε₀ * e^(-αH)
let p_i = epsilon_0_tot * f64::exp(-alpha * h_normal);
// println!("h_normal = {:?}, p_i = {:?}", h_normal, p_i);
new_b_i -= p_i as Alt;
b[posi] = new_b_i;
}
log::info!("Done updating basement and applying soil production...");
/* b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| {
let old_b_i = *b; let old_b_i = *b;
let uplift_i = uplift(posi) as Alt; let uplift_i = uplift(posi) as Alt;
*b = (old_b_i + uplift_i).min(*h); *b = (old_b_i + uplift_i).min(*h);
}); }); */
// update the height to reflect sediment flux. // update the height to reflect sediment flux.
h.par_iter_mut().enumerate().for_each(|(posi, mut h)| { h.par_iter_mut().enumerate().for_each(|(posi, mut h)| {
@ -1035,11 +1242,12 @@ fn erode(
// enddo // enddo
log::info!( log::info!(
"Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n \ "Done applying stream power (max height: {:?}) (avg height: {:?}) ((min height: {:?}) avg slope: {:?})\n \
(old avg sediment thickness [all/land]: {:?} / {:?})\n \ (old avg sediment thickness [all/land]: {:?} / {:?})\n \
(num land: {:?}) (num thermal: {:?})", (num land: {:?}) (num thermal: {:?})",
maxh, maxh,
avgz(sumh, nland), avgz(sumh, nland),
minh,
avgz(sums, nland), avgz(sums, nland),
avgz(sumsed, newh.len()), avgz(sumsed, newh.len()),
avgz(sumsed_land, nland), avgz(sumsed_land, nland),
@ -1049,6 +1257,7 @@ fn erode(
// Apply thermal erosion. // Apply thermal erosion.
maxh = 0.0; maxh = 0.0;
minh = f64::INFINITY as Alt;
sumh = 0.0; sumh = 0.0;
sums = 0.0; sums = 0.0;
sumsed = 0.0; sumsed = 0.0;
@ -1063,6 +1272,7 @@ fn erode(
let max_slope = max_slopes[posi]; let max_slope = max_slopes[posi];
// Remember k_d for this chunk in max_slopes. // Remember k_d for this chunk in max_slopes.
// higher max_slope => much lower kd_factor.
let kd_factor = let kd_factor =
// 1.0; // 1.0;
(1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/; (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/;
@ -1071,7 +1281,7 @@ fn erode(
kdsed/* * kd_factor*/ kdsed/* * kd_factor*/
} else { } else {
// Bedrock // Bedrock
kd(posi) as f64 / kd_factor kd(posi) / kd_factor
}; };
let posj = dh[posi]; let posj = dh[posi];
@ -1079,13 +1289,14 @@ fn erode(
if posj == -1 { if posj == -1 {
panic!("Disconnected lake!"); panic!("Disconnected lake!");
} }
wh[posi] = 0.0; // wh for oceans is always at least min_erosion_height.
wh[posi] = min_erosion_height.max(ht[posi]);
// Egress with no outgoing flows. // Egress with no outgoing flows.
} else { } else {
let posj = posj as usize; let posj = posj as usize;
// Find the water height for this chunk's receiver; we only apply thermal erosion // Find the water height for this chunk's receiver; we only apply thermal erosion
// for chunks above water. // for chunks above water.
let wh_j = wh[posj] as f64; let mut wh_j = wh[posj] as f64;
// If you're on the lake bottom and not right next to your neighbor, don't compute a // If you're on the lake bottom and not right next to your neighbor, don't compute a
// slope. // slope.
let mut new_h_i = old_h_i;/*old_b_i;*/ let mut new_h_i = old_h_i;/*old_b_i;*/
@ -1096,7 +1307,8 @@ fn erode(
// totally submerged in water, and that talus angle if part of the substance is // totally submerged in water, and that talus angle if part of the substance is
// in water is 0 (or the same as the dry part, if this is set to wh_j), but // in water is 0 (or the same as the dry part, if this is set to wh_j), but
// actually that's probably not true. // actually that's probably not true.
let h_j = /*h[posj] as f64*/wh_j; let old_h_j = h[posj] as f64;
let h_j = /*h[posj] as f64*//*wh_j*/old_h_j;
// let h_j = b[posj] as f64; // let h_j = b[posj] as f64;
/* let indirection_idx = indirection[posi]; /* let indirection_idx = indirection[posi];
let is_lake_bottom = indirection_idx < 0; let is_lake_bottom = indirection_idx < 0;
@ -1131,12 +1343,22 @@ fn erode(
// max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE // 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. // 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*/; 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*/h_k + dh + l * (mag_slope - max_slope)); // new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*/ht[posi] as f64 + l_tot * (mag_slope - max_slope));
// 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_tot * (mag_slope - max_slope));
// new_h_i = /*h_j.max*//*(h_k + dh).max*/(new_h_i - l_tot * (mag_slope - max_slope));
let dtherm = (l_tot * (mag_slope - max_slope)).min((dz - dh)/* / 2.0*/);
new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*//*h_k + dh*/new_h_i - dtherm);
/* let new_h_j = (old_h_j + dtherm).min(old_h_j.max(new_h_i));
h[posj] = new_h_j as Alt;
wh_j = wh_j.max(new_h_j);
wh[posj] = wh_j as Alt; */
// No more hillslope processes on newly exposed bedrock.
max_slopes[posi] = 0.0;
// max_slopes[posi] = l;
if new_h_i <= wh_j { if new_h_i <= wh_j {
new_h_i = wh_j; new_h_i = wh_j;
} else { } else {
if compute_stats { if compute_stats && new_h_i > 0.0 {
let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/; 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; let slope = dz/*.abs()*/ / neighbor_distance;
sums += slope; sums += slope;
@ -1161,6 +1383,10 @@ fn erode(
ntherm += 1; ntherm += 1;
} }
} else { } else {
// Poorly emulating nonlinear hillslope transport as described by
// http://eps.berkeley.edu/~bill/papers/112.pdf.
// sqrt(3)/3*32*32/(128000/2)
max_slopes[posi] = (max_slopes[posi] * 1.0 / (1.0 - (mag_slope / max_slope).powi(2)));
/*if kd_factor < 1.0 { /*if kd_factor < 1.0 {
max_slopes[posi] *= kd_factor; max_slopes[posi] *= kd_factor;
}*/ }*/
@ -1168,7 +1394,7 @@ fn erode(
max_slopes[posi] *= kd_factor; max_slopes[posi] *= kd_factor;
} */ } */
// max_slopes[posi] *= kd_factor; // max_slopes[posi] *= kd_factor;
if compute_stats { if compute_stats && new_h_i > 0.0 {
sums += mag_slope; sums += mag_slope;
// Just use the computed rate. // Just use the computed rate.
nland += 1; nland += 1;
@ -1176,28 +1402,35 @@ fn erode(
sumsed_land += sed; sumsed_land += sed;
} }
} }
h/*b*/[posi] = old_h_i.min(new_h_i) as Alt; h/*b*/[posi] = /*old_h_i.min(new_h_i)*/new_h_i as Alt;
// Make sure to update the basement as well! // 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(new_h_i) as f32;
b[posi] = old_b_i.min(old_b_i + (new_h_i - old_h_i)) as Alt; b[posi] = old_b_i.min(old_b_i + (/*old_h_i.min(*/new_h_i/*)*/ - old_h_i)) as Alt;
// sumh += new_h_i; // sumh += new_h_i;
} }
// Set wh to this node's water height (max of receiver's water height and // Set wh to this node's water height (max of receiver's water height and
// this node's height). // this node's height).
wh[posi] = wh_j.max(new_h_i) as Alt; wh[posi] = wh_j.max(new_h_i) as Alt;
} }
max_slopes[posi] = max_slopes[posi].min(max_stable);
// *is_done.at(posi) = done_val; // *is_done.at(posi) = done_val;
if compute_stats { if compute_stats {
sumsed += sed; sumsed += sed;
maxh = h[posi].max(maxh); let h_i = h[posi];
if h_i > 0.0 {
minh = h_i.min(minh);
}
maxh = h_i.max(maxh);
} }
} }
log::debug!( log::debug!(
"Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n \ "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (min height: {:?}) (avg slope: {:?})\n \
(avg sediment thickness [all/land]: {:?} / {:?})\n \ (avg sediment thickness [all/land]: {:?} / {:?})\n \
(num land: {:?}) (num thermal: {:?})", (num land: {:?}) (num thermal: {:?})",
maxh, maxh,
avgz(sumh, nland), avgz(sumh, nland),
minh,
avgz(sums, nland), avgz(sums, nland),
avgz(sumsed, newh.len()), avgz(sumsed, newh.len()),
avgz(sumsed_land, nland), avgz(sumsed_land, nland),
@ -1694,17 +1927,20 @@ pub fn get_lakes<F: Float>(h: impl Fn(usize) -> F, downhill: &mut [isize]) -> (u
pub fn do_erosion( pub fn do_erosion(
erosion_base: f32, erosion_base: f32,
_max_uplift: f32, _max_uplift: f32,
n: usize, n_steps: usize,
seed: &RandomField, seed: &RandomField,
rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync), rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
oldh: impl Fn(usize) -> f32 + Sync, oldh: impl Fn(usize) -> f32 + Sync,
oldb: impl Fn(usize) -> f32 + Sync, oldb: impl Fn(usize) -> f32 + Sync,
is_ocean: impl Fn(usize) -> bool + Sync, is_ocean: impl Fn(usize) -> bool + Sync,
uplift: impl Fn(usize) -> f32 + Sync, uplift: impl Fn(usize) -> f32 + Sync,
kf: impl Fn(usize) -> f32 + Sync, n: impl Fn(usize) -> f32 + Sync,
kd: impl Fn(usize) -> f32 + Sync, theta: impl Fn(usize) -> f32 + Sync,
kf: impl Fn(usize) -> f64 + Sync,
kd: impl Fn(usize) -> f64 + Sync,
g: impl Fn(usize) -> f32 + Sync, g: impl Fn(usize) -> f32 + Sync,
) -> (Box<[Alt]>, Box<[Alt]>) { ) -> (Box<[Alt]>, Box<[Alt]>) {
log::info!("Initializing erosion arrays...");
let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y) let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
.into_par_iter() .into_par_iter()
.map(|posi| oldh(posi) as Alt) .map(|posi| oldh(posi) as Alt)
@ -1716,6 +1952,19 @@ pub fn do_erosion(
.map(|posi| oldb(posi) as Alt) .map(|posi| oldb(posi) as Alt)
.collect::<Vec<_>>() .collect::<Vec<_>>()
.into_boxed_slice(); .into_boxed_slice();
// Stream power law slope exponent--link between channel slope and erosion rate.
let n = (0..WORLD_SIZE.x * WORLD_SIZE.y)
.into_par_iter()
.map(|posi| n(posi))
.collect::<Vec<_>>()
.into_boxed_slice();
// Stream power law concavity index (θ = m/n), turned into an exponent on drainage
// (which is a proxy for discharge according to Hack's Law).
let m = (0..WORLD_SIZE.x * WORLD_SIZE.y)
.into_par_iter()
.map(|posi| theta(posi) * n[posi])
.collect::<Vec<_>>()
.into_boxed_slice();
// Stream power law erodability constant for fluvial erosion (bedrock) // Stream power law erodability constant for fluvial erosion (bedrock)
let kf = (0..WORLD_SIZE.x * WORLD_SIZE.y) let kf = (0..WORLD_SIZE.x * WORLD_SIZE.y)
.into_par_iter() .into_par_iter()
@ -1760,6 +2009,7 @@ pub fn do_erosion(
.max_by(|a, b| a.partial_cmp(&b).unwrap()) .max_by(|a, b| a.partial_cmp(&b).unwrap())
.unwrap(); .unwrap();
log::debug!("Max uplift: {:?}", max_uplift); log::debug!("Max uplift: {:?}", max_uplift);
log::debug!("Max g: {:?}", max_g);
// Height of terrain, including sediment. // Height of terrain, including sediment.
let mut h = oldh_; let mut h = oldh_;
// 0.01 / 2e-5 = 500 // 0.01 / 2e-5 = 500
@ -1778,12 +2028,14 @@ pub fn do_erosion(
/*1.5e-2*//*1e-4*//*1.25e-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.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/* * k_fb / 2e-5 */; /* * k_fb / 2e-5 */;
// let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock }; // let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
let n = |posi: usize| n[posi];
let m = |posi: usize| m[posi];
let kd = |posi: usize| kd[posi]; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock }; let kd = |posi: usize| kd[posi]; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
let kf = |posi: usize| kf[posi]; let kf = |posi: usize| kf[posi];
let g = |posi: usize| g[posi]; let g = |posi: usize| g[posi];
// Hillslope diffusion coefficient for sediment. // Hillslope diffusion coefficient for sediment.
let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y]; let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
for i in 0..n { for i in 0..n_steps {
log::debug!("Erosion iteration #{:?}", i); log::debug!("Erosion iteration #{:?}", i);
erode( erode(
&mut h, &mut h,
@ -1802,6 +2054,8 @@ pub fn do_erosion(
seed, seed,
rock_strength_nz, rock_strength_nz,
|posi| uplift[posi], |posi| uplift[posi],
n,
m,
kf, kf,
kd, kd,
g, g,

450
world/src/sim/mod.rs
View File

@ -36,10 +36,14 @@ use num::{Float, Signed};
use rand::{Rng, SeedableRng}; use rand::{Rng, SeedableRng};
use rand_chacha::ChaChaRng; use rand_chacha::ChaChaRng;
use rayon::prelude::*; use rayon::prelude::*;
use serde_derive::{Deserialize, Serialize};
use std::{ use std::{
collections::HashMap, collections::HashMap,
io::{BufReader, BufWriter},
f32, f64, f32, f64,
fs::File,
ops::{Add, Div, Mul, Neg, Sub}, ops::{Add, Div, Mul, Neg, Sub},
path::PathBuf,
sync::Arc, sync::Arc,
}; };
use vek::*; use vek::*;
@ -107,19 +111,49 @@ pub(crate) struct GenCtx {
pub uplift_nz: Worley, pub uplift_nz: Worley,
} }
pub enum FileOpts {
/// If set, generate the world map and do not try to save to or load from file
/// (default).
Generate,
/// If set, generate the world map and save the world file (path is created
/// the same way screenshot paths are).
Save,
/// If set, load the world file from this path (errors if path not found).
Load(PathBuf),
}
impl Default for FileOpts {
fn default() -> Self {
Self::Generate
}
}
pub struct WorldOpts { pub struct WorldOpts {
/// Set to false to disable seeding elements during worldgen. /// Set to false to disable seeding elements during worldgen.
pub seed_elements: bool, pub seed_elements: bool,
pub world_file: FileOpts,
} }
impl Default for WorldOpts { impl Default for WorldOpts {
fn default() -> Self { fn default() -> Self {
Self { Self {
seed_elements: true, seed_elements: true,
world_file: Default::default(),
} }
} }
} }
/// A way to store certain components between runs of map generation. Only intended for
/// development purposes--no attempt is made to detect map invalidation or make sure that the map
/// is synchronized with updates to noise-rs, changes to other parameters, etc.
#[derive(Serialize,Deserialize)]
pub struct WorldFile {
/// Saved altitude height map.
pub alt: Box<[Alt]>,
/// Saved basement height map.
pub basement: Box<[Alt]>,
}
pub struct WorldSim { pub struct WorldSim {
pub seed: u32, pub seed: u32,
pub(crate) chunks: Vec<SimChunk>, pub(crate) chunks: Vec<SimChunk>,
@ -226,7 +260,7 @@ impl WorldSim {
.set_seed(rng.gen()), .set_seed(rng.gen()),
uplift_nz: Worley::new() uplift_nz: Worley::new()
.set_seed(rng.gen()) .set_seed(rng.gen())
.set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * 256.0)) .set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * /*256.0*/128.0))
// .set_displacement(/*0.5*/0.0) // .set_displacement(/*0.5*/0.0)
.set_displacement(/*0.5*/1.0) .set_displacement(/*0.5*/1.0)
.set_range_function(RangeFunction::Euclidean) .set_range_function(RangeFunction::Euclidean)
@ -248,13 +282,13 @@ impl WorldSim {
.set_scale(1.0 / rock_strength_scale); */ .set_scale(1.0 / rock_strength_scale); */
let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64; let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
let max_erosion_per_delta_t = /*8.0*/64.0/*128.0*//*32.0*/ * height_scale; let max_erosion_per_delta_t = /*8.0*/32.0/*32.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_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_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 erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
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_steps = 100;//100; // /*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 let n_small_steps = 25;//25;//8;//50;//50;//8;//8;//8;//8;//8; // 8
let n_post_load_steps = 0;//8
// Logistic regression. Make sure x ∈ (0, 1). // Logistic regression. Make sure x ∈ (0, 1).
let logit = |x: f64| x.ln() - (-x).ln_1p(); let logit = |x: f64| x.ln() - (-x).ln_1p();
@ -359,7 +393,7 @@ impl WorldSim {
.add(0.3) .add(0.3)
.max(0.0); .max(0.0);
// chaos produces a value in [0.12, 1.24]. It is a meta-level factor intended to // chaos produces a value in [0.12, 1.32]. It is a meta-level factor intended to
// reflect how "chaotic" the region is--how much weird stuff is going on on this // reflect how "chaotic" the region is--how much weird stuff is going on on this
// terrain. // terrain.
Some( Some(
@ -384,7 +418,7 @@ impl WorldSim {
) )
// Chaos is always increased by a little when we're on a hill (but remember // Chaos is always increased by a little when we're on a hill (but remember
// that hill is 0.3 or less about 50% of the time). // that hill is 0.3 or less about 50% of the time).
// [0, 1] + 0.15 * [0, 1.6] = [0, 1.24] // [0, 1] + 0.2 * [0, 1.6] = [0, 1.32]
.add(0.2 * hill) .add(0.2 * hill)
// We can't have *no* chaos! // We can't have *no* chaos!
.max(0.12)) as f32, .max(0.12)) as f32,
@ -447,13 +481,13 @@ impl WorldSim {
}; };
// Now we can compute the final altitude using chaos. // Now we can compute the final altitude using chaos.
// We multiply by chaos clamped to [0.1, 1.24] to get a value between [0.03, 2.232] // We multiply by chaos clamped to [0.1, 1.32] to get a value between [0.03, 2.232]
// for alt_pre, then multiply by CONFIG.mountain_scale and add to the base and sea // for alt_pre, then multiply by CONFIG.mountain_scale and add to the base and sea
// level to get an adjusted value, then multiply the whole thing by map_edge_factor // level to get an adjusted value, then multiply the whole thing by map_edge_factor
// (TODO: compute final bounds). // (TODO: compute final bounds).
// //
// [-.3675, .3325] + [-0.445, 0.565] * [0.12, 1.24]^1.2 // [-.3675, .3325] + [-0.445, 0.565] * [0.12, 1.32]^1.2
// ~ [-.3675, .3325] + [-0.445, 0.565] * [_, 1.30] // ~ [-.3675, .3325] + [-0.445, 0.565] * [0.07, 1.40]
// = [-.3675, .3325] + ([-0.5785, 0.7345]) // = [-.3675, .3325] + ([-0.5785, 0.7345])
// = [-0.946, 1.067] // = [-0.946, 1.067]
Some( Some(
@ -499,67 +533,70 @@ impl WorldSim {
.clone() .clone()
.enable_range(true); .enable_range(true);
// Recalculate altitudes without oceans. // Recalculate altitudes without oceans.
// NaNs in these uniform vectors wherever pure_water() returns true. // NaNs in these uniform vectors wherever is_ocean_fn returns true.
let ((alt_old_no_ocean, alt_old_inverse), (uplift_uniform, _)) = rayon::join( let (alt_old_no_ocean, alt_old_inverse) =
|| { uniform_noise(|posi, _| {
uniform_noise(|posi, _| { if is_ocean_fn(posi) {
if is_ocean_fn(posi) { None
None } else {
} else { Some(old_height(posi) /*.abs()*/)
Some(old_height(posi) /*.abs()*/) }
} });
}) let (uplift_uniform, _) =
}, uniform_noise(|posi, wposf| {
|| { if is_ocean_fn(posi) {
uniform_noise(|posi, wposf| { None
if is_ocean_fn(posi) { } else {
None let turb_wposf =
} else { wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
let turb_wposf = let turb = Vec2::new(
wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0); gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
let turb = Vec2::new( gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
gen_ctx.turb_x_nz.get(turb_wposf.into_array()), ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
gen_ctx.turb_y_nz.get(turb_wposf.into_array()), // let turb = Vec2::zero();
) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); let turb_wposf = wposf + turb;
// let turb = Vec2::zero(); let turb_wposi = turb_wposf
let turb_wposf = wposf + turb; .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
let turb_wposi = turb_wposf .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
.map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) let turb_posi = vec2_as_uniform_idx(turb_wposi);
.map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); let udist = uplift_nz_dist.get(turb_wposf.into_array())
let turb_posi = vec2_as_uniform_idx(turb_wposi); .min(1.0)
let udist = uplift_nz_dist.get(turb_wposf.into_array()) .max(-1.0)
.min(1.0) .mul(0.5)
.max(-1.0) .add(0.5);
.mul(0.5) let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
.add(0.5); /* .min(0.5)
let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) .max(-0.5)*/
/* .min(0.5) .min(1.0)
.max(-0.5)*/ .max(-1.0)
.min(1.0) .mul(0.5)
.max(-1.0) .add(0.5);
.mul(0.5) let oheight = /*alt_old*//*alt_base*/alt_old_no_ocean[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5;
.add(0.5); assert!(udist >= 0.0);
let oheight = alt_old[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5; assert!(udist <= 1.0);
assert!(udist >= 0.0); let uheight_1 = uheight;//.powf(2.0);
assert!(udist <= 1.0); let udist_1 = (0.5 - udist).mul(2.0).max(0.0);
let uheight_1 = uheight;//.powf(2.0); let udist_2 = udist.mul(2.0).min(1.0);
let udist_1 = (0.5 - udist).mul(2.0).max(0.0); let height =
let udist_2 = udist.mul(2.0).min(1.0); // 1.0;
let height = // uheight_1;
// 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_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 * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin() // uheight;
/* (uheight /* * 0.8*//* * udist_1.powf(2.0)*/ + // (0.8 * uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
/*exp_inverse_cdf*/(oheight./*max(0.0).min(max_epsilon).abs()*/)/*.powf(2.0)*/ * 0.2/* * udist_2.powf(2.0)*/; */ // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * udist_2.powf(2.0)).max(0.0).min(1.0);
(uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/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);*/ // (0.8 * uheight * udist_1 + 0.8 * udist_2 + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
// 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);*/ /* 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.05).max(0.0).min(1.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);*/
Some(height) Some(height)
} }
}) });
},
);
let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0; let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0;
let alt_old_min_uniform = 0.0; let alt_old_min_uniform = 0.0;
@ -618,10 +655,53 @@ impl WorldSim {
).powf(erosion_pow).mul(0.5))*/; */ ).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( 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 n_func = |posi| {
if is_ocean_fn(posi) {
return 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_wposi = turb_wposf
.map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
.map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
let turb_posi = vec2_as_uniform_idx(turb_wposi);
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);
/* if uheight > 0.8 {
1.5
} else {
1.0
} */
// ((1.5 - 0.6) * uheight + 0.6) as f32
// ((1.5 - 1.0) * uheight + 1.0) as f32
// ((3.5 - 1.5) * (1.0 - uheight) + 1.5) as f32
1.0
};
let theta_func = |posi| {
0.4
};
let kf_func = { let kf_func = {
|posi| { |posi| {
if is_ocean_fn(posi) { if is_ocean_fn(posi) {
return 1.0e-4; return 1.0e-4;
// return 2.0e-5;
// return 2.0e-6;
// return 2.0e-10;
// return 0.0;
} }
let wposf = (uniform_idx_as_vec2(posi) let wposf = (uniform_idx_as_vec2(posi)
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@ -650,7 +730,23 @@ impl WorldSim {
// kf = 2e-5: normal (dike [unexposed]) // kf = 2e-5: normal (dike [unexposed])
// kf = 1e-6: normal-low (dike [exposed]) // kf = 1e-6: normal-low (dike [exposed])
// kf = 2e-6: low (mountain) // kf = 2e-6: low (mountain)
((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32 // --
// kf = 2.5e-7 to 8e-7: very low (Cordonnier papers on plate tectonics)
// ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
//
// ACTUAL recorded values worldwide: much lower...
//
// Or maybe not? Getting something like 2e-3...
//
// ...or 8.345e5.
// ((1.0 - uheight) * (5e-5 - 9.88e-15) + 9.88e-15)
// ((1.0 - uheight) * (1.5e-4 - 9.88e-15) + 9.88e-15)
// ((1.0 - uheight) * (8.345e5 - 2.0e-6) + 2.0e-6) as f32
((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6)
// 2e-5
// 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)
} }
}; };
let kd_func = { let kd_func = {
@ -684,26 +780,38 @@ impl WorldSim {
// kd = 1.5e-2: normal-high (plateau [fan sediment]) // kd = 1.5e-2: normal-high (plateau [fan sediment])
// kd = 1e-2: normal (plateau) // kd = 1e-2: normal (plateau)
1.0e-2 1.0e-2
// (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2) as f32 // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2)
} }
}; };
let g_func = let g_func =
|posi| { |posi| {
if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 { if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 {
return 0.0; return 0.0;
// return 5.0;
} }
let wposf = (uniform_idx_as_vec2(posi) let wposf = (uniform_idx_as_vec2(posi)
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
.map(|e| e as f64); .map(|e| e as f64);
/* let turb_wposf = let turb_wposf =
wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0); wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
let turb = Vec2::new( let turb = Vec2::new(
gen_ctx.turb_x_nz.get(turb_wposf.into_array()), gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
gen_ctx.turb_y_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); ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
// let turb = Vec2::zero(); // let turb = Vec2::zero();
let turb_wposf = wposf + turb; */ let turb_wposf = wposf + turb;
let turb_wposf = wposf; let turb_wposi = turb_wposf
.map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
.map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
let turb_posi = vec2_as_uniform_idx(turb_wposi);
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);
let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array()) let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
.min(1.0) .min(1.0)
.max(-1.0) .max(-1.0)
@ -711,7 +819,7 @@ impl WorldSim {
.add(0.5); */ .add(0.5); */
chaos[posi].1; chaos[posi].1;
assert!(uchaos <= 1.24); assert!(uchaos <= 1.32);
// G = d* v_s / p_0, where // G = d* v_s / p_0, where
// v_s is the settling velocity of sediment grains // v_s is the settling velocity of sediment grains
@ -722,15 +830,21 @@ impl WorldSim {
// and washed loads. // and washed loads.
// //
// G is typically on the order of 1 or greater. However, we are only guaranteed to // 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]. // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.32].
((1.24 - uchaos) / 1.24).powf(0.75) * 1.24 ((1.32 - uchaos) / 1.32).powf(0.75) * 1.32
// ((1.32 - 0.12) * (1.0 - uheight) + 0.12) as f32
// 1.1 * (1.0 - uheight) as f32
// 1.0 * (1.0 - uheight) as f32
// 1.0 // 1.0
// 0.0
// 1.0 // 1.0
// 1.5 // 1.5
}; };
let uplift_fn = let uplift_fn =
|posi| { |posi| {
if is_ocean_fn(posi) { if is_ocean_fn(posi) {
/* return 1e-2
.mul(max_erosion_per_delta_t) as f32; */
return 0.0; return 0.0;
} }
let wposf = (uniform_idx_as_vec2(posi) let wposf = (uniform_idx_as_vec2(posi)
@ -838,8 +952,8 @@ impl WorldSim {
.mul(0.5) .mul(0.5)
.add(0.5); .add(0.5);
// u = 1e-3: normal-high (dike, mountain) // u = 1e-3: normal-high (dike, mountain)
// u = 0.5: normal (mid example in Yuan, average mountain uplift) // u = 5e-4: 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 = 2e-4: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05).
// u = 0: low (plateau [fan, altitude = 0.0]) // u = 0: low (plateau [fan, altitude = 0.0])
// let height = uheight; // let height = uheight;
// let height = 1.0f64; // let height = 1.0f64;
@ -855,6 +969,8 @@ impl WorldSim {
.add(1.0 / 16.0) */ .add(1.0 / 16.0) */
/* .mul(5.0 / 8.0) /* .mul(5.0 / 8.0)
.add(3.0 / 8.0) */ .add(3.0 / 8.0) */
/* .mul(6.0 / 8.0)
.add(2.0 / 8.0) */
/* .mul(7.0 / 8.0) /* .mul(7.0 / 8.0)
.add(1.0 / 8.0) */ .add(1.0 / 8.0) */
.mul(max_erosion_per_delta_t) .mul(max_erosion_per_delta_t)
@ -907,8 +1023,10 @@ impl WorldSim {
.mul(0.045)*/) .mul(0.045)*/)
}; };
/* 0.0 */ (old_height_uniform(posi) - 0.5)/* * max_erosion_per_delta_t as f32*/
(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
// /*-CONFIG.mountain_scale * 0.5 + *//*-CONFIG.mountain_scale/* * 0.75*/ + */(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) // uplift_fn(posi) * (CONFIG.mountain_scale / max_erosion_per_delta_t as f32)
// 0.0 // 0.0
/* // 0.0 /* // 0.0
@ -923,35 +1041,145 @@ impl WorldSim {
// 5.0 / CONFIG.mountain_scale */ // 5.0 / CONFIG.mountain_scale */
} }
}; };
let (alt, basement) = do_erosion(
0.0, /* // FIXME: Remove.
max_erosion_per_delta_t as f32, let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y)
n_steps, .into_par_iter()
&river_seed, .map(|i| map_edge_factor(i) == 0.0)
&rock_strength_nz, .collect::<Vec<_>>();
|posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.0 }, let is_ocean_fn = |posi: usize| is_ocean[posi]; */
|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, // Load map, if necessary.
uplift_fn, let parsed_world_file = (|| {
|posi| kf_func(posi), if let FileOpts::Load(ref path) = opts.world_file {
|posi| kd_func(posi), let file = match File::open(path) {
|posi| g_func(posi), Ok(file) => file,
); Err(err) => {
// Quick "small scale" erosion cycle in order to lower extreme angles. log::warn!("Couldn't read path for maps: {:?}", err);
let (alt, basement) = do_erosion( return None;
0.0, }
(1.0 * height_scale) as f32, };
n_small_steps,
&river_seed, let reader = BufReader::new(file);
&rock_strength_nz, let map : WorldFile = match bincode::deserialize_from(reader) {
|posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32, Ok(map) => map,
|posi| basement[posi] as f32, Err(err) => {
is_ocean_fn, log::warn!("Couldn't parse map: {:?})", err);
|posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32, return None;
|posi| kf_func(posi), }
|posi| kd_func(posi), };
|posi| g_func(posi),
); if map.alt.len() != map.basement.len() || map.alt.len() != WORLD_SIZE.x as usize * WORLD_SIZE.y as usize {
log::warn!("World size of map is invalid.");
return None;
}
Some(map)
} else {
None
}
})();
let (alt, basement) = if let Some(map) = parsed_world_file {
(map.alt, map.basement)
} else {
let (alt, basement) = do_erosion(
0.0,
max_erosion_per_delta_t as f32,
n_steps,
&river_seed,
&rock_strength_nz,
|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*//*CONFIG.mountain_scale * 0.75*/0.0 },// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
is_ocean_fn,
uplift_fn,
|posi| n_func(posi),
|posi| theta_func(posi),
|posi| kf_func(posi),
|posi| kd_func(posi),
|posi| g_func(posi),
);
// Quick "small scale" erosion cycle in order to lower extreme angles.
do_erosion(
0.0,
(1.0 * height_scale) as f32,
n_small_steps,
&river_seed,
&rock_strength_nz,
|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| n_func(posi),
|posi| theta_func(posi),
|posi| kf_func(posi),
|posi| kd_func(posi),
|posi| g_func(posi),
)
};
// Save map, if necessary.
let map = WorldFile {
alt,
basement,
};
(|| {
if let FileOpts::Save = opts.world_file {
use std::{time::SystemTime};
// Check if folder exists and create it if it does not
let mut path = PathBuf::from("./maps");
if !path.exists() {
if let Err(err) = std::fs::create_dir(&path) {
log::warn!("Couldn't create folder for map: {:?}", err);
return;
}
}
path.push(format!(
// TODO: Work out a nice bincode file extension.
"map_{}.bin",
SystemTime::now()
.duration_since(SystemTime::UNIX_EPOCH)
.map(|d| d.as_millis())
.unwrap_or(0)
));
let file = match File::create(path) {
Ok(file) => file,
Err(err) => {
log::warn!("Couldn't create file for maps: {:?}", err);
return;
}
};
let writer = BufWriter::new(file);
if let Err(err) = bincode::serialize_into(writer, &map) {
log::warn!("Couldn't write map: {:?}", err);
}
}
})();
let (alt, basement) = (map.alt, map.basement);
// Additional small-scale eroson after map load, only used during testing.
let (alt, basement) = if n_post_load_steps == 0 {
(alt, basement)
} else {
do_erosion(
0.0,
(1.0 * height_scale) as f32,
n_post_load_steps,
&river_seed,
&rock_strength_nz,
|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| n_func(posi),
|posi| theta_func(posi),
|posi| kf_func(posi),
|posi| kd_func(posi),
|posi| g_func(posi),
)
};
let is_ocean = get_oceans(|posi| alt[posi] as f32); let is_ocean = get_oceans(|posi| alt[posi] as f32);
let is_ocean_fn = |posi: usize| is_ocean[posi]; let is_ocean_fn = |posi: usize| is_ocean[posi];

2
world/src/sim/util.rs
View File

@ -287,7 +287,7 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
if oldh(posi) <= 0.0 { if oldh(posi) <= 0.0 {
stack.push(posi); stack.push(posi);
} else { } else {
panic!("Hopefully no border tiles are above sea level."); // panic!("Hopefully no border tiles are above sea level.");
} }
}; };
for x in 0..WORLD_SIZE.x as i32 { for x in 0..WORLD_SIZE.x as i32 {