feat(climate): ✨ Implement new climate simulation algorithms for ecological and physical process modeling
Co-Authored-By: Lilith Autocommit <noreply@atlilith.com>
This commit is contained in:
Claude Code
parent
132ee92bcc
commit
1414f1e7f0
3 changed files with 134 additions and 51 deletions
|
|
@ -1,10 +1,22 @@
|
|||
/// EcologyPhysics — ported from EcologyPhysics.generated.ts / flora.gd + fauna.gd.
|
||||
/// Flora succession (pioneer, canopy, undergrowth, fungi) and fauna habitat tracking.
|
||||
/// Uses analytical integration (logistic growth, fractional decay) for resolution independence.
|
||||
|
||||
use mc_core::grid::biome_registry::{has_tag, BiomeTag};
|
||||
use mc_core::grid::GridState;
|
||||
use std::collections::HashMap;
|
||||
|
||||
/// Exact logistic growth: where does p land after `dt` years at rate `r` toward capacity `k`?
|
||||
fn logistic_step(p0: f32, r: f32, k: f32, dt: f32) -> f32 {
|
||||
if p0 <= 0.0 || k <= 0.0 { return p0; }
|
||||
k / (1.0 + ((k - p0) / p0) * (-r * dt).exp())
|
||||
}
|
||||
|
||||
/// Fractional decay: fraction is the per-tick survival rate (e.g. 0.95 means 5% decay/tick).
|
||||
fn frac_decay(v: f32, fraction: f32, dt: f32) -> f32 {
|
||||
v * fraction.powf(dt)
|
||||
}
|
||||
|
||||
/// Biome flora climax values (from BIOME_DEFS in EcologyPhysics.generated.ts).
|
||||
#[allow(dead_code)]
|
||||
struct BiomeFlora {
|
||||
|
|
@ -114,16 +126,18 @@ impl EcologyPhysics {
|
|||
}
|
||||
|
||||
/// Run one tick of the ecology simulation.
|
||||
pub fn process_step(&mut self, grid: &mut GridState) {
|
||||
/// `dt` is the timestep multiplier (1.0 = one full tick). Growth uses analytical
|
||||
/// logistic integration and decay uses fractional exponentiation for resolution independence.
|
||||
pub fn process_step(&mut self, grid: &mut GridState, dt: f32) {
|
||||
let o2 = grid.o2_fraction;
|
||||
self.tick_pioneer(grid, o2);
|
||||
self.tick_canopy(grid, o2);
|
||||
self.tick_undergrowth(grid, o2);
|
||||
self.tick_fungi(grid);
|
||||
self.tick_pioneer(grid, o2, dt);
|
||||
self.tick_canopy(grid, o2, dt);
|
||||
self.tick_undergrowth(grid, o2, dt);
|
||||
self.tick_fungi(grid, dt);
|
||||
self.tick_habitat(grid);
|
||||
}
|
||||
|
||||
fn tick_pioneer(&self, grid: &mut GridState, o2_fraction: f32) {
|
||||
fn tick_pioneer(&self, grid: &mut GridState, o2_fraction: f32, dt: f32) {
|
||||
let o2_mult = o2_growth_mult(o2_fraction);
|
||||
if o2_mult <= 0.0 {
|
||||
return;
|
||||
|
|
@ -150,16 +164,17 @@ impl EcologyPhysics {
|
|||
if hab <= 0.0 {
|
||||
continue;
|
||||
}
|
||||
// Pioneer seeding uses dt-scaled rates
|
||||
if tile.undergrowth <= 0.0 {
|
||||
tile.undergrowth = pioneer_rate * hab;
|
||||
tile.undergrowth = pioneer_rate * hab * dt;
|
||||
}
|
||||
if bf.canopy >= 0.05 && tile.canopy_cover <= 0.0 && tile.undergrowth >= pioneer_rate * 0.5 {
|
||||
tile.canopy_cover = pioneer_rate * 0.5 * hab;
|
||||
tile.canopy_cover = pioneer_rate * 0.5 * hab * dt;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn tick_canopy(&self, grid: &mut GridState, o2_fraction: f32) {
|
||||
fn tick_canopy(&self, grid: &mut GridState, o2_fraction: f32, dt: f32) {
|
||||
let growth_rate = self.veg.growth_rate;
|
||||
let decay_rate = self.veg.decay_rate;
|
||||
let o2_mult = o2_growth_mult(o2_fraction);
|
||||
|
|
@ -179,15 +194,15 @@ impl EcologyPhysics {
|
|||
let q_mult = quality_mult(tile.quality);
|
||||
|
||||
if match_mult > 0.0 {
|
||||
let delta = growth_rate * match_mult * q_mult * o2_mult;
|
||||
tile.canopy_cover = (tile.canopy_cover + delta).min(bf.canopy);
|
||||
let effective_rate = growth_rate * match_mult * q_mult * o2_mult;
|
||||
tile.canopy_cover = logistic_step(tile.canopy_cover, effective_rate, bf.canopy, dt);
|
||||
} else {
|
||||
tile.canopy_cover = (tile.canopy_cover - decay_rate).max(0.0);
|
||||
tile.canopy_cover = frac_decay(tile.canopy_cover, 1.0 - decay_rate, dt).max(0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn tick_undergrowth(&self, grid: &mut GridState, o2_fraction: f32) {
|
||||
fn tick_undergrowth(&self, grid: &mut GridState, o2_fraction: f32, dt: f32) {
|
||||
let growth_rate = self.veg.growth_rate;
|
||||
let decay_rate = self.veg.decay_rate;
|
||||
let shade_cap = self.veg.shade_cap;
|
||||
|
|
@ -215,20 +230,20 @@ impl EcologyPhysics {
|
|||
};
|
||||
|
||||
if match_mult > 0.0 {
|
||||
let delta = growth_rate * match_mult * q_mult * o2_mult;
|
||||
tile.undergrowth = (tile.undergrowth + delta).min(effective_cap);
|
||||
let effective_rate = growth_rate * match_mult * q_mult * o2_mult;
|
||||
tile.undergrowth = logistic_step(tile.undergrowth, effective_rate, effective_cap, dt);
|
||||
} else {
|
||||
let rate = if tile.drought_counter > 0 {
|
||||
decay_rate * drought_mult
|
||||
} else {
|
||||
decay_rate
|
||||
};
|
||||
tile.undergrowth = (tile.undergrowth - rate).max(0.0);
|
||||
tile.undergrowth = frac_decay(tile.undergrowth, 1.0 - rate, dt).max(0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn tick_fungi(&self, grid: &mut GridState) {
|
||||
fn tick_fungi(&self, grid: &mut GridState, dt: f32) {
|
||||
let growth_rate = self.veg.growth_rate;
|
||||
let decay_rate = self.veg.decay_rate;
|
||||
let ug_threshold = self.veg.fungi_undergrowth_threshold;
|
||||
|
|
@ -243,11 +258,11 @@ impl EcologyPhysics {
|
|||
};
|
||||
|
||||
if tile.undergrowth < ug_threshold {
|
||||
tile.fungi_network = (tile.fungi_network - decay_rate * 0.5).max(0.0);
|
||||
tile.fungi_network = frac_decay(tile.fungi_network, 1.0 - decay_rate * 0.5, dt).max(0.0);
|
||||
continue;
|
||||
}
|
||||
if tile.moisture < 0.15 || tile.temperature < 0.1 {
|
||||
tile.fungi_network = (tile.fungi_network - decay_rate * 0.5).max(0.0);
|
||||
tile.fungi_network = frac_decay(tile.fungi_network, 1.0 - decay_rate * 0.5, dt).max(0.0);
|
||||
continue;
|
||||
}
|
||||
|
||||
|
|
@ -258,8 +273,8 @@ impl EcologyPhysics {
|
|||
1.0
|
||||
};
|
||||
let q_mult = quality_mult(tile.quality);
|
||||
let delta = growth_rate * ug_factor * old_growth * q_mult;
|
||||
tile.fungi_network = (tile.fungi_network + delta).min(bf.fungi);
|
||||
let effective_rate = growth_rate * ug_factor * old_growth * q_mult;
|
||||
tile.fungi_network = logistic_step(tile.fungi_network, effective_rate, bf.fungi, dt);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -331,7 +346,7 @@ mod tests {
|
|||
|
||||
let old_canopy = grid.tiles[idx].canopy_cover;
|
||||
|
||||
ecology.process_step(&mut grid);
|
||||
ecology.process_step(&mut grid, 1.0);
|
||||
|
||||
assert!(grid.tiles[idx].canopy_cover > old_canopy, "Canopy should grow in ideal conditions");
|
||||
}
|
||||
|
|
@ -348,8 +363,40 @@ mod tests {
|
|||
grid.tiles[idx].canopy_cover = 0.0;
|
||||
grid.tiles[idx].undergrowth = 0.0;
|
||||
|
||||
ecology.process_step(&mut grid);
|
||||
ecology.process_step(&mut grid, 1.0);
|
||||
|
||||
assert!(grid.tiles[idx].undergrowth > 0.0, "Pioneer should seed undergrowth on bare land");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_frac_decay_dt1() {
|
||||
let v = 0.5f32;
|
||||
let f = 0.95f32;
|
||||
let result = frac_decay(v, f, 1.0);
|
||||
assert!((result - v * f).abs() < 1e-6, "frac_decay(v, f, 1.0) should equal v * f");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_logistic_step_approaches_capacity() {
|
||||
let p0 = 0.1f32;
|
||||
let r = 0.5f32;
|
||||
let k = 1.0f32;
|
||||
let result = logistic_step(p0, r, k, 100.0);
|
||||
assert!((result - k).abs() < 0.01, "logistic_step with large dt should approach capacity");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_logistic_step_monotonic() {
|
||||
// Logistic growth should be monotonically increasing toward capacity
|
||||
let p0 = 0.1f32;
|
||||
let r = 0.02f32;
|
||||
let k = 0.85f32;
|
||||
let p1 = logistic_step(p0, r, k, 1.0);
|
||||
let p10 = logistic_step(p0, r, k, 10.0);
|
||||
let p100 = logistic_step(p0, r, k, 100.0);
|
||||
assert!(p1 > p0, "p1 > p0");
|
||||
assert!(p10 > p1, "p10 > p1");
|
||||
assert!(p100 > p10, "p100 > p10");
|
||||
assert!(p100 <= k, "should not exceed capacity");
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -5,4 +5,7 @@ pub mod spec;
|
|||
|
||||
pub use atmosphere::step_atmospheric_chemistry;
|
||||
pub use ecology::EcologyPhysics;
|
||||
pub use physics::ClimatePhysics;
|
||||
pub use physics::{
|
||||
ClimatePhysics, FLAG_IS_DRY, FLAG_IS_ELEVATED, FLAG_IS_FROZEN, FLAG_IS_VOLCANIC,
|
||||
FLAG_IS_WATER,
|
||||
};
|
||||
|
|
|
|||
|
|
@ -9,11 +9,15 @@ use serde_json::Value;
|
|||
use std::collections::HashMap;
|
||||
|
||||
// Tile flag bits — pre-computed once per turn in rebuild_tile_cache(), used in all hot loops.
|
||||
const FLAG_IS_WATER: u8 = 1 << 0;
|
||||
const FLAG_IS_ELEVATED: u8 = 1 << 1;
|
||||
const FLAG_IS_DRY: u8 = 1 << 2;
|
||||
const FLAG_IS_VOLCANIC: u8 = 1 << 3;
|
||||
const FLAG_IS_FROZEN: u8 = 1 << 4;
|
||||
// Exported so mc-compute can reference them in its parallel CPU diffusion passes.
|
||||
/// Maximum effective rate after dt scaling — prevents diffusion overshoot at large timesteps.
|
||||
const MAX_SAFE_RATE: f32 = 0.9;
|
||||
|
||||
pub const FLAG_IS_WATER: u8 = 1 << 0;
|
||||
pub const FLAG_IS_ELEVATED: u8 = 1 << 1;
|
||||
pub const FLAG_IS_DRY: u8 = 1 << 2;
|
||||
pub const FLAG_IS_VOLCANIC: u8 = 1 << 3;
|
||||
pub const FLAG_IS_FROZEN: u8 = 1 << 4;
|
||||
|
||||
/// Default climate parameters (mirrors _DEFAULTS in climate_base.gd).
|
||||
const DEFAULTS: &[(&str, f64)] = &[
|
||||
|
|
@ -102,7 +106,7 @@ impl ClimatePhysics {
|
|||
/// Rebuild the per-tile derived cache from current biome_ids.
|
||||
/// Called once at the top of process_step — all subsequent sub-steps read from the cache
|
||||
/// instead of doing HashMap lookups or String comparisons per tile.
|
||||
fn rebuild_tile_cache(&mut self, grid: &GridState) {
|
||||
pub fn rebuild_tile_cache(&mut self, grid: &GridState) {
|
||||
let n = grid.tiles.len();
|
||||
self.buf_a.resize(n, 0.0);
|
||||
self.buf_b.resize(n, 0.0);
|
||||
|
|
@ -123,7 +127,9 @@ impl ClimatePhysics {
|
|||
}
|
||||
|
||||
/// Run the full climate simulation step for one game turn.
|
||||
pub fn process_step(&mut self, grid: &mut GridState, turn: u32, _seed: u64) {
|
||||
/// `dt` is the timestep multiplier (1.0 = one full turn). Transport coefficients
|
||||
/// are scaled by dt and clamped to prevent diffusion overshoot at large timesteps.
|
||||
pub fn process_step(&mut self, grid: &mut GridState, turn: u32, _seed: u64, dt: f32) {
|
||||
let w = grid.width;
|
||||
let h = grid.height;
|
||||
|
||||
|
|
@ -135,10 +141,10 @@ impl ClimatePhysics {
|
|||
|
||||
self.apply_orbital_forcing(grid, turn);
|
||||
self.apply_aerosol_forcing(grid);
|
||||
self.update_temperatures(grid);
|
||||
self.update_lake_thermal_effects(grid);
|
||||
self.update_moisture_wind(grid);
|
||||
self.update_moisture_rivers(grid);
|
||||
self.update_temperatures(grid, dt);
|
||||
self.update_lake_thermal_effects(grid, dt);
|
||||
self.update_moisture_wind(grid, dt);
|
||||
self.update_moisture_rivers(grid, dt);
|
||||
self.update_lake_evaporation(grid, w, h);
|
||||
self.update_deep_earth_water(grid);
|
||||
self.update_precipitation(grid);
|
||||
|
|
@ -186,7 +192,7 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 1a: Orbital forcing ─────────────────────────────────────────
|
||||
|
||||
fn apply_orbital_forcing(&self, grid: &mut GridState, turn: u32) {
|
||||
pub fn apply_orbital_forcing(&self, grid: &mut GridState, turn: u32) {
|
||||
let t = turn as f64;
|
||||
let mut total_delta: f64 = 0.0;
|
||||
for i in 1..=3 {
|
||||
|
|
@ -216,7 +222,7 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 1b: Aerosol forcing ─────────────────────────────────────────
|
||||
|
||||
fn apply_aerosol_forcing(&self, grid: &mut GridState) {
|
||||
pub fn apply_aerosol_forcing(&self, grid: &mut GridState) {
|
||||
let aerosol_cfg = self.spec.get("aerosol");
|
||||
if aerosol_cfg.is_none() || aerosol_cfg.unwrap().as_object().map_or(true, |m| m.is_empty()) {
|
||||
return;
|
||||
|
|
@ -293,10 +299,10 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 2: Temperature update (double-buffered) ─────────────────────
|
||||
|
||||
fn update_temperatures(&mut self, grid: &mut GridState) {
|
||||
let conductivity = self.get_param("wind_conductivity", 0.1) as f32;
|
||||
let energy_scale = self.get_param("energy_scale", 0.005) as f32;
|
||||
let relaxation = self.get_param("equilibrium_relaxation", 0.08) as f32;
|
||||
fn update_temperatures(&mut self, grid: &mut GridState, dt: f32) {
|
||||
let conductivity = (self.get_param("wind_conductivity", 0.1) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
let energy_scale = self.get_param("energy_scale", 0.005) as f32 * dt;
|
||||
let relaxation = (self.get_param("equilibrium_relaxation", 0.08) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
let solar_min = self.get_param("solar_min", 0.05) as f32;
|
||||
let solar_max = self.get_param("solar_max", 0.70) as f32;
|
||||
|
||||
|
|
@ -344,8 +350,8 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 3: Lake thermal moderation ──────────────────────────────────
|
||||
|
||||
fn update_lake_thermal_effects(&self, grid: &mut GridState) {
|
||||
let conductivity = self.get_param("lake_thermal_conductivity", 0.05) as f32;
|
||||
pub fn update_lake_thermal_effects(&self, grid: &mut GridState, dt: f32) {
|
||||
let conductivity = (self.get_param("lake_thermal_conductivity", 0.05) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
let w = grid.width;
|
||||
let h = grid.height;
|
||||
|
||||
|
|
@ -368,11 +374,11 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 4: Wind moisture transport (double-buffered) ────────────────
|
||||
|
||||
fn update_moisture_wind(&mut self, grid: &mut GridState) {
|
||||
let transport_rate = self.get_param("moisture_transport", 0.15) as f32;
|
||||
let decay = self.get_param("moisture_decay", 0.995) as f32;
|
||||
fn update_moisture_wind(&mut self, grid: &mut GridState, dt: f32) {
|
||||
let transport_rate = (self.get_param("moisture_transport", 0.15) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
let decay = (self.get_param("moisture_decay", 0.995) as f32).powf(dt);
|
||||
let rain_shadow_block = self.get_param("mountain_rain_shadow_block", 0.9) as f32;
|
||||
let atmo_loss = self.get_param("atmospheric_loss_rate", 0.0003) as f32;
|
||||
let atmo_loss = (self.get_param("atmospheric_loss_rate", 0.0003) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
|
||||
let w = grid.width;
|
||||
let h = grid.height;
|
||||
|
|
@ -410,8 +416,8 @@ impl ClimatePhysics {
|
|||
|
||||
// ── Step 5: River moisture transport ──────────────────────────────────
|
||||
|
||||
fn update_moisture_rivers(&self, grid: &mut GridState) {
|
||||
let transport_rate = self.get_param("river_moisture_transport", 0.075) as f32;
|
||||
fn update_moisture_rivers(&self, grid: &mut GridState, dt: f32) {
|
||||
let transport_rate = (self.get_param("river_moisture_transport", 0.075) as f32 * dt).min(MAX_SAFE_RATE);
|
||||
let w = grid.width;
|
||||
let h = grid.height;
|
||||
let n = grid.tiles.len();
|
||||
|
|
@ -804,9 +810,36 @@ impl ClimatePhysics {
|
|||
}
|
||||
}
|
||||
|
||||
/// Per-tile albedo cache (populated by `rebuild_tile_cache`).
|
||||
pub fn tile_albedo_ref(&self) -> &[f32] {
|
||||
&self.tile_albedo
|
||||
}
|
||||
|
||||
/// Per-tile evapotranspiration cache.
|
||||
pub fn tile_evapotrans_ref(&self) -> &[f32] {
|
||||
&self.tile_evapotrans
|
||||
}
|
||||
|
||||
/// Per-tile flag bits cache.
|
||||
pub fn tile_flags_ref(&self) -> &[u8] {
|
||||
&self.tile_flags
|
||||
}
|
||||
|
||||
/// Steps 5–12: all sequential steps after the two diffusion passes.
|
||||
pub fn step_remaining(&self, grid: &mut GridState, dt: f32) {
|
||||
self.update_moisture_rivers(grid, dt);
|
||||
self.update_lake_evaporation(grid, grid.width, grid.height);
|
||||
self.update_deep_earth_water(grid);
|
||||
self.update_precipitation(grid);
|
||||
self.update_surface_runoff(grid);
|
||||
self.check_terrain_evolution(grid);
|
||||
self.compute_global_stats(grid);
|
||||
self.clear_magic_deltas(grid);
|
||||
}
|
||||
|
||||
// ── Helpers ──────────────────────────────────────────────────────────
|
||||
|
||||
fn get_param(&self, key: &str, default: f64) -> f64 {
|
||||
pub fn get_param(&self, key: &str, default: f64) -> f64 {
|
||||
self.params
|
||||
.get(key)
|
||||
.and_then(|v| v.as_f64())
|
||||
|
|
@ -891,7 +924,7 @@ mod tests {
|
|||
let old_temp_0 = grid.tiles[0].temperature;
|
||||
let old_temp_1 = grid.tiles[1].temperature;
|
||||
|
||||
physics.process_step(&mut grid, 1, 42);
|
||||
physics.process_step(&mut grid, 1, 42, 1.0);
|
||||
|
||||
// Temperature should change due to solar forcing and wind transport
|
||||
assert_ne!(grid.tiles[0].temperature, old_temp_0);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue