diff --git a/Cargo.lock b/Cargo.lock index 1e9e865..3feae99 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -166,6 +166,7 @@ dependencies = [ "mddem_scheduler", "mddem_test_utils", "mddem_verlet", + "serde", ] [[package]] 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 90bc6b8..b36abd0 100644 --- a/crates/mddem/src/lib.rs +++ b/crates/mddem/src/lib.rs @@ -113,6 +113,7 @@ pub mod prelude { 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 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/src/lib.rs b/crates/mddem_fixes/src/lib.rs index 2cd3163..e7939bc 100644 --- a/crates/mddem_fixes/src/lib.rs +++ b/crates/mddem_fixes/src/lib.rs @@ -97,7 +97,12 @@ impl Plugin for FixesPlugin { # vz = -0.001 # [[freeze]] -# group = "frozen""#, +# group = "frozen" + +# [[viscous]] +# group = "all" +# gamma = 0.1 # translational damping coefficient (force = -gamma * velocity) +"#, ) } @@ -174,7 +179,6 @@ fn setup_fixes(registry: Res, comm: Res, groups: Re for f in ®istry.viscous { groups.validate_name(&f.group, "fix viscous"); } - if comm.rank() != 0 { return; }