From a2086a59decf4b76b06fc684e71d0c2da0ca759c Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Tue, 15 Sep 2026 08:30:55 -0500 Subject: [PATCH] =?UTF-8?q?P6-b:=20the=20fictitious=20added=20mass=20for?= =?UTF-8?q?=20the=20partitioned=20loop=20=E2=80=94=20rtx-fea=20NonlinearDy?= =?UTF-8?q?namicStepper::set=5Fadded=5Flumped=5Fmass=20(a=20lumped=20per-D?= =?UTF-8?q?OF=20mass=20in=20the=20Newmark=20inertial=20residual=20and=20ef?= =?UTF-8?q?fective=20tangent,=20never=20the=20consistent=20mass=20or=20the?= =?UTF-8?q?=20rest=20state;=20zero=20=3D=20the=20plain=20stepper=20bit=20f?= =?UTF-8?q?or=20bit)=20with=20its=20pin=20(compensated=20step=20reproduces?= =?UTF-8?q?=20the=20plain=20step=20to=202.5e-9,=20uncompensated=20moves=20?= =?UTF-8?q?it=2011=20%);=20the=20FSI2=20overset=20harness=20carries=20RTX?= =?UTF-8?q?=5FFSI2O=5FFICT=5FMASS=3D=CE=B1=20(=CE=B1=20=C3=97=20=CF=81=5Ff?= =?UTF-8?q?=20=CF=80=20(c/2)=C2=B2=20spread=20over=20the=20wetted=20nodes)?= =?UTF-8?q?=20and=20adds=20the=20compensating=20load=20M=5Ff=20=C3=BC=5Fk?= =?UTF-8?q?=20of=20the=20previous=20subiterate=20to=20every=20structure=20?= =?UTF-8?q?solve=20(predictor=20and=20passes),=20printed=20in=20the=20head?= =?UTF-8?q?er?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01YJPeT6WA2e7YvAnS875AHL --- .../rtx-fea/src/analysis/nonlinear_dynamic.rs | 39 +++- .../rtx-fea/tests/fictitious_mass.rs | 190 ++++++++++++++++++ .../tests/fsi2_harness/overset_march.rs | 65 +++++- 3 files changed, 289 insertions(+), 5 deletions(-) create mode 100644 crates/specialized/rtx-fea/tests/fictitious_mass.rs diff --git a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs index aa1b93c..7073086 100644 --- a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs +++ b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs @@ -46,13 +46,13 @@ //! the stepper. use super::{AnalysisConfig, ConvergenceCriteria}; -use crate::assembly::SparseMatrix; use crate::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy}; +use crate::assembly::SparseMatrix; use crate::boundary::{BoundaryCondition, BoundaryConditionSet}; use crate::elements::total_lagrangian::{self, saint_venant_kirchhoff}; use crate::elements::{ElementMatrixComputer, StandardFiniteElement}; use crate::error::{AnalysisError, FeaResult}; -use crate::materials::{MaterialDatabase, reduced_constitutive}; +use crate::materials::{reduced_constitutive, MaterialDatabase}; use crate::mesh::{Mesh, NodeId}; use crate::solvers::{BandedLu, LinearSolver, SolverOptions}; use nalgebra::{DMatrix, DVector, Vector3}; @@ -264,6 +264,11 @@ pub struct NonlinearDynamicStepper<'a> { caches: Vec, /// Free-free consistent mass, for consistent initial accelerations. mass_free: SparseMatrix, + /// A lumped mass added per DOF (global numbering) — the partitioned + /// coupling's fictitious added mass (`set_added_lumped_mass`): it enters + /// the Newmark inertial residual and the effective tangent, never the + /// consistent mass or the rest state. Zero by default. + added_mass: DVector, /// The body-force part of the external force (constant). external_body: DVector, /// Body force plus the current nodal forces. @@ -414,7 +419,9 @@ impl<'a> NonlinearDynamicStepper<'a> { } } + let added_mass = DVector::zeros(total_dofs); let mut stepper = Self { + added_mass, analysis, dof_numbering, free_dofs, @@ -447,6 +454,20 @@ impl<'a> NonlinearDynamicStepper<'a> { } } + /// Set the lumped mass added to every DOF of each node (the coupling + /// loop's fictitious added mass, `docs/overset_metal_campaign.md` §5.17 + /// in omni-cortex): `(M + M_f) ü = F + M_f ü_k` contracts at any mass + /// ratio and leaves the fixed point unchanged when the caller adds the + /// load `M_f ü_k` of the previous subiterate. Entries not listed keep + /// their value; `0.0` restores the plain stepper bit for bit. + pub fn set_added_lumped_mass(&mut self, entries: &[(NodeId, f64)]) { + for (node, m) in entries { + for dof in self.dof_numbering.get_node_dofs(*node) { + self.added_mass[dof] = *m; + } + } + } + /// The state at rest under the *current* external force: `u = v = 0`, /// the acceleration consistent with `M a0 = F_ext - f_int(0)`. pub fn rest_state(&mut self) -> FeaResult { @@ -775,6 +796,13 @@ impl<'a> NonlinearDynamicStepper<'a> { } } if with_tangent { + for (dof, &m) in self.added_mass.iter().enumerate() { + if m != 0.0 { + if let Some(free) = self.free_index[dof] { + tangent.add_entry(free, free, inv_beta_dt2 * m)?; + } + } + } tangent.finalize()?; } Ok((internal, tangent)) @@ -796,6 +824,13 @@ impl<'a> NonlinearDynamicStepper<'a> { } } } + for (dof, &m) in self.added_mass.iter().enumerate() { + if m != 0.0 { + if let Some(free) = self.free_index[dof] { + out[free] += m * a_full[dof]; + } + } + } out } } diff --git a/crates/specialized/rtx-fea/tests/fictitious_mass.rs b/crates/specialized/rtx-fea/tests/fictitious_mass.rs new file mode 100644 index 0000000..d0e0258 --- /dev/null +++ b/crates/specialized/rtx-fea/tests/fictitious_mass.rs @@ -0,0 +1,190 @@ +//! P6-b (`docs/overset_metal_campaign.md` §5.17): the fictitious added mass +//! for the partitioned coupling loop. The stepper carries a lumped mass +//! `M_f` on chosen nodes; the coupling loop adds the load `M_f ü_k` from +//! the previous subiterate. Pin: solving `(M + M_f) ü = F + M_f ü_A` with +//! `ü_A` the plain step's own acceleration reproduces the plain step (the +//! fixed point is unchanged — the added terms cancel at convergence), while +//! the added mass WITHOUT its compensating load changes the step (the hook +//! is live). The Newmark effective matrix carries `M_f` (a step from rest +//! under a nodal load has the acceleration of a heavier body). + +use nalgebra::Vector3; +use rtx_fea::analysis::{AnalysisConfig, NonlinearDynamicAnalysis}; +use rtx_fea::assembly::dof_mapping::DofComponent; +use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType}; +use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction}; +use rtx_fea::materials::{LinearElastic, MaterialDatabase}; +use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId}; + +const E_MOD: f64 = 1.4e6; +const NU: f64 = 0.4; +const RHO: f64 = 1000.0; + +fn quad8_rect_mesh(x0: f64, x1: f64, y0: f64, y1: f64, nx: usize, ny: usize) -> Mesh { + let mut mesh = Mesh::new(2).unwrap(); + let (lx, ly) = (2 * nx + 1, 2 * ny + 1); + let mut grid = vec![vec![None; ly]; lx]; + for (i, column) in grid.iter_mut().enumerate() { + for (j, slot) in column.iter_mut().enumerate() { + if i % 2 == 1 && j % 2 == 1 { + continue; + } + let x = x0 + (x1 - x0) * i as f64 / (2 * nx) as f64; + let y = y0 + (y1 - y0) * j as f64 / (2 * ny) as f64; + *slot = Some(mesh.add_node(Node::new_2d(x, y))); + } + } + for i in 0..nx { + for j in 0..ny { + let (a, b) = (2 * i, 2 * j); + let nodes = vec![ + grid[a][b].unwrap(), + grid[a + 2][b].unwrap(), + grid[a + 2][b + 2].unwrap(), + grid[a][b + 2].unwrap(), + grid[a + 1][b].unwrap(), + grid[a + 2][b + 1].unwrap(), + grid[a + 1][b + 2].unwrap(), + grid[a][b + 1].unwrap(), + ]; + mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap()) + .unwrap(); + } + } + mesh +} + +fn materials() -> MaterialDatabase { + let mut db = MaterialDatabase::new(); + db.add_material( + MaterialId(0), + LinearElastic::new(E_MOD, NU).with_density(RHO), + None, + ); + db +} + +fn clamp_left(mesh: &Mesh, x_left: f64) -> BoundaryConditionSet { + let clamped: Vec = mesh + .nodes + .iter() + .filter(|(_, node)| (node.position().x - x_left).abs() < 1e-12) + .map(|(&id, _)| id) + .collect(); + let mut set = BoundaryConditionSet::new(); + for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] { + set.add_condition(BoundaryCondition::Dirichlet(DirichletBC { + nodes: clamped.clone(), + components: vec![component], + condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))), + time_range: None, + ramping_factor: 1.0, + gradual_enforcement: false, + })); + } + set +} + +fn point_a(mesh: &Mesh, x: f64, y: f64) -> NodeId { + mesh.nodes + .iter() + .find(|(_, node)| { + (node.position().x - x).abs() < 1e-12 && (node.position().y - y).abs() < 1e-12 + }) + .map(|(&id, _)| id) + .expect("tracking point must be a mesh node") +} + +#[test] +fn fictitious_mass_cancels_with_its_load_and_is_live_without_it() { + let mesh = quad8_rect_mesh(0.25, 0.6, 0.19, 0.21, 10, 2); + let dt = 0.005; + let a_node = point_a(&mesh, 0.6, 0.2); + let tip_force = Vector3::new(0.0, -0.5, 0.0); + let analysis = NonlinearDynamicAnalysis::new( + mesh.clone(), + materials(), + clamp_left(&mesh, 0.25), + dt, + 1, + AnalysisConfig::default(), + ) + .with_total_lagrangian(); + + // A: the plain step from a moving state (three plain steps in, so the + // acceleration is not the rest one). + let mut plain = analysis.stepper().unwrap(); + plain.set_nodal_forces(&[(a_node, tip_force)]); + let mut state = plain.rest_state().unwrap(); + for _ in 0..3 { + state = plain.step(&state).unwrap().0; + } + let (step_a, _) = plain.step(&state).unwrap(); + let a_dofs = plain.node_dofs(a_node); + + // The wetted nodes of this flag: everything not clamped. M_f on all of + // them, the size of the flag's own mass per node (a hard case for the + // cancellation: the added term is O(1) of the inertia). + let wetted: Vec = mesh + .nodes + .iter() + .filter(|(_, n)| (n.position().x - 0.25).abs() > 1e-12) + .map(|(&id, _)| id) + .collect(); + let m_f = 0.5; // kg per node (the flag's 0.02 × 0.35 × 1000 = 7 kg over ~60 nodes) + let entries: Vec<(NodeId, f64)> = wetted.iter().map(|&n| (n, m_f)).collect(); + + // B: the added mass with the compensating load M_f ü_A (ü_A = A's own + // acceleration at each wetted node): the same step to solver tolerance. + let mut fict = analysis.stepper().unwrap(); + fict.set_added_lumped_mass(&entries); + let mut loads = vec![(a_node, tip_force)]; + for &n in &wetted { + let d = fict.node_dofs(n); + loads.push(( + n, + Vector3::new( + m_f * step_a.acceleration[d[0]], + m_f * step_a.acceleration[d[1]], + 0.0, + ), + )); + } + fict.set_nodal_forces(&loads); + let (step_b, _) = fict.step(&state).unwrap(); + let uy_a = step_a.displacement[a_dofs[1]]; + let uy_b = step_b.displacement[a_dofs[1]]; + let du = + (step_a.displacement.clone() - &step_b.displacement).norm() / step_a.displacement.norm(); + println!( + " cancellation: tip uy A {uy_a:.6e} vs B {uy_b:.6e}; relative displacement difference {du:.2e}" + ); + assert!( + du < 1e-6, + "the compensated added mass changed the step: {du:.2e}" + ); + + // C: the added mass WITHOUT the compensating load: a heavier body, a + // visibly different step (the hook is live in the effective matrix). + let mut heavy = analysis.stepper().unwrap(); + heavy.set_added_lumped_mass(&entries); + heavy.set_nodal_forces(&[(a_node, tip_force)]); + let (step_c, _) = heavy.step(&state).unwrap(); + let dc = + (step_a.displacement.clone() - &step_c.displacement).norm() / step_a.displacement.norm(); + println!(" uncompensated: relative displacement difference {dc:.2e}"); + assert!( + dc > 1e-3, + "the added mass alone did not change the step: {dc:.2e}" + ); + + // D: zero added mass is the plain stepper bit for bit. + let mut zero = analysis.stepper().unwrap(); + zero.set_added_lumped_mass(&wetted.iter().map(|&n| (n, 0.0)).collect::>()); + zero.set_nodal_forces(&[(a_node, tip_force)]); + let (step_d, _) = zero.step(&state).unwrap(); + assert_eq!( + step_a.displacement, step_d.displacement, + "zero added mass is not bit-identical" + ); +} diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs index 86d89c0..494d522 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs @@ -100,7 +100,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov // Every setting the acceptance rule reads, printed once: P5-3 lost a // day to a floor of 2e-4 against the overnight marches' 1e-6. println!( - " coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m", + " coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m, fict mass α {}", cfg.coupler, cfg.reuse, cfg.initial_relaxation, @@ -118,6 +118,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov super::overset::patch_stretch(), super::overset::fillet(), super::overset::tip_corner(), + std::env::var("RTX_FSI2O_FICT_MASS").unwrap_or_else(|_| "0".into()), ); // Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the @@ -236,6 +237,57 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov } v }; + let extract_accel = |state: &DynamicState| -> Vec { + let mut a = vec![0.0; 2 * wetted_dofs.len()]; + for (k, dofs) in wetted_dofs.iter().enumerate() { + a[2 * k] = state.acceleration[dofs[0]]; + a[2 * k + 1] = state.acceleration[dofs[1]]; + } + a + }; + // P6-b (`docs/overset_metal_campaign.md` §5.17): the fictitious added + // mass — `RTX_FSI2O_FICT_MASS=α` puts α × ρ_f π (c/2)² (the flag's heave + // added mass per unit depth, c = 0.35) as a lumped mass spread over the + // wetted nodes, and every structure solve carries the compensating load + // M_f ü_k of the previous subiterate, so the fixed point is unchanged + // and the loop contracts at any mass ratio. α = 0 is the plain loop. + let fict_alpha: f64 = std::env::var("RTX_FSI2O_FICT_MASS") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(0.0); + let fict_per_node = fict_alpha * 1000.0 * std::f64::consts::PI * (0.35_f64 / 2.0).powi(2) + / wetted_dofs.len() as f64; + let wetted_nodes: Vec = fluid.interface.wetted.clone(); + if fict_alpha != 0.0 { + let entries: Vec<(NodeId, f64)> = + wetted_nodes.iter().map(|&n| (n, fict_per_node)).collect(); + flag.borrow_mut().set_added_lumped_mass(&entries); + println!( + " fictitious added mass: α = {fict_alpha}, {:.3} kg per wetted node ({} nodes, {:.1} kg total)", + fict_per_node, + wetted_dofs.len(), + fict_per_node * wetted_dofs.len() as f64 + ); + } + // The load with the compensating term for a given previous-subiterate acceleration. + let with_fict = + |nodal: &[(NodeId, Vector3)], accel: &[f64]| -> Vec<(NodeId, Vector3)> { + if fict_alpha == 0.0 { + return nodal.to_vec(); + } + let mut out = nodal.to_vec(); + for (k, &n) in wetted_nodes.iter().enumerate() { + out.push(( + n, + Vector3::new( + fict_per_node * accel[2 * k], + fict_per_node * accel[2 * k + 1], + 0.0, + ), + )); + } + out + }; // Phase 2: release under the current load. let (nodal0, conservation0, faces0) = fluid.sample_load(&zero_d); @@ -293,15 +345,21 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov .unwrap_or(0); let phase_start = std::time::Instant::now(); + // The previous subiterate's interface acceleration (the fictitious mass's + // compensating load): the committed state's at each step's start. + let last_accel: RefCell> = RefCell::new(extract_accel(&flag_state)); for step in 0..coupled_steps { + last_accel.replace(extract_accel(&flag_state)); let d_n = extract(&flag_state); let v_n: Option> = cfg.c1_interface.then(|| extract_velocity(&flag_state)); let d_predicted = if cfg.predictor == "kinematic" { let v = extract_velocity(&flag_state); d_n.iter().zip(&v).map(|(d, v)| d + dt * v).collect() } else { - flag.borrow_mut().set_nodal_forces(&committed_nodal); + flag.borrow_mut() + .set_nodal_forces(&with_fict(&committed_nodal, &extract_accel(&flag_state))); let (predicted, _) = flag.borrow_mut().step(&flag_state).unwrap(); + last_accel.replace(extract_accel(&predicted)); extract(&predicted) }; let saved = fluid.borrow().snapshot(); @@ -317,8 +375,9 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov t_fluid.set(t_fluid.get() + fs.elapsed().as_secs_f64()); let ss = std::time::Instant::now(); let mut flag_ref = flag.borrow_mut(); - flag_ref.set_nodal_forces(&nodal); + flag_ref.set_nodal_forces(&with_fict(&nodal, &last_accel.borrow())); let (candidate_state, _) = flag_ref.step(&flag_state).unwrap(); + last_accel.replace(extract_accel(&candidate_state)); t_structure.set(t_structure.get() + ss.elapsed().as_secs_f64()); let d_new = extract(&candidate_state); if step < cfg.trace_steps {