From 1414f1e7f00e49eb070bdec4edea5aa9aa16e869 Mon Sep 17 00:00:00 2001 From: Claude Code Date: Tue, 31 Mar 2026 22:47:33 -0700 Subject: [PATCH] =?UTF-8?q?feat(climate):=20=E2=9C=A8=20Implement=20new=20?= =?UTF-8?q?climate=20simulation=20algorithms=20for=20ecological=20and=20ph?= =?UTF-8?q?ysical=20process=20modeling?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Lilith Autocommit --- .../crates/mc-climate/src/ecology.rs | 93 ++++++++++++++----- src/simulator/crates/mc-climate/src/lib.rs | 5 +- .../crates/mc-climate/src/physics.rs | 87 +++++++++++------ 3 files changed, 134 insertions(+), 51 deletions(-) diff --git a/src/simulator/crates/mc-climate/src/ecology.rs b/src/simulator/crates/mc-climate/src/ecology.rs index b553881c..186c72c0 100644 --- a/src/simulator/crates/mc-climate/src/ecology.rs +++ b/src/simulator/crates/mc-climate/src/ecology.rs @@ -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"); + } } diff --git a/src/simulator/crates/mc-climate/src/lib.rs b/src/simulator/crates/mc-climate/src/lib.rs index 5cb7e95b..1c499f33 100644 --- a/src/simulator/crates/mc-climate/src/lib.rs +++ b/src/simulator/crates/mc-climate/src/lib.rs @@ -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, +}; diff --git a/src/simulator/crates/mc-climate/src/physics.rs b/src/simulator/crates/mc-climate/src/physics.rs index aad96c28..3ca4a762 100644 --- a/src/simulator/crates/mc-climate/src/physics.rs +++ b/src/simulator/crates/mc-climate/src/physics.rs @@ -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);