From f984a6881e9e1e933a88717db70621acc189ee6b Mon Sep 17 00:00:00 2001 From: M1thieu <18742831+M1thieu@users.noreply.github.com> Date: Thu, 12 Feb 2026 08:24:43 +0100 Subject: [PATCH 1/7] feat(utils): add adaptive neighbor index refactor(energy): deterministic pair loops refactor(forces): stage gravity pairs and document realtime limits --- crates/energy/src/electromagnetism/charges.rs | 35 +- crates/energy/src/thermodynamics/thermal.rs | 47 +- crates/forces/README.md | 1 + crates/forces/src/core/gravity.rs | 148 +++-- crates/utils/README.md | 26 +- crates/utils/src/lib.rs | 12 +- crates/utils/src/spatial/unified.rs | 568 ++++++++++++++++-- 7 files changed, 702 insertions(+), 135 deletions(-) diff --git a/crates/energy/src/electromagnetism/charges.rs b/crates/energy/src/electromagnetism/charges.rs index 6308762..25d3f3c 100644 --- a/crates/energy/src/electromagnetism/charges.rs +++ b/crates/energy/src/electromagnetism/charges.rs @@ -4,7 +4,7 @@ //! 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) use bevy::prelude::*; @@ -115,7 +115,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,7 +129,7 @@ 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). @@ -174,27 +174,34 @@ pub fn apply_coulomb_pairwise_forces( 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) { + // Deterministic outer iteration: sort entities by stable id + let mut sorted_entities: Vec = charge_data.keys().copied().collect(); + sorted_entities.sort_by_key(|e| e.to_bits()); + + // Iterate pairs via UnifiedSpatialIndex. + // TODO(LP-0 determinism mode): For exact replay builds, collect and sort emitted neighbors + // by entity id before accumulation. Default path stays allocation-free for realtime. + for entity_a in sorted_entities { + let (charge_a, pos_a, soft_a) = charge_data[&entity_a]; + // Find neighbors within cutoff using UnifiedSpatialIndex backend + index.for_each_neighbor_candidate_in_radius(pos_a, config.cutoff_radius, |entity_b| { // **Pair-once guarantee**: Only process pairs where B > A - if entity_b.index() <= entity_a.index() { - continue; + if entity_b.to_bits() <= entity_a.to_bits() { + return; } // Get data for entity B from staged map let Some((charge_b, pos_b, soft_b)) = charge_data.get(&entity_b) else { - continue; + return; }; - let r_vec = *pos_b - *pos_a; + let r_vec = *pos_b - pos_a; let r = r_vec.length(); // 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; + return; } // Coulomb force: F_on_A = -k·q₁·q₂/r² · r̂_AB @@ -210,7 +217,7 @@ pub fn apply_coulomb_pairwise_forces( 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) { + 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) { @@ -219,6 +226,6 @@ pub fn apply_coulomb_pairwise_forces( // **LP-0**: EM potential energy = 0 (force-only). // Future: Track U(r) = integral of switched force for energy conservation. - } + }); } } diff --git a/crates/energy/src/thermodynamics/thermal.rs b/crates/energy/src/thermodynamics/thermal.rs index 21b49a3..77c24d6 100644 --- a/crates/energy/src/thermodynamics/thermal.rs +++ b/crates/energy/src/thermodynamics/thermal.rs @@ -200,7 +200,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) @@ -292,28 +292,35 @@ pub(crate) fn compute_fourier_conduction( let cutoff_radius = config.cutoff_radius; let switch_on_radius = config.switch_on_radius; - // Iterate pairs via UnifiedSpatialIndex - for (entity_a, (pos_a, temp_a, k_a, capacity_a, radius_a)) in thermal_data.iter() { - // Find neighbors within cutoff using UnifiedSpatialIndex (O(N) average) - for entity_b in index.query_radius(*pos_a, cutoff_radius) { + // Deterministic outer iteration: sort entities by stable id + let mut sorted_entities: Vec = thermal_data.keys().copied().collect(); + sorted_entities.sort_by_key(|e| e.to_bits()); + + // Iterate pairs via UnifiedSpatialIndex. + // TODO(LP-0 determinism mode): For exact replay builds, collect and sort emitted neighbors + // by entity id before accumulation. Default path stays allocation-free for realtime. + for entity_a in sorted_entities { + let (pos_a, temp_a, k_a, capacity_a, radius_a) = thermal_data[&entity_a]; + // Find neighbors within cutoff using UnifiedSpatialIndex backend + index.for_each_neighbor_candidate_in_radius(pos_a, cutoff_radius, |entity_b| { // **Pair-once guarantee**: Only process pairs where B > A - if entity_b.index() <= entity_a.index() { - continue; + if entity_b.to_bits() <= entity_a.to_bits() { + return; } // Get data for entity B from staged map let Some((pos_b, temp_b, k_b, capacity_b, radius_b)) = thermal_data.get(&entity_b) else { - continue; + return; }; - let r_vec = *pos_b - *pos_a; + let r_vec = *pos_b - pos_a; let distance = r_vec.length(); // Property-based softening: use smaller radius as minimum distance let softening = radius_a.min(*radius_b); if distance < softening || distance >= cutoff_radius { - continue; + return; } let temp_diff = temp_a - temp_b; // +: A hotter, -: B hotter @@ -335,7 +342,7 @@ pub(crate) fn compute_fourier_conduction( let power = k_avg * contact_area * temp_diff / distance; // W if !power.is_finite() { - continue; // Skip non-finite power + return; // Skip non-finite power } // Apply C¹ force-switch for smooth cutoff @@ -359,18 +366,18 @@ pub(crate) fn compute_fourier_conduction( let delta_t_b = heat_energy / capacity_b; // Heating if heat_energy > 0 if !delta_t_a.is_finite() || !delta_t_b.is_finite() { - continue; // Skip non-finite temperature changes + return; // Skip non-finite temperature changes } - *temp_changes.entry(*entity_a).or_insert(0.0) += delta_t_a; + *temp_changes.entry(entity_a).or_insert(0.0) += delta_t_a; *temp_changes.entry(entity_b).or_insert(0.0) += delta_t_b; thermal_transfer_events.write(ThermalTransferEvent { - source: *entity_a, + source: entity_a, target: entity_b, heat_flow: power_switched.abs(), }); - } + }); } // Apply temperature changes @@ -448,15 +455,15 @@ fn check_thermal_stability( if dt > max_stable_dt { warn!( - "⚠️ Thermal diffusion may be UNSTABLE!\n\ + "Thermal diffusion may be UNSTABLE!\n\ Current timestep: dt = {:.6} s\n\ Stability limit: dt <= {:.6} s\n\ Thermal diffusivity: α = {:.2e} m²/s\n\ Grid spacing: dx = {:.2} m\n\ - → Simulation may diverge or explode. Consider:\n\ - • Smaller timestep (use fixed timestep plugin)\n\ - • Coarser grid (increase SpatialGrid cell_size)\n\ - • Implicit solver (future MPM work)", + Simulation may diverge or explode. Consider:\n\ + - Smaller timestep (use fixed timestep plugin)\n\ + - Coarser grid (increase SpatialGrid cell_size)\n\ + - Implicit solver (future MPM work)", dt, max_stable_dt, alpha, dx ); } diff --git a/crates/forces/README.md b/crates/forces/README.md index 775411f..080d318 100644 --- a/crates/forces/README.md +++ b/crates/forces/README.md @@ -11,6 +11,7 @@ Gravitational mechanics and Newton's laws for physics-based force systems. ## Scope & Limits - LP-0 integrates F = ma with explicit/symplectic Euler and optional acceleration clamps for stability; this is a numerical method, not a physical law. - Gravity defaults to a sim-tuned constant and softened inverse-square forces; Barnes-Hut is an approximation for large N. +- Mutual gravity mode is exact pairwise O(N^2); treat ~100 active sources as the LP-0 realtime comfort range. - Potential energy and work accounting are not yet tracked; conservation diagnostics are partial. ## Conservation status diff --git a/crates/forces/src/core/gravity.rs b/crates/forces/src/core/gravity.rs index 1091fd4..24f48ef 100644 --- a/crates/forces/src/core/gravity.rs +++ b/crates/forces/src/core/gravity.rs @@ -5,6 +5,8 @@ use bevy::prelude::*; // Simulation constants // NOTE: Sim-tuned G (not SI 6.67e-11). Pixel→meter mapping TBD; adjust after scale is fixed. pub const DEFAULT_GRAVITATIONAL_CONSTANT: f32 = 0.1; +/// Practical LP-0 guideline for exact mutual O(N^2) gravity in realtime. +pub const MUTUAL_REALTIME_BODY_LIMIT: usize = 100; /// Resource for gravity simulation parameters #[derive(Resource, Clone, Debug)] @@ -97,6 +99,50 @@ pub struct GravitySource; #[reflect(Component)] pub struct MassiveBody; +#[derive(Clone, Copy, Debug)] +struct GravityBody { + entity: Entity, + position: Vec3, + mass: f32, +} + +fn stage_gravity_bodies( + query: &Query<(Entity, &Transform, &Mass), With>, +) -> Vec { + query + .iter() + .map(|(entity, transform, mass)| GravityBody { + entity, + position: transform.translation, + mass: mass.value, + }) + .collect() +} + +fn pair_force_vector( + source_pos: Vec3, + source_mass: f32, + affected_pos: Vec3, + affected_mass: f32, + gravitational_constant: f32, + softening_squared: f32, +) -> Option { + let direction = source_pos - affected_pos; + let distance_squared = direction.length_squared(); + if distance_squared <= f32::EPSILON { + return None; + } + + let softened_distance_squared = distance_squared + softening_squared; + let force_magnitude = + gravitational_constant * source_mass * affected_mass / softened_distance_squared; + if !force_magnitude.is_finite() { + return None; + } + + Some(direction.normalize() * force_magnitude) +} + // Barnes-Hut spatial partitioning mod spatial { use bevy::prelude::*; @@ -328,40 +374,45 @@ pub fn calculate_mutual_gravitational_attraction( let softening_squared = gravity_params.softening * gravity_params.softening; let gravitational_constant = gravity_params.gravitational_constant; - // Collect force updates to avoid mutable aliasing - let mut force_updates: Vec<(Entity, Vec3)> = Vec::new(); - - // Read-only iteration to compute forces - let entities: Vec<_> = query + // Stage source data once (particular-style storage split). + let mut bodies: Vec<_> = query .iter() - .map(|(e, t, m, _)| (e, t.translation, m.value)) + .map(|(entity, transform, mass, _)| GravityBody { + entity, + position: transform.translation, + mass: mass.value, + }) .collect(); - - for i in 0..entities.len() { - for j in (i + 1)..entities.len() { - let (entity1, pos1, mass1) = entities[i]; - let (entity2, pos2, mass2) = entities[j]; - - let direction = pos2 - pos1; - let distance_squared = direction.length_squared(); - if distance_squared <= f32::EPSILON { + // Stable ordering keeps accumulation deterministic for replay/debug. + bodies.sort_by_key(|body| body.entity.to_bits()); + let mut accumulated_forces = vec![Vec3::ZERO; bodies.len()]; + // TODO(LP-1): For mutual mode above MUTUAL_REALTIME_BODY_LIMIT, add an approximate + // path (e.g. Barnes-Hut variant with symmetry corrections) to keep realtime budgets. + + // Pair-once pass, add opposite forces by index. + for i in 0..bodies.len() { + let body_a = bodies[i]; + for j in (i + 1)..bodies.len() { + let body_b = bodies[j]; + let Some(force_on_a) = pair_force_vector( + body_b.position, + body_b.mass, + body_a.position, + body_a.mass, + gravitational_constant, + softening_squared, + ) else { continue; - } - - let softened_distance_squared = distance_squared + softening_squared; - let force_magnitude = - gravitational_constant * mass1 * mass2 / softened_distance_squared; - let force_vector = direction.normalize() * force_magnitude; - - force_updates.push((entity1, force_vector)); - force_updates.push((entity2, -force_vector)); + }; + accumulated_forces[i] += force_on_a; + accumulated_forces[j] -= force_on_a; // Newton's 3rd law } } - // Apply accumulated forces - for (entity, force) in force_updates { - if let Ok((_, _, _, mut applied_force)) = query.get_mut(entity) { - applied_force.force += force; + // Single writeback pass to ECS. + for (index, body) in bodies.iter().enumerate() { + if let Ok((_, _, _, mut applied_force)) = query.get_mut(body.entity) { + applied_force.force += accumulated_forces[index]; } } } @@ -377,31 +428,29 @@ pub fn calculate_gravitational_attraction( let softening_squared = gravity_params.softening * gravity_params.softening; let gravitational_constant = gravity_params.gravitational_constant; - let sources: Vec<(Entity, Vec3, f32)> = query - .iter() - .map(|(e, t, m)| (e, t.translation, m.value)) - .collect(); + let sources = stage_gravity_bodies(&query); + // Safe to parallelize: each worker mutates a distinct affected entity while sources stay read-only. affected_query.par_iter_mut().for_each( |(affected_entity, affected_transform, affected_mass, mut force)| { let affected_pos = affected_transform.translation; - for &(source_entity, source_pos, source_mass) in &sources { - if source_entity == affected_entity { + for source in &sources { + if source.entity == affected_entity { continue; } - let direction = source_pos - affected_pos; - let distance_squared = direction.length_squared(); - let softened_distance_squared = distance_squared + softening_squared; - let force_magnitude = gravitational_constant * source_mass * affected_mass.value - / softened_distance_squared; - - if !force_magnitude.is_finite() { - continue; // Skip non-finite gravity force - } - - force.force += direction.normalize() * force_magnitude; + let Some(force_vector) = pair_force_vector( + source.position, + source.mass, + affected_pos, + affected_mass.value, + gravitational_constant, + softening_squared, + ) else { + continue; + }; + force.force += force_vector; } }, ); @@ -422,13 +471,14 @@ pub fn calculate_barnes_hut_attraction( return; } - let bodies: Vec<(Entity, Vec3, f32)> = query + let bodies = stage_gravity_bodies(&query); + + let body_data: Vec<_> = bodies .iter() - .map(|(e, t, m)| (e, t.translation, m.value)) + .map(|body| (body.entity, body.position, body.mass)) .collect(); - let quadtree = spatial::Quadtree::from_bodies( - &bodies, + &body_data, gravity_params.barnes_hut_max_depth, gravity_params.barnes_hut_max_bodies_per_node, ); diff --git a/crates/utils/README.md b/crates/utils/README.md index 5ec5d2c..653a7fc 100644 --- a/crates/utils/README.md +++ b/crates/utils/README.md @@ -4,13 +4,15 @@ Shared utilities and optimizations. ## Modules -- `spatial::SpatialGrid` - Sparse spatial hash grid for neighbor queries +- `spatial::UnifiedSpatialIndex` - Hybrid spatial index (grid hash + optional tree backend) +- `spatial::SpatialGrid` - Sparse spatial hash grid used by UnifiedSpatialIndex - `pool::EntityPool` - Entity recycling for spawn/despawn scenarios ## Scope & Limits -- SpatialGrid and UnifiedSpatialIndex are 2D, point-entity hash grids; large bodies or 3D need other structures. -- Radius queries return candidate sets from grid cells; exact distance filtering remains a caller responsibility. +- UnifiedSpatialIndex is 2D and point-entity based; large bodies or 3D need other structures. +- Radius queries return candidate sets from the active backend; exact distance filtering remains a caller responsibility. - ECS-based helpers are tuned for LP-0 scale; high-N systems will move to SoA/MPM paths. +- LP-0 uses 2D neighbor search. If LP moves core continuum simulation to 3D MPM, this layer needs a 3D upgrade path. ## Usage @@ -22,9 +24,10 @@ app.add_plugins(utils::UtilsPlugin); ### Current Implementation -**SpatialGrid:** Incremental updates for ECS queries (thermal, waves, AI) -- Handles 100-10k entities efficiently -- Prefer stable cell sizes to avoid churn; reuse a single grid per query domain to reduce allocations. +**UnifiedSpatialIndex:** Runtime backend policy for ECS queries (thermal, EM, AI) +- `UniformCellField`: fastest for dense, highly dynamic domains. +- `HierarchicalEnvelopeField`: useful for sparse/non-uniform distributions. +- `Adaptive`: switches backend from observed density with cooldown to avoid thrashing. **EntityPool (Hard Pool):** Strips components on release - Use for moderate spawn rates (<100/frame) @@ -46,7 +49,16 @@ struct ParticleSystem { } ``` +### LP-0 Spatial TODO (Non-Overkill) + +These are the next practical upgrades from the current hybrid index work: + +- Add `for_each_in_aabb` query path for batch neighborhood pulls (important for later MPM coupling). +- Add an optional strict-determinism query mode that emits neighbors sorted by entity id (debug/replay builds). +- Tune tree split/rebuild heuristics using LP scene profiles before adding new backend complexity. +- Keep scope 2D for LP-0; only move to 3D index primitives when MPM integration starts. + ## Architecture -- **Mind (ECS):** Behavior, AI -> SpatialGrid + EntityPool +- **Mind (ECS):** Behavior, AI -> UnifiedSpatialIndex + EntityPool - **Body (SoA):** Physics, chemistry -> custom structures diff --git a/crates/utils/src/lib.rs b/crates/utils/src/lib.rs index 8a04f06..9fd5ef5 100644 --- a/crates/utils/src/lib.rs +++ b/crates/utils/src/lib.rs @@ -5,7 +5,8 @@ pub mod units; use bevy::prelude::*; use spatial::unified::{ - attach_spatial_cells, remove_from_index_on_marker_removed, update_spatial_index, + attach_spatial_cells, refresh_spatial_index_policy, remove_from_index_on_marker_removed, + update_spatial_index, }; /// Plugin for registering shared utility components @@ -16,8 +17,11 @@ impl Plugin for UtilsPlugin { app // Spatial indexing .init_resource::() + .init_resource::() .register_type::() .register_type::() + .register_type::() + .register_type::() .register_type::() .register_type::() // Units and scale @@ -39,6 +43,7 @@ impl Plugin for UtilsPlugin { attach_spatial_cells, update_spatial_index, remove_from_index_on_marker_removed, + refresh_spatial_index_policy, ) .chain() .in_set(spatial::unified::SpatialIndexSet::Maintain), @@ -53,5 +58,8 @@ impl Plugin for UtilsPlugin { pub use cutoff::force_switch; pub use pool::{EntityPool, Pooled}; pub use spatial::grid::{GridCell, SpatialGrid}; -pub use spatial::unified::{SpatialCell, SpatialIndexSet, SpatiallyIndexed, UnifiedSpatialIndex}; +pub use spatial::unified::{ + NeighborSearchConfig, NeighborSearchMode, SpatialCell, SpatialIndexSet, SpatiallyIndexed, + UnifiedSpatialIndex, +}; pub use units::{PhysicsScale, physics_to_render, render_to_physics}; diff --git a/crates/utils/src/spatial/unified.rs b/crates/utils/src/spatial/unified.rs index 31d1aaf..14c9be6 100644 --- a/crates/utils/src/spatial/unified.rs +++ b/crates/utils/src/spatial/unified.rs @@ -2,16 +2,17 @@ //! //! - No crate cycles: this module knows nothing about Charge/Temperature/Mass //! - Entities opt-in via `SpatiallyIndexed` marker (physics crates inject) -//! - Membership tracking: entity→cell map is authoritative +//! - Membership tracking: entity->position map is authoritative use bevy::prelude::*; +use std::cmp::Ordering; use std::collections::HashMap; use super::grid::SpatialGrid; /// System sets for spatial index maintenance. /// -/// **Execution order**: InjectMarkers → Maintain (in PreUpdate) +/// **Execution order**: InjectMarkers -> Maintain (in PreUpdate) #[derive(SystemSet, Debug, Hash, PartialEq, Eq, Clone)] pub enum SpatialIndexSet { /// Marker injection (physics crates insert SpatiallyIndexed) @@ -34,84 +35,460 @@ pub struct SpatialCell { pub cell: (i32, i32), } +/// Neighbor-search backend mode. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Reflect)] +pub enum NeighborSearchMode { + /// Uniform cell field broadphase (hash grid). + UniformCellField, + /// LP-native tree index (AABB tree with bulk rebuild from tracked points). + HierarchicalEnvelopeField, + /// Runtime backend selection from observed sparsity. + Adaptive, +} + +/// Runtime configuration for unified spatial indexing. +#[derive(Resource, Debug, Clone, Reflect)] +#[reflect(Resource)] +pub struct NeighborSearchConfig { + /// Default cell size in meters (used for grid indexing and density estimation). + pub cell_size_meters: f32, + /// Backend selection policy. + pub mode: NeighborSearchMode, + /// Minimum entity count before Adaptive mode considers switching to tree mode. + pub adaptive_min_entities_for_hierarchy: usize, + /// If entities/cell is lower than this threshold, Adaptive mode prefers tree mode. + pub adaptive_sparse_entities_per_cell_threshold: f32, + /// Cooldown to avoid backend thrashing in Adaptive mode. + pub adaptive_switch_cooldown_frames: u32, + /// Maximum points per tree leaf. + pub hierarchy_leaf_capacity: usize, +} + +impl Default for NeighborSearchConfig { + fn default() -> Self { + Self { + cell_size_meters: 50.0, + mode: NeighborSearchMode::UniformCellField, + adaptive_min_entities_for_hierarchy: 1000, + adaptive_sparse_entities_per_cell_threshold: 0.35, + adaptive_switch_cooldown_frames: 120, + hierarchy_leaf_capacity: 24, + } + } +} + +#[derive(Default, Clone, Copy, Debug, PartialEq, Eq)] +enum BackendStorage { + #[default] + Grid, + Tree, +} + +#[derive(Clone, Copy, Debug)] +struct TreeAabb { + min: Vec2, + max: Vec2, +} + +impl TreeAabb { + fn from_points(points: &[(Entity, Vec2)]) -> Self { + let mut min = Vec2::splat(f32::MAX); + let mut max = Vec2::splat(f32::MIN); + for (_, p) in points { + min.x = min.x.min(p.x); + min.y = min.y.min(p.y); + max.x = max.x.max(p.x); + max.y = max.y.max(p.y); + } + Self { min, max } + } + + fn merge(a: Self, b: Self) -> Self { + Self { + min: Vec2::new(a.min.x.min(b.min.x), a.min.y.min(b.min.y)), + max: Vec2::new(a.max.x.max(b.max.x), a.max.y.max(b.max.y)), + } + } + + fn distance2_to_point(self, p: Vec2) -> f32 { + let dx = if p.x < self.min.x { + self.min.x - p.x + } else if p.x > self.max.x { + p.x - self.max.x + } else { + 0.0 + }; + let dy = if p.y < self.min.y { + self.min.y - p.y + } else if p.y > self.max.y { + p.y - self.max.y + } else { + 0.0 + }; + dx * dx + dy * dy + } + + fn widest_axis(self) -> usize { + let ext = self.max - self.min; + if ext.x >= ext.y { 0 } else { 1 } + } +} + +#[derive(Clone, Debug)] +enum TreeNode { + Leaf { + aabb: TreeAabb, + points: Vec<(Entity, Vec2)>, + }, + Branch { + aabb: TreeAabb, + left: Box, + right: Box, + }, +} + +impl TreeNode { + fn aabb(&self) -> TreeAabb { + match self { + TreeNode::Leaf { aabb, .. } | TreeNode::Branch { aabb, .. } => *aabb, + } + } + + fn for_each_in_radius(&self, center: Vec2, radius2: f32, emit: &mut impl FnMut(Entity)) { + if self.aabb().distance2_to_point(center) > radius2 { + return; + } + + match self { + TreeNode::Leaf { points, .. } => { + for (entity, p) in points { + if p.distance_squared(center) <= radius2 { + emit(*entity); + } + } + } + TreeNode::Branch { left, right, .. } => { + left.for_each_in_radius(center, radius2, emit); + right.for_each_in_radius(center, radius2, emit); + } + } + } +} + +#[derive(Default)] +struct SpatialTreeIndex { + root: Option>, + leaf_size: usize, +} + +impl SpatialTreeIndex { + fn new(leaf_size: usize) -> Self { + Self { + root: None, + leaf_size: leaf_size.max(4), + } + } + + fn set_leaf_size(&mut self, leaf_size: usize) { + self.leaf_size = leaf_size.max(4); + } + + fn rebuild_from_positions(&mut self, positions: &HashMap) { + if positions.is_empty() { + self.root = None; + return; + } + + let mut points: Vec<(Entity, Vec2)> = positions.iter().map(|(e, p)| (*e, *p)).collect(); + self.root = Some(Self::build_node(&mut points, self.leaf_size)); + } + + fn build_node(points: &mut Vec<(Entity, Vec2)>, leaf_size: usize) -> Box { + let aabb = TreeAabb::from_points(points); + if points.len() <= leaf_size { + return Box::new(TreeNode::Leaf { + aabb, + points: std::mem::take(points), + }); + } + + let axis = aabb.widest_axis(); + points.sort_by(|a, b| { + let va = if axis == 0 { a.1.x } else { a.1.y }; + let vb = if axis == 0 { b.1.x } else { b.1.y }; + match va.partial_cmp(&vb).unwrap_or(Ordering::Equal) { + Ordering::Equal => a.0.to_bits().cmp(&b.0.to_bits()), + ord => ord, + } + }); + + let mid = points.len() / 2; + let mut right_points = points.split_off(mid); + let mut left_points = std::mem::take(points); + + let left = Self::build_node(&mut left_points, leaf_size); + let right = Self::build_node(&mut right_points, leaf_size); + let branch_aabb = TreeAabb::merge(left.aabb(), right.aabb()); + + Box::new(TreeNode::Branch { + aabb: branch_aabb, + left, + right, + }) + } + + fn for_each_in_radius(&self, center: Vec2, radius: f32, mut emit: impl FnMut(Entity)) { + if let Some(root) = &self.root { + root.for_each_in_radius(center, radius * radius, &mut emit); + } + } +} + /// Unified spatial index with correct membership tracking. /// /// **Correctness guarantees**: -/// - On insert/update: entity in exactly one cell -/// - On move: remove from old cell, insert into new cell -/// - No stale membership, no duplicates +/// - On insert/update: entity has exactly one stored position +/// - Query returns candidate entities only; callers apply exact physical filtering +/// - Backend switching preserves all tracked entities #[derive(Resource)] pub struct UnifiedSpatialIndex { grid: SpatialGrid, - entity_cells: HashMap, + tree: SpatialTreeIndex, + backend: BackendStorage, + entity_positions: HashMap, + config: NeighborSearchConfig, + frames_since_switch: u32, + tree_dirty: bool, } impl UnifiedSpatialIndex { pub fn new(cell_size_meters: f32) -> Self { + let config = NeighborSearchConfig { + cell_size_meters, + ..NeighborSearchConfig::default() + }; + Self::from_config(config) + } + + pub fn from_config(config: NeighborSearchConfig) -> Self { Self { - grid: SpatialGrid::new(cell_size_meters), - entity_cells: HashMap::new(), + grid: SpatialGrid::new(config.cell_size_meters), + tree: SpatialTreeIndex::new(config.hierarchy_leaf_capacity), + backend: BackendStorage::Grid, + entity_positions: HashMap::new(), + config, + frames_since_switch: 0, + tree_dirty: false, + } + } + + pub fn set_config(&mut self, config: &NeighborSearchConfig) { + let cell_size_changed = + (self.config.cell_size_meters - config.cell_size_meters).abs() > f32::EPSILON; + self.config = config.clone(); + self.tree.set_leaf_size(self.config.hierarchy_leaf_capacity); + + if cell_size_changed { + self.grid = SpatialGrid::new(self.config.cell_size_meters); + if matches!(self.backend, BackendStorage::Grid) { + for (entity, position) in &self.entity_positions { + self.grid.insert(*entity, *position); + } + } } } - /// Insert entity and return computed cell. + pub fn active_mode(&self) -> NeighborSearchMode { + match self.backend { + BackendStorage::Grid => NeighborSearchMode::UniformCellField, + BackendStorage::Tree => NeighborSearchMode::HierarchicalEnvelopeField, + } + } + + /// Insert entity and return computed grid cell for compatibility with existing callers. pub fn insert(&mut self, entity: Entity, position: Vec2) -> (i32, i32) { - let cell = self.grid.world_to_grid(position); + let old = self.entity_positions.insert(entity, position); - // Handle duplicate insertion (idempotent) - if let Some(old) = self.entity_cells.insert(entity, cell) { - if old != cell { - self.grid.remove_from_cell(entity, old); + match self.backend { + BackendStorage::Grid => { + if let Some(old_pos) = old { + let old_cell = self.grid.world_to_grid(old_pos); + self.grid.remove_from_cell(entity, old_cell); + } + self.grid.insert(entity, position); + } + BackendStorage::Tree => { + self.tree_dirty = true; } } - self.grid.insert_in_cell(entity, cell); - cell + self.grid.world_to_grid(position) } - /// Update entity position and return current cell. + /// Update entity position and return current grid cell. pub fn update(&mut self, entity: Entity, position: Vec2) -> (i32, i32) { - let new_cell = self.grid.world_to_grid(position); - - match self.entity_cells.get(&entity).copied() { - Some(old_cell) if old_cell == new_cell => new_cell, - Some(old_cell) => { - self.grid.move_entity(entity, old_cell, new_cell); - self.entity_cells.insert(entity, new_cell); - new_cell - } - None => { - // First-time tracking (or missed attach) - self.grid.insert_in_cell(entity, new_cell); - self.entity_cells.insert(entity, new_cell); - new_cell + let old = self.entity_positions.insert(entity, position); + + match self.backend { + BackendStorage::Grid => match old { + Some(old_pos) => { + let old_cell = self.grid.world_to_grid(old_pos); + let new_cell = self.grid.world_to_grid(position); + self.grid.move_entity(entity, old_cell, new_cell); + } + None => self.grid.insert(entity, position), + }, + BackendStorage::Tree => { + self.tree_dirty = true; } } + + self.grid.world_to_grid(position) } pub fn remove(&mut self, entity: Entity) { - if let Some(cell) = self.entity_cells.remove(&entity) { - self.grid.remove_from_cell(entity, cell); + if let Some(old_pos) = self.entity_positions.remove(&entity) { + match self.backend { + BackendStorage::Grid => { + let old_cell = self.grid.world_to_grid(old_pos); + self.grid.remove_from_cell(entity, old_cell); + } + BackendStorage::Tree => { + self.tree_dirty = true; + } + } + } + } + + /// Query candidate entities within radius. + /// + /// Returned entities are candidates only; exact distance checks remain a caller responsibility. + pub fn for_each_neighbor_candidate_in_radius( + &self, + position: Vec2, + radius: f32, + mut emit: impl FnMut(Entity), + ) { + debug_assert!(radius >= 0.0, "Negative radius query is undefined"); + + match self.backend { + BackendStorage::Grid => { + for entity in self.grid.get_entities_in_radius(position, radius) { + emit(entity); + } + } + BackendStorage::Tree => self.tree.for_each_in_radius(position, radius, emit), } } - /// Query entities in radius. Domain filtering via staged maps. - pub fn query_radius(&self, position: Vec2, radius: f32) -> impl Iterator + '_ { - self.grid.get_entities_in_radius(position, radius) + /// Query candidate entities within radius into a newly allocated vector. + /// + /// Use `for_each_neighbor_candidate_in_radius` in hot paths to avoid per-query allocations. + pub fn query_radius(&self, position: Vec2, radius: f32) -> Vec { + let mut out = Vec::new(); + self.for_each_neighbor_candidate_in_radius(position, radius, |entity| out.push(entity)); + out } /// Get the cell size in meters. - /// - /// **Phase A2**: Exposed for CFL stability checking in thermal system. pub fn cell_size(&self) -> f32 { - self.grid.cell_size + self.config.cell_size_meters + } + + fn rebuild_grid(&mut self) { + self.grid.clear(); + for (entity, position) in &self.entity_positions { + self.grid.insert(*entity, *position); + } + } + + fn rebuild_tree_if_needed(&mut self) { + if matches!(self.backend, BackendStorage::Tree) && self.tree_dirty { + self.tree.rebuild_from_positions(&self.entity_positions); + self.tree_dirty = false; + } + } + + fn switch_backend(&mut self, backend: BackendStorage) { + self.backend = backend; + match self.backend { + BackendStorage::Grid => self.rebuild_grid(), + BackendStorage::Tree => { + self.tree_dirty = true; + self.rebuild_tree_if_needed(); + } + } + self.frames_since_switch = 0; + } + + fn estimate_entities_per_cell(&self) -> f32 { + let n = self.entity_positions.len(); + if n == 0 { + return 0.0; + } + + let mut min = Vec2::splat(f32::MAX); + let mut max = Vec2::splat(f32::MIN); + for position in self.entity_positions.values() { + min.x = min.x.min(position.x); + min.y = min.y.min(position.y); + max.x = max.x.max(position.x); + max.y = max.y.max(position.y); + } + + let dx = (max.x - min.x).max(self.config.cell_size_meters); + let dy = (max.y - min.y).max(self.config.cell_size_meters); + let cells_x = (dx / self.config.cell_size_meters).ceil().max(1.0); + let cells_y = (dy / self.config.cell_size_meters).ceil().max(1.0); + let estimated_cells = (cells_x * cells_y).max(1.0); + n as f32 / estimated_cells + } + + fn preferred_backend_for_auto(&self) -> BackendStorage { + let n = self.entity_positions.len(); + if n < self.config.adaptive_min_entities_for_hierarchy { + return BackendStorage::Grid; + } + let entities_per_cell = self.estimate_entities_per_cell(); + if entities_per_cell <= self.config.adaptive_sparse_entities_per_cell_threshold { + BackendStorage::Tree + } else { + BackendStorage::Grid + } + } + + /// Estimated local packing for adaptive broadphase policy. + pub fn estimated_entities_per_cell(&self) -> f32 { + self.estimate_entities_per_cell() + } + + /// Sync backend mode and rebuild backend storage if needed. + pub fn prepare_for_queries(&mut self) { + self.frames_since_switch = self.frames_since_switch.saturating_add(1); + + let target_backend = match self.config.mode { + NeighborSearchMode::UniformCellField => BackendStorage::Grid, + NeighborSearchMode::HierarchicalEnvelopeField => BackendStorage::Tree, + NeighborSearchMode::Adaptive => { + if self.frames_since_switch < self.config.adaptive_switch_cooldown_frames { + self.backend + } else { + self.preferred_backend_for_auto() + } + } + }; + + if target_backend != self.backend { + self.switch_backend(target_backend); + } + + self.rebuild_tree_if_needed(); } } impl Default for UnifiedSpatialIndex { fn default() -> Self { - // Explicit default, performance parameter - Self::new(50.0) + Self::from_config(NeighborSearchConfig::default()) } } @@ -123,7 +500,7 @@ pub fn attach_spatial_cells( ) { for (e, t) in q.iter() { let pos = t.translation.truncate(); - let cell = index.insert(e, pos); // Returns cell + let cell = index.insert(e, pos); commands.entity(e).insert(SpatialCell { cell }); } } @@ -138,7 +515,7 @@ pub fn update_spatial_index( ) { for (e, t, mut cell) in q.iter_mut() { let pos = t.translation.truncate(); - let new_cell = index.update(e, pos); // Returns cell + let new_cell = index.update(e, pos); cell.cell = new_cell; } } @@ -152,3 +529,108 @@ pub fn remove_from_index_on_marker_removed( index.remove(e); } } + +/// Sync runtime config and prepare active backend for query systems. +pub fn refresh_spatial_index_policy( + mut index: ResMut, + config: Res, +) { + index.set_config(&config); + index.prepare_for_queries(); +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn grid_mode_radius_query_returns_candidates() { + let mut world = World::new(); + let mut index = UnifiedSpatialIndex::default(); + index.set_config(&NeighborSearchConfig { + mode: NeighborSearchMode::UniformCellField, + ..Default::default() + }); + index.prepare_for_queries(); + + let a = world.spawn_empty().id(); + let b = world.spawn_empty().id(); + index.insert(a, Vec2::new(0.0, 0.0)); + index.insert(b, Vec2::new(500.0, 500.0)); + index.prepare_for_queries(); + + let results = index.query_radius(Vec2::new(0.0, 0.0), 25.0); + assert!(results.contains(&a)); + assert!(!results.contains(&b)); + } + + #[test] + fn tree_mode_radius_query_returns_candidates() { + let mut world = World::new(); + let mut index = UnifiedSpatialIndex::default(); + index.set_config(&NeighborSearchConfig { + mode: NeighborSearchMode::HierarchicalEnvelopeField, + ..Default::default() + }); + index.prepare_for_queries(); + + let a = world.spawn_empty().id(); + let b = world.spawn_empty().id(); + index.insert(a, Vec2::new(0.0, 0.0)); + index.insert(b, Vec2::new(500.0, 500.0)); + index.prepare_for_queries(); + + let results = index.query_radius(Vec2::new(0.0, 0.0), 25.0); + assert!(results.contains(&a)); + assert!(!results.contains(&b)); + } + + #[test] + fn adaptive_prefers_hierarchy_for_sparse_large_sets() { + let mut world = World::new(); + let mut index = UnifiedSpatialIndex::default(); + index.set_config(&NeighborSearchConfig { + mode: NeighborSearchMode::Adaptive, + cell_size_meters: 10.0, + adaptive_min_entities_for_hierarchy: 64, + adaptive_sparse_entities_per_cell_threshold: 0.20, + adaptive_switch_cooldown_frames: 1, + ..Default::default() + }); + + for i in 0..128 { + let e = world.spawn_empty().id(); + index.insert(e, Vec2::new(i as f32 * 100.0, 0.0)); + } + + index.prepare_for_queries(); + assert_eq!( + index.active_mode(), + NeighborSearchMode::HierarchicalEnvelopeField + ); + } + + #[test] + fn adaptive_prefers_uniform_cells_for_dense_sets() { + let mut world = World::new(); + let mut index = UnifiedSpatialIndex::default(); + index.set_config(&NeighborSearchConfig { + mode: NeighborSearchMode::Adaptive, + cell_size_meters: 10.0, + adaptive_min_entities_for_hierarchy: 64, + adaptive_sparse_entities_per_cell_threshold: 0.20, + adaptive_switch_cooldown_frames: 1, + ..Default::default() + }); + + for i in 0..128 { + let x = (i % 16) as f32; + let y = (i / 16) as f32; + let e = world.spawn_empty().id(); + index.insert(e, Vec2::new(x, y)); + } + + index.prepare_for_queries(); + assert_eq!(index.active_mode(), NeighborSearchMode::UniformCellField); + } +} From 98190ad773a2de640e35a906029654d51921e9d1 Mon Sep 17 00:00:00 2001 From: M1thieu <18742831+M1thieu@users.noreply.github.com> Date: Sun, 15 Feb 2026 03:24:12 +0100 Subject: [PATCH 2/7] refactor(core): unify pairwise staging and add safety scaffolds - energy: shared pairwise helpers used by charges/thermal\n- energy: reusable Local staging contexts + O(N) thermal sanity checks\n- forces: SoA staging path for one-way gravity hot loop\n- save_system: atomic writes + checkpoint/rollback world APIs\n- information: baseline entropy/MI estimator traits + deterministic tests\n\nPhysics behavior unchanged; focused on determinism/perf hygiene and observability. feat(energy): add optional strict pairwise neighbor determinism - add PairwiseDeterminismConfig (default off)\n- add shared ordered-neighbor emission helper\n- wire charges/thermal pairwise loops to strict mode\n- keep realtime fast path when strict mode is disabled\n\nThis closes inner-neighbor ordering nondeterminism for replay/campaign runs. --- crates/energy/src/electromagnetism/charges.rs | 121 +++++--- crates/energy/src/lib.rs | 5 + crates/energy/src/pairwise.rs | 68 +++++ crates/energy/src/thermodynamics/thermal.rs | 275 +++++++++++++----- crates/forces/src/core/gravity.rs | 40 ++- crates/information/src/measures/estimators.rs | 71 +++++ crates/information/src/measures/mod.rs | 5 + crates/systems/save_system/src/save_system.rs | 162 ++++++++++- 8 files changed, 616 insertions(+), 131 deletions(-) create mode 100644 crates/energy/src/pairwise.rs create mode 100644 crates/information/src/measures/estimators.rs diff --git a/crates/energy/src/electromagnetism/charges.rs b/crates/energy/src/electromagnetism/charges.rs index 25d3f3c..f28c2ab 100644 --- a/crates/energy/src/electromagnetism/charges.rs +++ b/crates/energy/src/electromagnetism/charges.rs @@ -12,6 +12,11 @@ use forces::core::newton_laws::AppliedForce; use std::collections::HashMap; use utils::{SpatiallyIndexed, UnifiedSpatialIndex, force_switch}; +use crate::pairwise::{ + PairwiseDeterminismConfig, for_each_neighbor_candidate, is_forward_entity_pair, + prepare_sorted_entities_from_keys, prepare_staging_map, +}; + /// Electric charge component. /// /// Units: Coulombs (C) @@ -58,6 +63,13 @@ pub struct SofteningLength { pub value: f32, } +#[derive(Default)] +pub(crate) struct CoulombComputeContext { + charge_data: HashMap, + sorted_entities: Vec, + neighbor_candidates: Vec, +} + impl Default for SofteningLength { fn default() -> Self { Self { value: 0.01 } @@ -133,7 +145,7 @@ pub fn mark_charged_entities_spatially_indexed( /// /// **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 +155,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,61 +186,71 @@ 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)); } // Deterministic outer iteration: sort entities by stable id - let mut sorted_entities: Vec = charge_data.keys().copied().collect(); - sorted_entities.sort_by_key(|e| e.to_bits()); + let staged_entities: Vec = ctx.charge_data.keys().copied().collect(); + prepare_sorted_entities_from_keys(&mut ctx.sorted_entities, staged_entities); // Iterate pairs via UnifiedSpatialIndex. - // TODO(LP-0 determinism mode): For exact replay builds, collect and sort emitted neighbors - // by entity id before accumulation. Default path stays allocation-free for realtime. - for entity_a in sorted_entities { + 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 - index.for_each_neighbor_candidate_in_radius(pos_a, config.cutoff_radius, |entity_b| { - // **Pair-once guarantee**: Only process pairs where B > A - if entity_b.to_bits() <= entity_a.to_bits() { - return; - } + // 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; + } - // Get data for entity B from staged map - let Some((charge_b, pos_b, soft_b)) = charge_data.get(&entity_b) else { - return; - }; + // Get data for entity B from staged map + let Some((charge_b, pos_b, soft_b)) = charge_data.get(&entity_b) else { + return; + }; - let r_vec = *pos_b - pos_a; - let r = r_vec.length(); + let r_vec = *pos_b - pos_a; + let r = r_vec.length(); - // 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 { - 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 { + return; + } - // 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 - } + // 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 + } - // **LP-0**: EM potential energy = 0 (force-only). - // Future: Track U(r) = integral of switched force for energy conservation. - }); + // **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/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 77c24d6..423e8de 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. @@ -222,8 +262,10 @@ pub(crate) fn compute_fourier_conduction( mut commands: Commands, index: Res, config: Res, + determinism: Res, time: Res