From 18caa3ce07af2ab9255bf644e992911b4dd8ba87 Mon Sep 17 00:00:00 2001 From: elizabeth-suehr Date: Mon, 16 Mar 2026 16:46:13 -0700 Subject: [PATCH 1/2] Add rotational viscous damping fix: torque = -gamma * omega Implements angular velocity-proportional damping in mddem_fixes, mirroring the existing translational viscous fix. Registered at PostForce, accesses DemAtom via AtomDataRegistry for omega and torque fields. Includes 4 tests covering torque opposition, zero-at-rest, group filtering, and linear scaling. Co-Authored-By: Claude Opus 4.6 --- Cargo.lock | 1 + crates/mddem/src/lib.rs | 2 +- crates/mddem_fixes/Cargo.toml | 1 + crates/mddem_fixes/src/lib.rs | 192 +++++++++++++++++++++++++++++++++- 4 files changed, 192 insertions(+), 4 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 1e9e865..255a83e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -530,6 +530,7 @@ dependencies = [ name = "mddem_fixes" version = "0.1.3" dependencies = [ + "dem_atom", "mddem_app", "mddem_core", "mddem_scheduler", diff --git a/crates/mddem/src/lib.rs b/crates/mddem/src/lib.rs index 90bc6b8..61f78eb 100644 --- a/crates/mddem/src/lib.rs +++ b/crates/mddem/src/lib.rs @@ -112,7 +112,7 @@ pub mod prelude { pub use md_polymer::{ChainStatsConfig, ChainStatsPlugin, PolymerConfig, PolymerInitConfig, PolymerInitPlugin, PolymerPlugin}; pub use md_thermostat::{LangevinConfig, LangevinPlugin, LangevinState, NoseHooverPlugin, NoseHooverState, ThermostatConfig}; pub use mddem_deform::{DeformConfig, DeformPlugin, DeformState}; - pub use mddem_fixes::{AddForceDef, FixesPlugin, FixesRegistry, FreezeDef, MoveLinearDef, SetForceDef, ViscousDef}; + pub use mddem_fixes::{AddForceDef, FixesPlugin, FixesRegistry, FreezeDef, MoveLinearDef, RotationalViscousDef, SetForceDef, ViscousDef}; pub use mddem_app::prelude::*; pub use mddem_derive::StageEnum; pub use mddem_core::*; diff --git a/crates/mddem_fixes/Cargo.toml b/crates/mddem_fixes/Cargo.toml index 99993da..b7084b2 100644 --- a/crates/mddem_fixes/Cargo.toml +++ b/crates/mddem_fixes/Cargo.toml @@ -12,6 +12,7 @@ categories = ["simulation", "science"] mddem_core = { path = "../mddem_core", version = "0.1.3", default-features = false } mddem_app = { path = "../mddem_app", version = "0.1.3" } mddem_scheduler = { path = "../mddem_scheduler", version = "0.1.3" } +dem_atom = { path = "../dem_atom", version = "0.1.3" } serde = { version = "1", features = ["derive"] } [dev-dependencies] diff --git a/crates/mddem_fixes/src/lib.rs b/crates/mddem_fixes/src/lib.rs index 2cd3163..51e5d77 100644 --- a/crates/mddem_fixes/src/lib.rs +++ b/crates/mddem_fixes/src/lib.rs @@ -4,7 +4,7 @@ use mddem_app::prelude::*; use mddem_scheduler::prelude::*; use serde::Deserialize; -use mddem_core::{Atom, CommResource, Config, GroupRegistry}; +use mddem_core::{Atom, AtomDataRegistry, CommResource, Config, GroupRegistry}; // ── Config structs ───────────────────────────────────────────────────────── @@ -61,6 +61,13 @@ pub struct ViscousDef { pub gamma: f64, } +#[derive(Deserialize, Clone, Debug)] +#[serde(deny_unknown_fields)] +pub struct RotationalViscousDef { + pub group: String, + pub gamma: f64, +} + // ── Registry ─────────────────────────────────────────────────────────────── pub struct FixesRegistry { @@ -69,6 +76,7 @@ pub struct FixesRegistry { pub move_linears: Vec, pub freezes: Vec, pub viscous: Vec, + pub rotational_viscous: Vec, } // ── Plugin ───────────────────────────────────────────────────────────────── @@ -97,7 +105,15 @@ impl Plugin for FixesPlugin { # vz = -0.001 # [[freeze]] -# group = "frozen""#, +# group = "frozen" + +# [[viscous]] +# group = "all" +# gamma = 0.1 # translational damping coefficient (force = -gamma * velocity) + +# [[rotational_viscous]] +# group = "all" +# gamma = 0.001 # angular damping coefficient (torque = -gamma * omega)"#, ) } @@ -112,6 +128,7 @@ impl Plugin for FixesPlugin { move_linears: config.parse_array::("move_linear"), freezes: config.parse_array::("freeze"), viscous: config.parse_array::("viscous"), + rotational_viscous: config.parse_array::("rotational_viscous"), }; drop(config); @@ -120,7 +137,8 @@ impl Plugin for FixesPlugin { || !registry.set_forces.is_empty() || !registry.move_linears.is_empty() || !registry.freezes.is_empty() - || !registry.viscous.is_empty(); + || !registry.viscous.is_empty() + || !registry.rotational_viscous.is_empty(); if !has_any { app.add_resource(registry); @@ -132,6 +150,7 @@ impl Plugin for FixesPlugin { let has_set = !registry.set_forces.is_empty(); let has_freeze = !registry.freezes.is_empty(); let has_viscous = !registry.viscous.is_empty(); + let has_rotational_viscous = !registry.rotational_viscous.is_empty(); app.add_resource(registry) .add_setup_system(setup_fixes, ScheduleSetupSet::PostSetup); @@ -152,6 +171,9 @@ impl Plugin for FixesPlugin { if has_viscous { app.add_update_system(apply_viscous, ScheduleSet::PostForce); } + if has_rotational_viscous { + app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); + } } } @@ -174,6 +196,9 @@ fn setup_fixes(registry: Res, comm: Res, groups: Re for f in ®istry.viscous { groups.validate_name(&f.group, "fix viscous"); } + for f in ®istry.rotational_viscous { + groups.validate_name(&f.group, "fix rotational_viscous"); + } if comm.rank() != 0 { return; @@ -202,6 +227,12 @@ fn setup_fixes(registry: Res, comm: Res, groups: Re for f in ®istry.viscous { println!("Fix viscous: group='{}', gamma={}", f.group, f.gamma); } + for f in ®istry.rotational_viscous { + println!( + "Fix rotational_viscous: group='{}', gamma={}", + f.group, f.gamma + ); + } } /// Set velocity = constant BEFORE Verlet position update, so positions move at prescribed rate. @@ -322,6 +353,28 @@ fn apply_viscous( } } +/// Apply angular-velocity-proportional damping: τ = -gamma * ω. +fn apply_rotational_viscous( + atoms: Res, + atom_data: Res, + registry: Res, + groups: Res, +) { + let mut dem = atom_data.expect_mut::("apply_rotational_viscous"); + let nlocal = atoms.nlocal as usize; + for def in ®istry.rotational_viscous { + let group = groups.expect(&def.group); + let gamma = def.gamma; + for i in 0..nlocal { + if group.mask[i] { + dem.torque[i][0] -= gamma * dem.omega[i][0]; + dem.torque[i][1] -= gamma * dem.omega[i][1]; + dem.torque[i][2] -= gamma * dem.omega[i][2]; + } + } + } +} + // ── Gravity ───────────────────────────────────────────────────────────────── #[derive(Deserialize, Clone)] @@ -403,6 +456,7 @@ mod tests { move_linears: vec![], freezes: vec![], viscous: vec![], + rotational_viscous: vec![], }; // Set some initial force @@ -444,6 +498,7 @@ mod tests { move_linears: vec![], freezes: vec![], viscous: vec![], + rotational_viscous: vec![], }; let mut app = App::new(); @@ -479,6 +534,7 @@ mod tests { group: "frozen".to_string(), }], viscous: vec![], + rotational_viscous: vec![], }; let mut app = App::new(); @@ -513,6 +569,7 @@ mod tests { }], freezes: vec![], viscous: vec![], + rotational_viscous: vec![], }; // Pre step: sets velocity @@ -551,6 +608,7 @@ mod tests { group: "all".to_string(), gamma: 0.1, }], + rotational_viscous: vec![], }; let mut app = App::new(); @@ -580,6 +638,7 @@ mod tests { group: "all".to_string(), gamma: 0.1, }], + rotational_viscous: vec![], }; let mut app = App::new(); @@ -596,6 +655,133 @@ mod tests { assert!((a.force[0][2]).abs() < 1e-15); } + // ── Rotational viscous tests ─────────────────────────────────────────── + + fn make_rotational_viscous_app( + n: usize, + group_name: &str, + mask: Vec, + gamma: f64, + omegas: Vec<[f64; 3]>, + ) -> App { + let atoms = make_atoms(n); + let groups = make_group_registry(group_name, mask); + + let mut dem = dem_atom::DemAtom::new(); + for omega in &omegas { + dem.radius.push(0.5); + dem.density.push(2500.0); + dem.inv_inertia.push(1.0); + dem.quaternion.push([1.0, 0.0, 0.0, 0.0]); + dem.omega.push(*omega); + dem.ang_mom.push([0.0; 3]); + dem.torque.push([0.0; 3]); + } + + let mut registry = AtomDataRegistry::new(); + registry.register(dem); + + let fixes = FixesRegistry { + add_forces: vec![], + set_forces: vec![], + move_linears: vec![], + freezes: vec![], + viscous: vec![], + rotational_viscous: vec![RotationalViscousDef { + group: group_name.to_string(), + gamma, + }], + }; + + let mut app = App::new(); + app.add_resource(atoms); + app.add_resource(groups); + app.add_resource(registry); + app.add_resource(fixes); + app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); + app.organize_systems(); + app.run(); + app + } + + #[test] + fn test_rotational_viscous_opposes_angular_velocity() { + let app = make_rotational_viscous_app( + 2, + "all", + vec![true, true], + 0.1, + vec![[1.0, -2.0, 0.5], [0.0; 3]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + // torque = -gamma * omega + assert!((dem.torque[0][0] - (-0.1)).abs() < 1e-12, "tx = -gamma*wx"); + assert!((dem.torque[0][1] - 0.2).abs() < 1e-12, "ty = -gamma*wy"); + assert!( + (dem.torque[0][2] - (-0.05)).abs() < 1e-12, + "tz = -gamma*wz" + ); + } + + #[test] + fn test_rotational_viscous_zero_at_rest() { + let app = make_rotational_viscous_app( + 2, + "all", + vec![true, true], + 0.1, + vec![[0.0; 3], [0.0; 3]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!((dem.torque[0][0]).abs() < 1e-15); + assert!((dem.torque[0][1]).abs() < 1e-15); + assert!((dem.torque[0][2]).abs() < 1e-15); + } + + #[test] + fn test_rotational_viscous_group_filtering() { + let app = make_rotational_viscous_app( + 3, + "subset", + vec![true, false, true], + 0.5, + vec![[2.0, 0.0, 0.0], [3.0, 0.0, 0.0], [4.0, 0.0, 0.0]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + // Atom 0 in group: torque = -0.5 * 2.0 = -1.0 + assert!((dem.torque[0][0] - (-1.0)).abs() < 1e-12); + // Atom 1 NOT in group: torque unchanged (0) + assert!((dem.torque[1][0]).abs() < 1e-15); + // Atom 2 in group: torque = -0.5 * 4.0 = -2.0 + assert!((dem.torque[2][0] - (-2.0)).abs() < 1e-12); + } + + #[test] + fn test_rotational_viscous_scales_with_gamma_and_omega() { + // gamma=0.2, omega=[3.0, 0.0, 0.0] -> torque = -0.6 + let app = make_rotational_viscous_app( + 1, + "all", + vec![true], + 0.2, + vec![[3.0, 0.0, 0.0]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!( + (dem.torque[0][0] - (-0.6)).abs() < 1e-12, + "torque should scale linearly: -0.2*3.0 = -0.6, got {}", + dem.torque[0][0] + ); + } + // ── Gravity tests ────────────────────────────────────────────────────── fn make_gravity_atom(mass: f64) -> Atom { From 2f4943790b158312e4d4628cadef1a7764e7e4ae Mon Sep 17 00:00:00 2001 From: elizabeth-suehr Date: Mon, 16 Mar 2026 17:02:41 -0700 Subject: [PATCH 2/2] Move rotational viscous damping from mddem_fixes to dem_granular Fixes architectural layering: mddem_* (shared infra) should not depend on dem_* (DEM-specific) crates. Rotational viscous damping needs DemAtom for omega/torque access, so it belongs in dem_granular. - Remove dem_atom dependency from mddem_fixes - Add RotationalViscousPlugin to dem_granular::rotational_viscous - Re-export from mddem prelude - All 4 rotational viscous tests preserved and passing Co-Authored-By: Claude Opus 4.6 --- Cargo.lock | 2 +- crates/dem_granular/Cargo.toml | 1 + crates/dem_granular/src/lib.rs | 2 + crates/dem_granular/src/rotational_viscous.rs | 212 ++++++++++++++++++ crates/mddem/src/lib.rs | 3 +- crates/mddem_fixes/Cargo.toml | 1 - crates/mddem_fixes/src/lib.rs | 188 +--------------- 7 files changed, 221 insertions(+), 188 deletions(-) create mode 100644 crates/dem_granular/src/rotational_viscous.rs diff --git a/Cargo.lock b/Cargo.lock index 255a83e..3feae99 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -166,6 +166,7 @@ dependencies = [ "mddem_scheduler", "mddem_test_utils", "mddem_verlet", + "serde", ] [[package]] @@ -530,7 +531,6 @@ dependencies = [ name = "mddem_fixes" version = "0.1.3" dependencies = [ - "dem_atom", "mddem_app", "mddem_core", "mddem_scheduler", diff --git a/crates/dem_granular/Cargo.toml b/crates/dem_granular/Cargo.toml index 8199206..3b9f39a 100644 --- a/crates/dem_granular/Cargo.toml +++ b/crates/dem_granular/Cargo.toml @@ -16,6 +16,7 @@ mddem_app = { path = "../mddem_app", version = "0.1.3" } mddem_scheduler = { path = "../mddem_scheduler", version = "0.1.3" } dem_atom = { path = "../dem_atom", version = "0.1.3" } mddem_verlet = { path = "../mddem_verlet", version = "0.1.3" } +serde = { version = "1", features = ["derive"] } [dev-dependencies] mddem_test_utils = { path = "../mddem_test_utils", version = "0.1.3" } diff --git a/crates/dem_granular/src/lib.rs b/crates/dem_granular/src/lib.rs index 78dda3f..83c367f 100644 --- a/crates/dem_granular/src/lib.rs +++ b/crates/dem_granular/src/lib.rs @@ -3,11 +3,13 @@ pub mod granular_temp; pub mod normal; pub mod rotational; +pub mod rotational_viscous; pub mod tangential; pub use granular_temp::GranularTempPlugin; pub use normal::HertzNormalForcePlugin; pub use rotational::RotationalDynamicsPlugin; +pub use rotational_viscous::RotationalViscousPlugin; pub use tangential::MindlinTangentialForcePlugin; pub mod contact; diff --git a/crates/dem_granular/src/rotational_viscous.rs b/crates/dem_granular/src/rotational_viscous.rs new file mode 100644 index 0000000..985c282 --- /dev/null +++ b/crates/dem_granular/src/rotational_viscous.rs @@ -0,0 +1,212 @@ +//! Rotational viscous damping: torque = -gamma * omega. + +use mddem_app::prelude::*; +use mddem_scheduler::prelude::*; +use serde::Deserialize; + +use dem_atom::DemAtom; +use mddem_core::{Atom, AtomDataRegistry, CommResource, Config, GroupRegistry}; + +#[derive(Deserialize, Clone, Debug)] +#[serde(deny_unknown_fields)] +pub struct RotationalViscousDef { + pub group: String, + pub gamma: f64, +} + +pub struct RotationalViscousRegistry { + pub defs: Vec, +} + +pub struct RotationalViscousPlugin; + +impl Plugin for RotationalViscousPlugin { + fn default_config(&self) -> Option<&str> { + Some( + r#"# [[rotational_viscous]] +# group = "all" +# gamma = 0.001 # angular damping coefficient (torque = -gamma * omega)"#, + ) + } + + fn build(&self, app: &mut App) { + let config = app + .get_resource_ref::() + .expect("Config resource must exist before RotationalViscousPlugin"); + + let defs = config.parse_array::("rotational_viscous"); + let has_any = !defs.is_empty(); + drop(config); + + let registry = RotationalViscousRegistry { defs }; + app.add_resource(registry); + + if has_any { + app.add_setup_system(setup_rotational_viscous, ScheduleSetupSet::PostSetup); + app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); + } + } +} + +fn setup_rotational_viscous( + registry: Res, + comm: Res, + groups: Res, +) { + for f in ®istry.defs { + groups.validate_name(&f.group, "fix rotational_viscous"); + } + + if comm.rank() != 0 { + return; + } + for f in ®istry.defs { + println!( + "Fix rotational_viscous: group='{}', gamma={}", + f.group, f.gamma + ); + } +} + +/// Apply angular-velocity-proportional damping: τ = -gamma * ω. +fn apply_rotational_viscous( + atoms: Res, + atom_data: Res, + registry: Res, + groups: Res, +) { + let mut dem = atom_data.expect_mut::("apply_rotational_viscous"); + let nlocal = atoms.nlocal as usize; + for def in ®istry.defs { + let group = groups.expect(&def.group); + let gamma = def.gamma; + for i in 0..nlocal { + if group.mask[i] { + dem.torque[i][0] -= gamma * dem.omega[i][0]; + dem.torque[i][1] -= gamma * dem.omega[i][1]; + dem.torque[i][2] -= gamma * dem.omega[i][2]; + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use mddem_test_utils::{make_atoms, make_group_registry}; + + fn make_rotational_viscous_app( + n: usize, + group_name: &str, + mask: Vec, + gamma: f64, + omegas: Vec<[f64; 3]>, + ) -> App { + let atoms = make_atoms(n); + let groups = make_group_registry(group_name, mask); + + let mut dem = DemAtom::new(); + for omega in &omegas { + dem.radius.push(0.5); + dem.density.push(2500.0); + dem.inv_inertia.push(1.0); + dem.quaternion.push([1.0, 0.0, 0.0, 0.0]); + dem.omega.push(*omega); + dem.ang_mom.push([0.0; 3]); + dem.torque.push([0.0; 3]); + } + + let mut atom_registry = AtomDataRegistry::new(); + atom_registry.register(dem); + + let registry = RotationalViscousRegistry { + defs: vec![RotationalViscousDef { + group: group_name.to_string(), + gamma, + }], + }; + + let mut app = App::new(); + app.add_resource(atoms); + app.add_resource(groups); + app.add_resource(atom_registry); + app.add_resource(registry); + app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); + app.organize_systems(); + app.run(); + app + } + + #[test] + fn test_opposes_angular_velocity() { + let app = make_rotational_viscous_app( + 2, + "all", + vec![true, true], + 0.1, + vec![[1.0, -2.0, 0.5], [0.0; 3]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!((dem.torque[0][0] - (-0.1)).abs() < 1e-12, "tx = -gamma*wx"); + assert!((dem.torque[0][1] - 0.2).abs() < 1e-12, "ty = -gamma*wy"); + assert!( + (dem.torque[0][2] - (-0.05)).abs() < 1e-12, + "tz = -gamma*wz" + ); + } + + #[test] + fn test_zero_at_rest() { + let app = make_rotational_viscous_app( + 2, + "all", + vec![true, true], + 0.1, + vec![[0.0; 3], [0.0; 3]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!((dem.torque[0][0]).abs() < 1e-15); + assert!((dem.torque[0][1]).abs() < 1e-15); + assert!((dem.torque[0][2]).abs() < 1e-15); + } + + #[test] + fn test_group_filtering() { + let app = make_rotational_viscous_app( + 3, + "subset", + vec![true, false, true], + 0.5, + vec![[2.0, 0.0, 0.0], [3.0, 0.0, 0.0], [4.0, 0.0, 0.0]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!((dem.torque[0][0] - (-1.0)).abs() < 1e-12); + assert!((dem.torque[1][0]).abs() < 1e-15); + assert!((dem.torque[2][0] - (-2.0)).abs() < 1e-12); + } + + #[test] + fn test_scales_with_gamma_and_omega() { + let app = make_rotational_viscous_app( + 1, + "all", + vec![true], + 0.2, + vec![[3.0, 0.0, 0.0]], + ); + + let registry = app.get_resource_ref::().unwrap(); + let dem = registry.expect::("test"); + assert!( + (dem.torque[0][0] - (-0.6)).abs() < 1e-12, + "torque should scale linearly: -0.2*3.0 = -0.6, got {}", + dem.torque[0][0] + ); + } +} diff --git a/crates/mddem/src/lib.rs b/crates/mddem/src/lib.rs index 61f78eb..b36abd0 100644 --- a/crates/mddem/src/lib.rs +++ b/crates/mddem/src/lib.rs @@ -112,7 +112,8 @@ pub mod prelude { pub use md_polymer::{ChainStatsConfig, ChainStatsPlugin, PolymerConfig, PolymerInitConfig, PolymerInitPlugin, PolymerPlugin}; pub use md_thermostat::{LangevinConfig, LangevinPlugin, LangevinState, NoseHooverPlugin, NoseHooverState, ThermostatConfig}; pub use mddem_deform::{DeformConfig, DeformPlugin, DeformState}; - pub use mddem_fixes::{AddForceDef, FixesPlugin, FixesRegistry, FreezeDef, MoveLinearDef, RotationalViscousDef, SetForceDef, ViscousDef}; + pub use mddem_fixes::{AddForceDef, FixesPlugin, FixesRegistry, FreezeDef, MoveLinearDef, SetForceDef, ViscousDef}; + pub use dem_granular::rotational_viscous::{RotationalViscousDef, RotationalViscousPlugin, RotationalViscousRegistry}; pub use mddem_app::prelude::*; pub use mddem_derive::StageEnum; pub use mddem_core::*; diff --git a/crates/mddem_fixes/Cargo.toml b/crates/mddem_fixes/Cargo.toml index b7084b2..99993da 100644 --- a/crates/mddem_fixes/Cargo.toml +++ b/crates/mddem_fixes/Cargo.toml @@ -12,7 +12,6 @@ categories = ["simulation", "science"] mddem_core = { path = "../mddem_core", version = "0.1.3", default-features = false } mddem_app = { path = "../mddem_app", version = "0.1.3" } mddem_scheduler = { path = "../mddem_scheduler", version = "0.1.3" } -dem_atom = { path = "../dem_atom", version = "0.1.3" } serde = { version = "1", features = ["derive"] } [dev-dependencies] diff --git a/crates/mddem_fixes/src/lib.rs b/crates/mddem_fixes/src/lib.rs index 51e5d77..e7939bc 100644 --- a/crates/mddem_fixes/src/lib.rs +++ b/crates/mddem_fixes/src/lib.rs @@ -4,7 +4,7 @@ use mddem_app::prelude::*; use mddem_scheduler::prelude::*; use serde::Deserialize; -use mddem_core::{Atom, AtomDataRegistry, CommResource, Config, GroupRegistry}; +use mddem_core::{Atom, CommResource, Config, GroupRegistry}; // ── Config structs ───────────────────────────────────────────────────────── @@ -61,13 +61,6 @@ pub struct ViscousDef { pub gamma: f64, } -#[derive(Deserialize, Clone, Debug)] -#[serde(deny_unknown_fields)] -pub struct RotationalViscousDef { - pub group: String, - pub gamma: f64, -} - // ── Registry ─────────────────────────────────────────────────────────────── pub struct FixesRegistry { @@ -76,7 +69,6 @@ pub struct FixesRegistry { pub move_linears: Vec, pub freezes: Vec, pub viscous: Vec, - pub rotational_viscous: Vec, } // ── Plugin ───────────────────────────────────────────────────────────────── @@ -110,10 +102,7 @@ impl Plugin for FixesPlugin { # [[viscous]] # group = "all" # gamma = 0.1 # translational damping coefficient (force = -gamma * velocity) - -# [[rotational_viscous]] -# group = "all" -# gamma = 0.001 # angular damping coefficient (torque = -gamma * omega)"#, +"#, ) } @@ -128,7 +117,6 @@ impl Plugin for FixesPlugin { move_linears: config.parse_array::("move_linear"), freezes: config.parse_array::("freeze"), viscous: config.parse_array::("viscous"), - rotational_viscous: config.parse_array::("rotational_viscous"), }; drop(config); @@ -137,8 +125,7 @@ impl Plugin for FixesPlugin { || !registry.set_forces.is_empty() || !registry.move_linears.is_empty() || !registry.freezes.is_empty() - || !registry.viscous.is_empty() - || !registry.rotational_viscous.is_empty(); + || !registry.viscous.is_empty(); if !has_any { app.add_resource(registry); @@ -150,7 +137,6 @@ impl Plugin for FixesPlugin { let has_set = !registry.set_forces.is_empty(); let has_freeze = !registry.freezes.is_empty(); let has_viscous = !registry.viscous.is_empty(); - let has_rotational_viscous = !registry.rotational_viscous.is_empty(); app.add_resource(registry) .add_setup_system(setup_fixes, ScheduleSetupSet::PostSetup); @@ -171,9 +157,6 @@ impl Plugin for FixesPlugin { if has_viscous { app.add_update_system(apply_viscous, ScheduleSet::PostForce); } - if has_rotational_viscous { - app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); - } } } @@ -196,10 +179,6 @@ fn setup_fixes(registry: Res, comm: Res, groups: Re for f in ®istry.viscous { groups.validate_name(&f.group, "fix viscous"); } - for f in ®istry.rotational_viscous { - groups.validate_name(&f.group, "fix rotational_viscous"); - } - if comm.rank() != 0 { return; } @@ -227,12 +206,6 @@ fn setup_fixes(registry: Res, comm: Res, groups: Re for f in ®istry.viscous { println!("Fix viscous: group='{}', gamma={}", f.group, f.gamma); } - for f in ®istry.rotational_viscous { - println!( - "Fix rotational_viscous: group='{}', gamma={}", - f.group, f.gamma - ); - } } /// Set velocity = constant BEFORE Verlet position update, so positions move at prescribed rate. @@ -353,28 +326,6 @@ fn apply_viscous( } } -/// Apply angular-velocity-proportional damping: τ = -gamma * ω. -fn apply_rotational_viscous( - atoms: Res, - atom_data: Res, - registry: Res, - groups: Res, -) { - let mut dem = atom_data.expect_mut::("apply_rotational_viscous"); - let nlocal = atoms.nlocal as usize; - for def in ®istry.rotational_viscous { - let group = groups.expect(&def.group); - let gamma = def.gamma; - for i in 0..nlocal { - if group.mask[i] { - dem.torque[i][0] -= gamma * dem.omega[i][0]; - dem.torque[i][1] -= gamma * dem.omega[i][1]; - dem.torque[i][2] -= gamma * dem.omega[i][2]; - } - } - } -} - // ── Gravity ───────────────────────────────────────────────────────────────── #[derive(Deserialize, Clone)] @@ -456,7 +407,6 @@ mod tests { move_linears: vec![], freezes: vec![], viscous: vec![], - rotational_viscous: vec![], }; // Set some initial force @@ -498,7 +448,6 @@ mod tests { move_linears: vec![], freezes: vec![], viscous: vec![], - rotational_viscous: vec![], }; let mut app = App::new(); @@ -534,7 +483,6 @@ mod tests { group: "frozen".to_string(), }], viscous: vec![], - rotational_viscous: vec![], }; let mut app = App::new(); @@ -569,7 +517,6 @@ mod tests { }], freezes: vec![], viscous: vec![], - rotational_viscous: vec![], }; // Pre step: sets velocity @@ -608,7 +555,6 @@ mod tests { group: "all".to_string(), gamma: 0.1, }], - rotational_viscous: vec![], }; let mut app = App::new(); @@ -638,7 +584,6 @@ mod tests { group: "all".to_string(), gamma: 0.1, }], - rotational_viscous: vec![], }; let mut app = App::new(); @@ -655,133 +600,6 @@ mod tests { assert!((a.force[0][2]).abs() < 1e-15); } - // ── Rotational viscous tests ─────────────────────────────────────────── - - fn make_rotational_viscous_app( - n: usize, - group_name: &str, - mask: Vec, - gamma: f64, - omegas: Vec<[f64; 3]>, - ) -> App { - let atoms = make_atoms(n); - let groups = make_group_registry(group_name, mask); - - let mut dem = dem_atom::DemAtom::new(); - for omega in &omegas { - dem.radius.push(0.5); - dem.density.push(2500.0); - dem.inv_inertia.push(1.0); - dem.quaternion.push([1.0, 0.0, 0.0, 0.0]); - dem.omega.push(*omega); - dem.ang_mom.push([0.0; 3]); - dem.torque.push([0.0; 3]); - } - - let mut registry = AtomDataRegistry::new(); - registry.register(dem); - - let fixes = FixesRegistry { - add_forces: vec![], - set_forces: vec![], - move_linears: vec![], - freezes: vec![], - viscous: vec![], - rotational_viscous: vec![RotationalViscousDef { - group: group_name.to_string(), - gamma, - }], - }; - - let mut app = App::new(); - app.add_resource(atoms); - app.add_resource(groups); - app.add_resource(registry); - app.add_resource(fixes); - app.add_update_system(apply_rotational_viscous, ScheduleSet::PostForce); - app.organize_systems(); - app.run(); - app - } - - #[test] - fn test_rotational_viscous_opposes_angular_velocity() { - let app = make_rotational_viscous_app( - 2, - "all", - vec![true, true], - 0.1, - vec![[1.0, -2.0, 0.5], [0.0; 3]], - ); - - let registry = app.get_resource_ref::().unwrap(); - let dem = registry.expect::("test"); - // torque = -gamma * omega - assert!((dem.torque[0][0] - (-0.1)).abs() < 1e-12, "tx = -gamma*wx"); - assert!((dem.torque[0][1] - 0.2).abs() < 1e-12, "ty = -gamma*wy"); - assert!( - (dem.torque[0][2] - (-0.05)).abs() < 1e-12, - "tz = -gamma*wz" - ); - } - - #[test] - fn test_rotational_viscous_zero_at_rest() { - let app = make_rotational_viscous_app( - 2, - "all", - vec![true, true], - 0.1, - vec![[0.0; 3], [0.0; 3]], - ); - - let registry = app.get_resource_ref::().unwrap(); - let dem = registry.expect::("test"); - assert!((dem.torque[0][0]).abs() < 1e-15); - assert!((dem.torque[0][1]).abs() < 1e-15); - assert!((dem.torque[0][2]).abs() < 1e-15); - } - - #[test] - fn test_rotational_viscous_group_filtering() { - let app = make_rotational_viscous_app( - 3, - "subset", - vec![true, false, true], - 0.5, - vec![[2.0, 0.0, 0.0], [3.0, 0.0, 0.0], [4.0, 0.0, 0.0]], - ); - - let registry = app.get_resource_ref::().unwrap(); - let dem = registry.expect::("test"); - // Atom 0 in group: torque = -0.5 * 2.0 = -1.0 - assert!((dem.torque[0][0] - (-1.0)).abs() < 1e-12); - // Atom 1 NOT in group: torque unchanged (0) - assert!((dem.torque[1][0]).abs() < 1e-15); - // Atom 2 in group: torque = -0.5 * 4.0 = -2.0 - assert!((dem.torque[2][0] - (-2.0)).abs() < 1e-12); - } - - #[test] - fn test_rotational_viscous_scales_with_gamma_and_omega() { - // gamma=0.2, omega=[3.0, 0.0, 0.0] -> torque = -0.6 - let app = make_rotational_viscous_app( - 1, - "all", - vec![true], - 0.2, - vec![[3.0, 0.0, 0.0]], - ); - - let registry = app.get_resource_ref::().unwrap(); - let dem = registry.expect::("test"); - assert!( - (dem.torque[0][0] - (-0.6)).abs() < 1e-12, - "torque should scale linearly: -0.2*3.0 = -0.6, got {}", - dem.torque[0][0] - ); - } - // ── Gravity tests ────────────────────────────────────────────────────── fn make_gravity_atom(mass: f64) -> Atom {