diff --git a/Cargo.toml b/Cargo.toml index 07b6d6b..d006db6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -34,8 +34,6 @@ utils = { path = "crates/utils" } bevy = { version = "0.18", features = ["dynamic_linking"] } -glam = "0.29.2" - # Serde for serialization and deserialization serde = "1.0" @@ -44,3 +42,6 @@ serde_json = "1.0" # rand for picking random values rand = "0.9" + +# bevy_egui for UI controls +bevy_egui = "0.39" diff --git a/crates/energy/README.md b/crates/energy/README.md index b486663..d3b51b5 100644 --- a/crates/energy/README.md +++ b/crates/energy/README.md @@ -2,16 +2,18 @@ Energy conservation, thermodynamics, electromagnetism, and wave mechanics for LP physics simulations. -## Scope and units +## Core Principles + - Tracks energy in joules via an accounting ledger; conservation is enforced where modeled. - Thermodynamics, EM, and waves modules are present but remain partial implementations. - Uses consistent sim units (SI-style); couple to time/space steps used by the broader sim. ## Scope & Limits + - LP-0 thermodynamics and EM use explicit approximations: pairwise interactions, cutoffs, and quasi-static assumptions for performance. - Wave solvers use finite differences and simplified damping models; energy coupling is partial. - Thermal energy uses constant heat capacity and clamps temperatures at 0 K; phase change and EOS are deferred. -## Not yet implemented -- **Wave/EM energy accounting**: Wave damping exists but not coupled to energy ledger; EM fields exist but no work/energy tracking. -- **EOS, convection, latent heat, phase transitions**: Deferred to matter/MPM coupling (see tmp/GD/Systems_Overview.md). +## Status + +Early-stage. Partially integrated with forces crate; wave/EM energy accounting not yet coupled to ledger. Blocked on matter/MPM stabilization for phase transitions and EOS. diff --git a/crates/energy/src/conservation.rs b/crates/energy/src/conservation.rs index 6f3e6f3..50913ee 100644 --- a/crates/energy/src/conservation.rs +++ b/crates/energy/src/conservation.rs @@ -341,9 +341,9 @@ impl Plugin for EnergyConservationPlugin { .init_resource::() // Add event channel .add_message::() - // Add systems with explicit ordering + // Track energy in FixedUpdate to match physics integration schedule .add_systems( - Update, + FixedUpdate, ( initialize_energy_balance, ApplyDeferred, diff --git a/crates/energy/src/electromagnetism/charges.rs b/crates/energy/src/electromagnetism/charges.rs index 6308762..acbc3ba 100644 --- a/crates/energy/src/electromagnetism/charges.rs +++ b/crates/energy/src/electromagnetism/charges.rs @@ -4,14 +4,25 @@ //! Future: Will be replaced by grid-based Poisson solve (ρ → φ → E). //! //! Physics: F = k·q₁·q₂/r² (Coulomb's law) -//! Complexity: O(N) with SpatialGrid +//! Complexity: Candidate lookup via UnifiedSpatialIndex backend (uniform cell field or hierarchy) //! Conservation: Force-only (EM potential energy = 0 for LP-0) +//! +//! FIXME: Coulomb singularity at r, + sorted_entities: Vec, + neighbor_candidates: Vec, +} + impl Default for SofteningLength { fn default() -> Self { Self { value: 0.01 } @@ -69,9 +87,10 @@ impl Default for SofteningLength { pub struct CoulombConfig { /// Coulomb constant k = 1/(4πε₀). /// - /// **UNITS**: Newtons × meters² / Coulombs² (N·m²/C²) - /// **Default**: 8.99×10⁹ N·m²/C² (vacuum permittivity ε₀ = 8.854×10⁻¹² F/m) - /// **IRL PHYSICS**: k = 8.987551792×10⁹ N·m²/C² (exact) + /// **UNITS**: Simulation units (not IRL N·m²/C²) + /// **Default**: 1.0 (scaled for simulation, tunable via UI controls) + /// **IRL PHYSICS**: k = 8.987551792×10⁹ N·m²/C² (not used in LP) + /// **NOTE**: LP uses simulation-scaled values to keep EM forces visible but non-dominant vs gravity pub coulomb_constant: f32, /// Cutoff radius for Coulomb interactions. @@ -93,7 +112,9 @@ impl Default for CoulombConfig { fn default() -> Self { let cutoff = 20.0; Self { - coulomb_constant: 8.99e9, // N⋅m²/C² + // Simulation-scaled constant (not IRL 8.99e9) + // Scaled to produce visible but non-dominant EM forces vs gravity + coulomb_constant: 1.0, // Simulation units (tune via egui controls) cutoff_radius: cutoff, switch_on_radius: 0.8 * cutoff, } @@ -115,7 +136,7 @@ pub fn mark_charged_entities_spatially_indexed( /// Apply Coulomb electrostatic forces between charged particles. /// -/// **LP-0 SCAFFOLDING**: Pairwise particle-particle Coulomb interactions (O(N) via spatial hash grid). +/// **LP-0 SCAFFOLDING**: Pairwise particle-particle Coulomb interactions. /// **TEMPORARY**: Will be replaced by grid-based Poisson solve (ρ → φ → E) in LP-1. /// /// **PHYSICS**: Coulomb's Law F = k·q₁·q₂/r² (Newtons) @@ -129,11 +150,11 @@ pub fn mark_charged_entities_spatially_indexed( /// - Cutoff radius: 20m default (performance hack, IRL Coulomb has infinite range in vacuum) /// - Softening: 0.01m default (singularity avoidance for r→0) /// - Potential energy: Not tracked (force-only, PE = 0 in LP-0) -/// - Pair-once guarantee: Only processes pairs where entity_b.index() > entity_a.index() to avoid double-counting +/// - Pair-once guarantee: Only processes pairs where entity_b.id > entity_a.id to avoid double-counting /// /// **CONSERVATION**: Momentum conserved (F_ab = -F_ba, Newton's 3rd law). /// Energy NOT conserved (PE missing from accounting). -pub fn apply_coulomb_pairwise_forces( +pub(crate) fn apply_coulomb_pairwise_forces( mut charges: Query<( Entity, &Charge, @@ -143,12 +164,15 @@ pub fn apply_coulomb_pairwise_forces( )>, index: Res, config: Res, + determinism: Res, + mut ctx: Local, ) { // **LP-0 SCAFFOLDING**: Pairwise particle-particle Coulomb forces. // Future: Grid-based Poisson solve (ρ → φ → E). - // Stage charges into map to avoid nested query - let mut charge_data: HashMap = HashMap::new(); + // Reuse staging buffers across frames to avoid per-frame allocation churn. + let estimated = charges.iter().len(); + prepare_staging_map(&mut ctx.charge_data, estimated); for (entity, charge, trans, softening, _) in charges.iter() { let pos = trans.translation.truncate(); @@ -171,54 +195,81 @@ pub fn apply_coulomb_pairwise_forces( } }; - charge_data.insert(entity, (charge.value, pos, soft.value)); + ctx.charge_data + .insert(entity, (charge.value, pos, soft.value)); } - // Iterate pairs via UnifiedSpatialIndex - for (entity_a, (charge_a, pos_a, soft_a)) in charge_data.iter() { - // Find neighbors within cutoff using UnifiedSpatialIndex (O(N) average) - for entity_b in index.query_radius(*pos_a, config.cutoff_radius) { - // **Pair-once guarantee**: Only process pairs where B > A - if entity_b.index() <= entity_a.index() { - continue; - } + // Deterministic outer iteration: sort entities by stable id + let staged_entities: Vec = ctx.charge_data.keys().copied().collect(); + prepare_sorted_entities_from_keys(&mut ctx.sorted_entities, staged_entities); - // Get data for entity B from staged map - let Some((charge_b, pos_b, soft_b)) = charge_data.get(&entity_b) else { - continue; - }; + // Iterate pairs via UnifiedSpatialIndex. + let sorted_entities = std::mem::take(&mut ctx.sorted_entities); + let charge_data = std::mem::take(&mut ctx.charge_data); + for &entity_a in &sorted_entities { + let (charge_a, pos_a, soft_a) = charge_data[&entity_a]; + // Find neighbors within cutoff using UnifiedSpatialIndex backend. + for_each_neighbor_candidate( + &index, + pos_a, + config.cutoff_radius, + determinism.strict_neighbor_order, + &mut ctx.neighbor_candidates, + |entity_b| { + // **Pair-once guarantee**: Only process pairs where B > A + if !is_forward_entity_pair(entity_a, entity_b) { + return; + } - let r_vec = *pos_b - *pos_a; - let r = r_vec.length(); + // Get data for entity B from staged map + let Some((charge_b, pos_b, soft_b)) = charge_data.get(&entity_b) else { + return; + }; - // Property-based softening: use max of the two particles' softening lengths - let softening = soft_a.max(*soft_b); - if r < softening || r >= config.cutoff_radius { - continue; - } + let r_vec = *pos_b - pos_a; + let r = r_vec.length(); - // Coulomb force: F_on_A = -k·q₁·q₂/r² · r̂_AB - // Same-sign charges (k_qq > 0) → repulsive (force pushes A away from B) - // Opposite-sign charges (k_qq < 0) → attractive (force pulls A toward B) - let k_qq = config.coulomb_constant * charge_a * charge_b; - let force_magnitude = k_qq / r.powi(2); - let force_bare = -(force_magnitude / r) * r_vec; - - // Apply C¹ force-switch for smooth cutoff - let switch = force_switch(r, config.switch_on_radius, config.cutoff_radius); - let force_2d = force_bare * switch; - let force = force_2d.extend(0.0); // Convert to Vec3 for AppliedForce - - // Apply forces symmetrically (Newton's 3rd law) - if let Ok((_, _, _, _, mut force_a)) = charges.get_mut(*entity_a) { - force_a.force += force; - } - if let Ok((_, _, _, _, mut force_b)) = charges.get_mut(entity_b) { - force_b.force -= force; // F_ba = -F_ab - } + // Skip if beyond cutoff (but apply softened force within softening zone) + if r >= config.cutoff_radius { + return; + } - // **LP-0**: EM potential energy = 0 (force-only). - // Future: Track U(r) = integral of switched force for energy conservation. - } + // Property-based softening: use max of the two particles' softening lengths + let softening = soft_a.max(*soft_b); + + // Coulomb force: F_on_A = -k·q₁·q₂/r² · r̂_AB + // Same-sign charges (k_qq > 0) → repulsive (force pushes A away from B) + // Opposite-sign charges (k_qq < 0) → attractive (force pulls A toward B) + let k_qq = config.coulomb_constant * charge_a * charge_b; + + // Apply softened force using 1/(r² + softening²)^1.5 near singularity + // This smoothly reduces force as r→0 instead of disabling it + let r_softened = (r * r + softening * softening).sqrt(); + let force_magnitude = k_qq / (r_softened * r_softened * r_softened); + let force_bare = if r > 1e-6 { + -(force_magnitude / r) * r_vec + } else { + Vec2::ZERO // Avoid division by zero at r≈0 + }; + + // Apply C¹ force-switch for smooth cutoff + let switch = force_switch(r, config.switch_on_radius, config.cutoff_radius); + let force_2d = force_bare * switch; + let force = force_2d.extend(0.0); // Convert to Vec3 for AppliedForce + + // Apply forces symmetrically (Newton's 3rd law) + if let Ok((_, _, _, _, mut force_a)) = charges.get_mut(entity_a) { + force_a.force += force; + } + if let Ok((_, _, _, _, mut force_b)) = charges.get_mut(entity_b) { + force_b.force -= force; // F_ba = -F_ab + } + + // **LP-0**: EM potential energy = 0 (force-only). + // Future: Track U(r) = integral of switched force for energy conservation. + }, + ); } + ctx.charge_data = charge_data; + ctx.sorted_entities = sorted_entities; } diff --git a/crates/energy/src/electromagnetism/mod.rs b/crates/energy/src/electromagnetism/mod.rs index da7c450..b8fe654 100644 --- a/crates/energy/src/electromagnetism/mod.rs +++ b/crates/energy/src/electromagnetism/mod.rs @@ -27,9 +27,9 @@ impl Plugin for ElectromagnetismPlugin { charges::mark_charged_entities_spatially_indexed .in_set(SpatialIndexSet::InjectMarkers), ) - // Coulomb forces in Update, in force accumulation, after gravity + // Coulomb forces in FixedUpdate, in force accumulation, after gravity .add_systems( - Update, + FixedUpdate, charges::apply_coulomb_pairwise_forces .in_set(PhysicsSet::AccumulateForces) .after(GravitySet::NBodyGravity), diff --git a/crates/energy/src/lib.rs b/crates/energy/src/lib.rs index b4b0036..001a0f6 100644 --- a/crates/energy/src/lib.rs +++ b/crates/energy/src/lib.rs @@ -1,5 +1,6 @@ pub mod conservation; pub mod electromagnetism; +pub(crate) mod pairwise; pub mod thermodynamics; pub mod waves; @@ -7,6 +8,7 @@ use bevy::prelude::*; pub use conservation::EnergyConservationPlugin; pub use electromagnetism::ElectromagnetismPlugin; +pub use pairwise::PairwiseDeterminismConfig; pub use thermodynamics::ThermodynamicsPlugin; pub use waves::WavesPlugin; @@ -102,6 +104,8 @@ pub struct EnergyPlugin; impl Plugin for EnergyPlugin { fn build(&self, app: &mut App) { app.register_type::() + .init_resource::() + .register_type::() .add_plugins(EnergyConservationPlugin) .add_plugins(ThermodynamicsPlugin) .add_plugins(ElectromagnetismPlugin) @@ -111,6 +115,7 @@ impl Plugin for EnergyPlugin { pub mod prelude { pub use super::{EnergySystem, EnergyTransferError}; + pub use crate::PairwiseDeterminismConfig; pub use crate::conservation::{ EnergyBalance, EnergyConservationPlugin, EnergyConservationTracker, EnergyDriftMonitor, diff --git a/crates/energy/src/pairwise.rs b/crates/energy/src/pairwise.rs new file mode 100644 index 0000000..5dc6ea6 --- /dev/null +++ b/crates/energy/src/pairwise.rs @@ -0,0 +1,68 @@ +use bevy::prelude::*; +use std::collections::HashMap; +use std::hash::Hash; +use utils::UnifiedSpatialIndex; + +/// Controls determinism level for pairwise neighbor emission. +#[derive(Resource, Debug, Clone, Reflect)] +#[reflect(Resource)] +pub struct PairwiseDeterminismConfig { + /// When true, neighbor candidates are collected and sorted by entity id before processing. + /// This enables exact replay determinism at additional per-frame cost. + pub strict_neighbor_order: bool, +} + +impl Default for PairwiseDeterminismConfig { + fn default() -> Self { + Self { + strict_neighbor_order: false, + } + } +} + +/// Prepare a staging map for a new frame without reallocating from scratch. +pub(crate) fn prepare_staging_map( + map: &mut HashMap, + estimated_items: usize, +) { + map.clear(); + map.reserve(estimated_items); +} + +/// Build deterministic entity iteration order from keys. +pub(crate) fn prepare_sorted_entities_from_keys( + sorted_entities: &mut Vec, + keys: impl IntoIterator, +) { + sorted_entities.clear(); + sorted_entities.extend(keys); + sorted_entities.sort_by_key(|e| e.to_bits()); +} + +/// Pair-once check: returns true only if `b` should be processed when iterating from `a`. +pub(crate) fn is_forward_entity_pair(a: Entity, b: Entity) -> bool { + b.to_bits() > a.to_bits() +} + +/// Emit neighbor candidates in either fast path order or strict deterministic order. +pub(crate) fn for_each_neighbor_candidate( + index: &UnifiedSpatialIndex, + position: Vec2, + radius: f32, + strict_neighbor_order: bool, + scratch: &mut Vec, + mut emit: impl FnMut(Entity), +) { + if strict_neighbor_order { + scratch.clear(); + index.for_each_neighbor_candidate_in_radius(position, radius, |entity| { + scratch.push(entity); + }); + scratch.sort_by_key(|e| e.to_bits()); + for &entity in scratch.iter() { + emit(entity); + } + } else { + index.for_each_neighbor_candidate_in_radius(position, radius, emit); + } +} diff --git a/crates/energy/src/thermodynamics/thermal.rs b/crates/energy/src/thermodynamics/thermal.rs index 21b49a3..b1e3f37 100644 --- a/crates/energy/src/thermodynamics/thermal.rs +++ b/crates/energy/src/thermodynamics/thermal.rs @@ -4,6 +4,10 @@ use std::collections::HashMap; use utils::{SpatialIndexSet, SpatiallyIndexed, UnifiedSpatialIndex, force_switch}; use crate::conservation::{EnergyQuantity, EnergyType}; +use crate::pairwise::{ + PairwiseDeterminismConfig, for_each_neighbor_candidate, is_forward_entity_pair, + prepare_sorted_entities_from_keys, prepare_staging_map, +}; /// Configuration for Fourier conduction system. /// @@ -44,6 +48,34 @@ impl Default for ThermalConductionConfig { } } +/// Lightweight runtime sanity checks for thermal state. +/// +/// These checks are O(N) and intended for realtime runs. +/// They are not full conservation diagnostics. +#[derive(Resource, Debug, Clone, Reflect)] +#[reflect(Resource)] +pub struct ThermalSanityConfig { + /// Enable/disable sanity checks. + pub enabled: bool, + /// Warn if total thermal energy drops below this floor (J). + pub min_total_thermal_energy_j: f32, + /// Warn if total thermal energy grows more than this factor frame-to-frame. + pub max_frame_growth_factor: f32, + /// Warn if any temperature exceeds this bound (K). + pub max_temperature_k: f32, +} + +impl Default for ThermalSanityConfig { + fn default() -> Self { + Self { + enabled: true, + min_total_thermal_energy_j: 1e-3, + max_frame_growth_factor: 2.0, + max_temperature_k: 1e8, + } + } +} + // Physical constants pub const STEFAN_BOLTZMANN: f32 = 5.67e-8; // W/(m²·K⁴) @@ -185,6 +217,14 @@ pub struct ThermalTransferEvent { pub heat_flow: f32, } +#[derive(Default)] +pub(crate) struct ThermalComputeContext { + thermal_data: HashMap, + temp_changes: HashMap, + sorted_entities: Vec, + neighbor_candidates: Vec, +} + /// Mark thermal entities for spatial indexing. /// /// **Phase A2**: Inject SpatiallyIndexed marker for UnifiedSpatialIndex. @@ -200,7 +240,7 @@ fn mark_thermal_entities_spatially_indexed( /// Compute thermal transfer via Fourier's Law of conduction. /// -/// **LP-0 SCAFFOLDING**: Pairwise particle-particle thermal conduction (O(N) via spatial hash grid). +/// **LP-0 SCAFFOLDING**: Pairwise particle-particle thermal conduction. /// **TEMPORARY**: Will be replaced by grid-based diffusion PDE (∇·(k∇T) = ρc_p ∂T/∂t) in LP-1. /// /// **PHYSICS**: Fourier's Law q = k·A·ΔT/d (W), Q = P·dt (J) @@ -222,8 +262,10 @@ pub(crate) fn compute_fourier_conduction( mut commands: Commands, index: Res, config: Res, + determinism: Res, time: Res