rtx-fea: reduced Newmark (mor::dynamic) + the phase-4a offline replay — the ≥10x gate is REFUTED by measurement at the validated resolution

The dynamic layer over ReducedNonlinearModel: reduced consistent mass
V'MV (full element sum, never ECSW-sampled — ECSW weights are trained
on internal-force virtual work and would conserve the wrong inertia),
reduced_force_and_jacobian exposed (solve() refactored onto it), and
ReducedNewmark mirroring NonlinearDynamicStepper::newmark_newton in
reduced coordinates (same predictor, residual, tangent shape; no
rescue ladder by design — a reduced Newton death is a finding).

TDD (tests/reduced_newmark.rs): identity-basis march reproduces the
full stepper to 2.4e-14 over 15 steps (both Newton loops tightened to
1e-10 so only solver rounding separates them); rigid-translation
reduced mass = rho*A to 1e-9; a 6-mode POD basis tracks its training
trajectory at 4.2e-4 rms against a 1.0e-4 projection floor.

Phase 4a (fsi3_ecsw_offline.rs, fsi3_reduced_newmark_replay,
env-gated): reduced Newmark replay of the harvested FSI3 trajectory at
record cadence (dt_rec = 5x march dt), driven by the recorded
end-of-step loads. Measured, m=12/20:

- COST (dt-independent, the verdict): 3,068/3,580 us/step at 4.6/5.0
  Newton iters — 2.0-2.3x the banded full-order structural step
  (7,200 us/pass, bandedlu_fsi3_ny62_t85). The >=10x gate needs
  <=720 us/step; one reduced eval alone costs ~640 us because phase 2
  refuted hyperreduction (every eval loops all 70 elements). The gate
  arithmetic is closed: reduced Newton needs >=2 evals, capping the
  ROM at ~5x. THE CAMPAIGN GATE (pinned cycle bands at >=10x
  structural speedup) CANNOT BE MET at the validated resolution.
- TRACKING at record cadence diverges in the release transient (dies
  t=4.35-4.45) — and the RTX_REPLAY_IDENTITY control dies EARLIER
  (t=4.13) in the exact subspace: the death is the 5x-coarse
  integration + aliased loads, NOT the reduction. The record-cadence
  replay cannot judge subspace dynamics; the projection floor
  (1.1e-3 at m=12) remains the honest subspace statement.

Campaign verdict to be recorded in omni-cortex in the pre-registered
words.

Co-Authored-By: Claude Fable 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

crates/specialized/rtx-fea/src/mor/dynamic.rs

@@ -0,0 +1,152 @@
+//! Reduced Newmark: the dynamic layer over [`ReducedNonlinearModel`] —
+//! phase 4 of the ECSW×FSI campaign.
+//!
+//! Mirrors `NonlinearDynamicStepper::newmark_newton` in reduced
+//! coordinates: Newmark's displacement predictor, Newton on the
+//! end-of-step reduced displacement with the tangent
+//! `V'KV + V'MV/(β Δt²)`, velocities and accelerations updated per
+//! Newmark. The reduced mass is projected once at construction (it is
+//! constant); the m×m tangent is factorized dense per iteration — at
+//! m = 12–20 that cost is noise next to the element-loop assembly.
+//!
+//! Deliberately NOT carried over from the full stepper: the rescue
+//! ladder (line search, substepping). The reduced model's first duty is
+//! the offline replay measurement; if reduced plain Newton dies at a
+//! load reversal, that is a phase-4 finding to record, not to paper
+//! over silently.
+
+use super::reduced::ReducedNonlinearModel;
+use crate::analysis::ConvergenceCriteria;
+use crate::error::{AnalysisError, FeaError, FeaResult};
+use nalgebra::{DMatrix, DVector};
+
+/// Reduced-coordinate dynamic state (`q`, `q̇`, `q̈`), each of length
+/// `modes`.
+#[derive(Debug, Clone)]
+pub struct ReducedState {
+    pub q: DVector<f64>,
+    pub q_dot: DVector<f64>,
+    pub q_ddot: DVector<f64>,
+}
+
+/// Newmark time stepping in the reduced coordinates of a
+/// [`ReducedNonlinearModel`].
+pub struct ReducedNewmark<'m, 'a> {
+    model: &'m ReducedNonlinearModel<'a>,
+    /// `V' M V`, full element sum (never ECSW-sampled).
+    mass: DMatrix<f64>,
+    dt: f64,
+    gamma: f64,
+    beta: f64,
+    criteria: ConvergenceCriteria,
+}
+
+impl<'m, 'a> ReducedNewmark<'m, 'a> {
+    /// Average-acceleration Newmark (γ = 1/2, β = 1/4) over `model`,
+    /// with the reduced mass projected here once.
+    pub fn new(model: &'m ReducedNonlinearModel<'a>, dt: f64) -> FeaResult<Self> {
+        Ok(Self {
+            model,
+            mass: model.reduced_mass()?,
+            dt,
+            gamma: 0.5,
+            beta: 0.25,
+            criteria: ConvergenceCriteria::default(),
+        })
+    }
+
+    /// Newmark parameters (default γ = 1/2, β = 1/4).
+    #[must_use]
+    pub fn with_newmark_parameters(mut self, gamma: f64, beta: f64) -> Self {
+        self.gamma = gamma;
+        self.beta = beta;
+        self
+    }
+
+    /// Convergence criteria for the per-step Newton loop (default:
+    /// [`ConvergenceCriteria::default`], as on the full analysis).
+    #[must_use]
+    pub fn with_convergence_criteria(mut self, criteria: ConvergenceCriteria) -> Self {
+        self.criteria = criteria;
+        self
+    }
+
+    /// The state at rest under `external_force_free`: `q = q̇ = 0`, the
+    /// acceleration consistent with `M_r q̈ = V'(F_ext) − f_int_r(0)` —
+    /// the reduced mirror of the full stepper's `rest_state`.
+    pub fn rest_state(&self, external_force_free: &DVector<f64>) -> FeaResult<ReducedState> {
+        let modes = self.model.modes();
+        let q = DVector::zeros(modes);
+        let (internal, _) = self.model.reduced_force_and_jacobian(&q)?;
+        let residual = self.model.reduce_vector(external_force_free) - internal;
+        let q_ddot = self
+            .mass
+            .clone()
+            .lu()
+            .solve(&residual)
+            .ok_or_else(|| FeaError::InvalidInput("singular reduced mass".to_string()))?;
+        Ok(ReducedState {
+            q,
+            q_dot: DVector::zeros(modes),
+            q_ddot,
+        })
+    }
+
+    /// One Newmark step of `dt` from `state` under the end-of-step load
+    /// `external_force_free` (free-DOF space; projected here). Returns
+    /// the end-of-step state and the Newton iteration count. Pure
+    /// function of `(state, load)` — commits nothing.
+    pub fn step(
+        &self,
+        state: &ReducedState,
+        external_force_free: &DVector<f64>,
+    ) -> FeaResult<(ReducedState, usize)> {
+        let dt = self.dt;
+        let (gamma, beta) = (self.gamma, self.beta);
+        let inv_beta_dt2 = 1.0 / (beta * dt * dt);
+
+        let reduced_external = self.model.reduce_vector(external_force_free);
+        let force_scale = reduced_external.norm().max(1.0);
+
+        let q_pred = &state.q + dt * &state.q_dot + dt * dt * (0.5 - beta) * &state.q_ddot;
+
+        let mut q_iter = q_pred.clone();
+        let mut converged = false;
+        let mut iterations = 0usize;
+        for _ in 0..self.criteria.max_iterations {
+            let a_new = inv_beta_dt2 * (&q_iter - &q_pred);
+            let (internal, jacobian) = self.model.reduced_force_and_jacobian(&q_iter)?;
+            let residual = &reduced_external - internal - &self.mass * &a_new;
+            if residual.norm() < self.criteria.force_tolerance * force_scale {
+                converged = true;
+                break;
+            }
+            iterations += 1;
+            let tangent = jacobian + inv_beta_dt2 * &self.mass;
+            let delta = tangent
+                .lu()
+                .solve(&residual)
+                .ok_or_else(|| FeaError::InvalidInput("singular reduced tangent".to_string()))?;
+            let delta_norm = delta.norm();
+            q_iter += delta;
+            if delta_norm < self.criteria.displacement_tolerance * q_iter.norm().max(1.0) {
+                converged = true;
+                break;
+            }
+        }
+        if !converged {
+            return Err(AnalysisError::ConvergenceFailed { iterations }.into());
+        }
+
+        let q_ddot_new = inv_beta_dt2 * (&q_iter - &q_pred);
+        let q_dot_new = &state.q_dot + dt * ((1.0 - gamma) * &state.q_ddot + gamma * &q_ddot_new);
+        Ok((
+            ReducedState {
+                q: q_iter,
+                q_dot: q_dot_new,
+                q_ddot: q_ddot_new,
+            },
+            iterations,
+        ))
+    }
+}

crates/specialized/rtx-fea/src/mor/mod.rs

@@ -21,11 +21,13 @@
 //! `NonlinearDynamicAnalysis::with_total_lagrangian` — required when the
 //! snapshots come from a total-Lagrangian trajectory (the FSI flag).
 
+pub mod dynamic;
 pub mod ecsw;
 pub mod nnls;
 pub mod pod;
 pub mod reduced;
 
+pub use dynamic::{ReducedNewmark, ReducedState};
 pub use ecsw::{EcswModel, ecsw_residual, train_ecsw, train_ecsw_formulated};
 pub use nnls::nnls;
 pub use pod::pod_basis;

crates/specialized/rtx-fea/src/mor/reduced.rs

@@ -160,6 +160,37 @@ impl ElementOperator {
         }
     }
 
+    /// `V_e' M_e V_e` — this element's contribution to the reduced
+    /// consistent mass. The scalar consistent mass is expanded to vector
+    /// DOFs exactly as `NonlinearDynamicStepper` builds `mass_free`.
+    pub(crate) fn reduced_mass(
+        &self,
+        materials: &MaterialDatabase,
+        spatial_dim: usize,
+    ) -> FeaResult<DMatrix<f64>> {
+        let material = materials
+            .get_material(self.material_index)
+            .expect("checked at build time");
+        let density = material.properties().density;
+        let scalar = ElementMatrixComputer::compute_consistent_mass_matrix(
+            &self.finite_element,
+            &self.node_coords,
+            density,
+            None,
+        )?;
+        let nodes = self.node_coords.len();
+        let mut mass = DMatrix::zeros(nodes * spatial_dim, nodes * spatial_dim);
+        for a in 0..nodes {
+            for b in 0..nodes {
+                let m = scalar.matrix[(a, b)];
+                for d in 0..spatial_dim {
+                    mass[(a * spatial_dim + d, b * spatial_dim + d)] = m;
+                }
+            }
+        }
+        Ok(self.local_basis.transpose() * mass * &self.local_basis)
+    }
+
     /// `V_e' f_e(u)` for a full free-DOF displacement — used by training.
     pub(crate) fn reduced_internal_force(
         &self,
@@ -245,6 +276,58 @@ impl<'a> ReducedNonlinearModel<'a> {
         self.active.len()
     }
 
+    /// Number of modes (reduced coordinates).
+    pub fn modes(&self) -> usize {
+        self.basis.ncols()
+    }
+
+    /// `V' v` — project a free-DOF vector into reduced coordinates.
+    pub fn reduce_vector(&self, free_vector: &DVector<f64>) -> DVector<f64> {
+        self.basis.transpose() * free_vector
+    }
+
+    /// `V q` — expand reduced coordinates to the free-DOF space.
+    pub fn expand(&self, q: &DVector<f64>) -> DVector<f64> {
+        &self.basis * q
+    }
+
+    /// The reduced consistent mass `V' M V`, summed over EVERY element —
+    /// never the ECSW sample: ECSW weights are trained on internal-force
+    /// virtual work only, and reweighting the mass with them would
+    /// conserve the wrong inertia.
+    pub fn reduced_mass(&self) -> FeaResult<DMatrix<f64>> {
+        let modes = self.basis.ncols();
+        let spatial_dim = self.mesh.spatial_dimension;
+        let mut mass = DMatrix::zeros(modes, modes);
+        for operator in &self.operators {
+            mass += operator.reduced_mass(self.materials, spatial_dim)?;
+        }
+        Ok(mass)
+    }
+
+    /// The reduced internal force `Σ w_e V_e' f_e(V q)` and tangent
+    /// `Σ w_e V_e' K_e V_e` over the active element set, at reduced
+    /// coordinates `q` — one Newton iteration's assembly, exposed for
+    /// the dynamic (Newmark) driver.
+    pub fn reduced_force_and_jacobian(
+        &self,
+        q: &DVector<f64>,
+    ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
+        let modes = self.basis.ncols();
+        let spatial_dim = self.mesh.spatial_dimension;
+        let free_displacement = &self.basis * q;
+        let mut force = DVector::zeros(modes);
+        let mut jacobian = DMatrix::zeros(modes, modes);
+        for &(index, weight) in &self.active {
+            let operator = &self.operators[index];
+            let local = operator.gather(&free_displacement);
+            let (f, k) = operator.force_and_tangent(self.materials, &local, spatial_dim)?;
+            force += operator.local_basis.transpose() * f * weight;
+            jacobian += operator.local_basis.transpose() * k * &operator.local_basis * weight;
+        }
+        Ok((force, jacobian))
+    }
+
     /// Assemble the weighted reduced internal force `Σ w_e V_e' f_e(u)`
     /// over the active element set at a full free-DOF displacement — the
     /// per-Newton-iteration cost ECSW reduces, exposed for wall-clock
@@ -277,25 +360,13 @@ impl<'a> ReducedNonlinearModel<'a> {
         max_iterations: usize,
     ) -> FeaResult<DVector<f64>> {
         let modes = self.basis.ncols();
-        let spatial_dim = self.mesh.spatial_dimension;
         let reduced_external = self.basis.transpose() * external_force_free;
         let scale = reduced_external.norm().max(1.0);
 
         let mut q = DVector::zeros(modes);
         for _iteration in 0..max_iterations {
-            let free_displacement = &self.basis * &q;
+            let (internal, jacobian) = self.reduced_force_and_jacobian(&q)?;
-
+            let residual = &reduced_external - internal;
-            let mut residual = reduced_external.clone();
-            let mut jacobian = DMatrix::zeros(modes, modes);
-            for &(index, weight) in &self.active {
-                let operator = &self.operators[index];
-                let local = operator.gather(&free_displacement);
-                let (force, tangent) =
-                    operator.force_and_tangent(self.materials, &local, spatial_dim)?;
-                residual -= operator.local_basis.transpose() * force * weight;
-                jacobian +=
-                    operator.local_basis.transpose() * tangent * &operator.local_basis * weight;
-            }
 
             if residual.norm() < force_tolerance * scale {
                 return Ok(&self.basis * &q);

crates/specialized/rtx-fea/tests/reduced_newmark.rs

@@ -0,0 +1,387 @@
+//! The reduced Newmark driver (`mor::dynamic`) against the full-order
+//! stepper — phase 4's machinery, verified before it touches the flag.
+//!
+//! 1. Identity-basis equivalence: with `V = I` over the free DOFs the
+//!    reduced Newmark IS the full Newmark (same predictor, same
+//!    residual, same tangent, different linear-solver rounding), so a
+//!    march must reproduce the full stepper's trajectory to near
+//!    machine precision when both Newton loops are run tight.
+//! 2. The reduced mass carries the right physics: a rigid-translation
+//!    vector must see exactly the total mass ρ·A.
+//! 3. A truncated POD basis built from full-order snapshots must track
+//!    the full trajectory it was trained on to within a band far above
+//!    its projection error but far below any wrong-dynamics answer.
+
+use nalgebra::{DVector, Vector3};
+use rtx_fea::analysis::{AnalysisConfig, ConvergenceCriteria, NonlinearDynamicAnalysis};
+use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
+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};
+use rtx_fea::mor::{Formulation, ReducedNewmark, ReducedNonlinearModel, pod_basis};
+
+const E_MOD: f64 = 1.4e6;
+const NU: f64 = 0.4;
+const RHO: f64 = 1000.0;
+
+/// `nx` by `ny` Quad8 serendipity mesh of `[x0, x1] x [y0, y1]` (as in
+/// `newton_rescue.rs` and the FSI harness).
+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<NodeId> = 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
+}
+
+/// The numbering the analysis uses internally, rebuilt identically
+/// (Sequential strategy, same clamp criterion) — the
+/// `fsi3_ecsw_offline` pattern, self-checked there against the dump.
+fn clamped_numbering(mesh: &Mesh, x_left: f64) -> AdvancedDofNumbering {
+    let mut numbering =
+        AdvancedDofNumbering::displacement_only(mesh, DofMappingStrategy::Sequential).unwrap();
+    for (&node_id, node) in &mesh.nodes {
+        if (node.position().x - x_left).abs() < 1e-12 {
+            for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
+                let dof = numbering.get_dof(node_id, component).unwrap();
+                numbering.constrain_dof(dof).unwrap();
+            }
+        }
+    }
+    numbering
+}
+
+fn tip_node(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("tip node")
+}
+
+/// A nodal force as a free-DOF vector under `numbering`.
+fn free_force(numbering: &AdvancedDofNumbering, loads: &[(NodeId, Vector3<f64>)]) -> DVector<f64> {
+    let mut free_index = vec![None; numbering.total_dofs];
+    for (i, &dof) in numbering.free_dofs.iter().enumerate() {
+        free_index[dof] = Some(i);
+    }
+    let mut force = DVector::zeros(numbering.free_dofs.len());
+    for (node, f) in loads {
+        for (component, &dof) in numbering.get_node_dofs(*node).iter().enumerate() {
+            if let Some(free) = free_index[dof] {
+                force[free] += f[component];
+            }
+        }
+    }
+    force
+}
+
+/// Tight Newton on both sides so the converged states differ only by
+/// linear-solver rounding, not by the stopping tolerance.
+fn tight_criteria() -> ConvergenceCriteria {
+    ConvergenceCriteria {
+        force_tolerance: 1e-10,
+        displacement_tolerance: 1e-12,
+        energy_tolerance: 1e-14,
+        max_iterations: 50,
+    }
+}
+
+#[test]
+fn identity_basis_matches_full_stepper() {
+    let (x0, x1) = (0.0, 0.35);
+    let mesh = quad8_rect_mesh(x0, x1, 0.0, 0.02, 4, 1);
+    let tip = tip_node(&mesh, x1, 0.01);
+    let dt = 1e-3;
+
+    // Full order.
+    let analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        materials(),
+        clamp_left(&mesh, x0),
+        dt,
+        1,
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian()
+    .with_convergence_criteria(tight_criteria());
+    let mut stepper = analysis.stepper().unwrap();
+    let load = [(tip, Vector3::new(0.0, -40.0, 0.0))];
+    stepper.set_nodal_forces(&load);
+    let mut full_state = stepper.rest_state().unwrap();
+
+    // Reduced with V = I over the free DOFs.
+    let db = materials();
+    let numbering = clamped_numbering(&mesh, x0);
+    let n_free = numbering.free_dofs.len();
+    let basis = nalgebra::DMatrix::identity(n_free, n_free);
+    let model = ReducedNonlinearModel::new_formulated(
+        &mesh,
+        &db,
+        &numbering,
+        basis,
+        Formulation::TotalLagrangian,
+    )
+    .unwrap();
+    let newmark = ReducedNewmark::new(&model, dt)
+        .unwrap()
+        .with_convergence_criteria(tight_criteria());
+    let external = free_force(&numbering, &load);
+    let mut reduced_state = newmark.rest_state(&external).unwrap();
+
+    // The consistent initial accelerations must already agree.
+    let a_full: DVector<f64> = DVector::from_iterator(
+        n_free,
+        numbering
+            .free_dofs
+            .iter()
+            .map(|&d| full_state.acceleration[d]),
+    );
+    let a0_err = (&a_full - &reduced_state.q_ddot).norm() / a_full.norm();
+    assert!(
+        a0_err < 1e-9,
+        "rest-state accelerations differ: {a0_err:.3e}"
+    );
+
+    let mut worst = 0.0f64;
+    for step in 0..15 {
+        let (next_full, _) = stepper.step(&full_state).expect("full step");
+        let (next_reduced, _) = newmark
+            .step(&reduced_state, &external)
+            .expect("reduced step");
+        full_state = next_full;
+        reduced_state = next_reduced;
+
+        let d_full: DVector<f64> = DVector::from_iterator(
+            n_free,
+            numbering
+                .free_dofs
+                .iter()
+                .map(|&d| full_state.displacement[d]),
+        );
+        let d_reduced = model.expand(&reduced_state.q);
+        let rel = (&d_full - &d_reduced).norm() / d_full.norm().max(1e-30);
+        worst = worst.max(rel);
+        assert!(
+            rel < 1e-9,
+            "step {step}: identity-basis reduced Newmark left the full trajectory \
+             (rel {rel:.3e})"
+        );
+    }
+    println!("  identity-basis march: worst per-step rel deviation {worst:.3e} over 15 steps");
+}
+
+#[test]
+fn reduced_mass_is_total_mass_on_rigid_translation() {
+    // Unconstrained 2x1 mesh: every DOF free, identity basis — the
+    // reduced mass IS the assembled mass. A rigid x-translation stores
+    // the total mass ρ·A; any quadrature or expansion slip breaks the
+    // total while keeping the matrix symmetric positive definite.
+    let (w, h) = (2.0, 0.5);
+    let mesh = quad8_rect_mesh(0.0, w, 0.0, h, 2, 1);
+    let db = materials();
+    let numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+    let n_free = numbering.free_dofs.len();
+    let basis = nalgebra::DMatrix::identity(n_free, n_free);
+    let model = ReducedNonlinearModel::new_formulated(
+        &mesh,
+        &db,
+        &numbering,
+        basis,
+        Formulation::TotalLagrangian,
+    )
+    .unwrap();
+    let mass = model.reduced_mass().unwrap();
+
+    // x-translation: 1 on every x DOF (free DOF order follows the
+    // numbering; even/odd split by component index within a node).
+    let mut e_x = DVector::zeros(n_free);
+    let mut free_index = vec![None; numbering.total_dofs];
+    for (i, &dof) in numbering.free_dofs.iter().enumerate() {
+        free_index[dof] = Some(i);
+    }
+    for (&node_id, _) in &mesh.nodes {
+        let dofs = numbering.get_node_dofs(node_id);
+        if let Some(free) = free_index[dofs[0]] {
+            e_x[free] = 1.0;
+        }
+    }
+    let total = (e_x.transpose() * &mass * &e_x)[(0, 0)];
+    let expected = RHO * w * h;
+    let rel = (total - expected).abs() / expected;
+    assert!(
+        rel < 1e-9,
+        "rigid x-translation sees {total:.6} against total mass {expected:.6} (rel {rel:.3e})"
+    );
+}
+
+#[test]
+fn truncated_basis_tracks_its_training_trajectory() {
+    let (x0, x1) = (0.0, 0.35);
+    let mesh = quad8_rect_mesh(x0, x1, 0.0, 0.02, 4, 1);
+    let tip = tip_node(&mesh, x1, 0.01);
+    let dt = 1e-3;
+    let steps = 60;
+
+    // Full-order march under a smoothly ramped tip load; snapshot every
+    // step.
+    let analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        materials(),
+        clamp_left(&mesh, x0),
+        dt,
+        1,
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian()
+    .with_convergence_criteria(tight_criteria());
+    let mut stepper = analysis.stepper().unwrap();
+    let db = materials();
+    let numbering = clamped_numbering(&mesh, x0);
+    let n_free = numbering.free_dofs.len();
+    let load_at = |k: usize| {
+        let t = (k as f64) * dt;
+        [(tip, Vector3::new(0.0, -30.0 * (35.0 * t).sin(), 0.0))]
+    };
+
+    stepper.set_nodal_forces(&load_at(0));
+    let mut full_state = stepper.rest_state().unwrap();
+    let mut snapshots: Vec<DVector<f64>> = Vec::with_capacity(steps);
+    let mut full_trajectory: Vec<DVector<f64>> = Vec::with_capacity(steps);
+    for k in 0..steps {
+        stepper.set_nodal_forces(&load_at(k + 1));
+        let (next, _) = stepper.step(&full_state).expect("full step");
+        full_state = next;
+        let d: DVector<f64> = DVector::from_iterator(
+            n_free,
+            numbering
+                .free_dofs
+                .iter()
+                .map(|&d| full_state.displacement[d]),
+        );
+        snapshots.push(d.clone());
+        full_trajectory.push(d);
+    }
+
+    // POD basis from the trajectory, truncated hard.
+    let basis_full = pod_basis(&snapshots, 1e-14).unwrap();
+    let modes = basis_full.ncols().min(6);
+    let basis = basis_full.columns(0, modes).into_owned();
+    // Same normalization as the tracking metric below (absolute error
+    // rms over the max displacement scale) so the two are comparable.
+    let scale = snapshots
+        .iter()
+        .map(nalgebra::DVector::norm)
+        .fold(0.0f64, f64::max);
+    let projection_rms = {
+        let mut sum = 0.0;
+        for d in &snapshots {
+            let err = (d - &basis * (basis.transpose() * d)).norm();
+            sum += err * err;
+        }
+        (sum / snapshots.len() as f64).sqrt() / scale
+    };
+
+    // Reduced march under the same loads from the same rest state.
+    let model = ReducedNonlinearModel::new_formulated(
+        &mesh,
+        &db,
+        &numbering,
+        basis,
+        Formulation::TotalLagrangian,
+    )
+    .unwrap();
+    let newmark = ReducedNewmark::new(&model, dt)
+        .unwrap()
+        .with_convergence_criteria(tight_criteria());
+    let mut reduced_state = newmark
+        .rest_state(&free_force(&numbering, &load_at(0)))
+        .unwrap();
+    let mut sum_sq = 0.0;
+    let mut scale_sq = 0.0f64;
+    for (k, d_full) in full_trajectory.iter().enumerate() {
+        let external = free_force(&numbering, &load_at(k + 1));
+        let (next, _) = newmark
+            .step(&reduced_state, &external)
+            .expect("reduced step");
+        reduced_state = next;
+        let err = (d_full - model.expand(&reduced_state.q)).norm();
+        sum_sq += err * err;
+        scale_sq = scale_sq.max(d_full.norm_squared());
+    }
+    let tracking_rms = (sum_sq / full_trajectory.len() as f64).sqrt() / scale_sq.sqrt();
+    println!(
+        "  {modes}-mode reduced march: tracking rms {tracking_rms:.3e} \
+         (projection rms {projection_rms:.3e})"
+    );
+    // The reduced march may legitimately exceed pure projection error
+    // (closure: the dynamics leave the subspace and come back), but a
+    // wrong mass, wrong formulation, or wrong Newmark constant lands
+    // orders of magnitude higher.
+    assert!(
+        tracking_rms < 1e-2,
+        "reduced march does not track its own training trajectory: rms {tracking_rms:.3e} \
+         (projection floor {projection_rms:.3e})"
+    );
+}

crates/specialized/rtx-fsi/tests/fsi3_ecsw_offline.rs

@@ -26,15 +26,23 @@ use nalgebra::{DMatrix, DVector};
 use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
 use rtx_fea::materials::{LinearElastic, MaterialDatabase};
 use rtx_fea::mesh::MaterialId;
+use rtx_fea::mesh::NodeId;
 use rtx_fea::mor::{
-    Formulation, ReducedNonlinearModel, ecsw_residual, pod_basis, train_ecsw_formulated,
+    Formulation, ReducedNewmark, ReducedNonlinearModel, ecsw_residual, pod_basis,
+    train_ecsw_formulated,
 };
 
 /// One FSNP record (the phase-1 dump format; see `march::write_snapshot`).
 struct Snapshot {
-    #[allow(dead_code)]
     t: f64,
     displacement: Vec<f64>,
+    /// Full-DOF velocity and acceleration — the phase-4 replay's initial
+    /// state.
+    velocity: Vec<f64>,
+    acceleration: Vec<f64>,
+    /// The committed sparse nodal load `(node id, [fx, fy, fz])` — the
+    /// external force the structural step was fed.
+    loads: Vec<(usize, [f64; 3])>,
 }
 
 fn read_fsnp(path: &str) -> (usize, Vec<Snapshot>) {
@@ -50,11 +58,34 @@ fn read_fsnp(path: &str) -> (usize, Vec<Snapshot>) {
     while off < data.len() {
         let t = f64_at(&data, off);
         off += 8;
-        let displacement: Vec<f64> = (0..n_dofs).map(|i| f64_at(&data, off + 8 * i)).collect();
+        let mut series = |off: &mut usize| -> Vec<f64> {
-        off += 8 * n_dofs * 3; // skip velocity + acceleration for this phase
+            let v: Vec<f64> = (0..n_dofs).map(|i| f64_at(&data, *off + 8 * i)).collect();
+            *off += 8 * n_dofs;
+            v
+        };
+        let displacement = series(&mut off);
+        let velocity = series(&mut off);
+        let acceleration = series(&mut off);
         let n_forces = u64::from_le_bytes(data[off..off + 8].try_into().unwrap()) as usize;
-        off += 8 + n_forces * 32;
+        off += 8;
-        records.push(Snapshot { t, displacement });
+        let mut loads = Vec::with_capacity(n_forces);
+        for _ in 0..n_forces {
+            let node = u64::from_le_bytes(data[off..off + 8].try_into().unwrap()) as usize;
+            let f = [
+                f64_at(&data, off + 8),
+                f64_at(&data, off + 16),
+                f64_at(&data, off + 24),
+            ];
+            off += 32;
+            loads.push((node, f));
+        }
+        records.push(Snapshot {
+            t,
+            displacement,
+            velocity,
+            acceleration,
+            loads,
+        });
     }
     assert_eq!(off, data.len(), "trailing bytes");
     (n_dofs, records)
@@ -274,3 +305,192 @@ fn fsi3_ecsw_offline_study() {
         println!("    assembly speedup {:.1}x", t_full / t_sampled);
     }
 }
+
+/// ECSW×FSI campaign, phase 4a: the OFFLINE reduced Newmark replay.
+///
+/// The reduced dynamic model (POD basis, total-Lagrangian operators,
+/// projected consistent mass, plain reduced Newton — `mor::dynamic`,
+/// identity-basis-verified against the full stepper at 2.4e-14) marches
+/// the harvested trajectory's horizon at the RECORD cadence, driven by
+/// the recorded end-of-step nodal loads, from the recorded initial
+/// state. Measured, pre-registered:
+///
+/// 1. Tracking error vs the recorded full-order displacement, release
+///    window and settled cycle separately, against the projection floor
+///    (the best any model in this subspace can do).
+/// 2. Wall-clock per reduced step (the gate's numerator): the banded
+///    full-order structural step measured 7.2 ms/pass on the FSI3
+///    study march (2026-08-30, `bandedlu_fsi3_ny62_t85`), so ≥10×
+///    demands ≤0.72 ms/step here.
+///
+/// Caveat, pre-named: the replay integrates at the record spacing
+/// (5× the march dt), so tracking error conflates subspace closure
+/// with integrator-dt difference; the cost number does not care, and a
+/// model that tracks at THIS dt would only track better at the march's.
+/// A reduced Newton death (no rescue ladder) is a finding, printed and
+/// not papered over.
+#[test]
+fn fsi3_reduced_newmark_replay() {
+    let Ok(snap_path) = std::env::var("RTX_ECSW_SNAP") else {
+        println!("  RTX_ECSW_SNAP not set — reduced Newmark replay skipped");
+        return;
+    };
+    let (n_dofs, records) = read_fsnp(&snap_path);
+    let mesh = flag_mesh(35, 2);
+    let mut materials = MaterialDatabase::new();
+    materials.add_material(
+        MaterialId(0),
+        LinearElastic::new(FSI3.e_s, FSI3.nu_s).with_density(FSI3.rho_s),
+        None,
+    );
+    let dof_numbering = numbering(&mesh);
+    let n_free = dof_numbering.free_dofs.len();
+    assert_eq!(dof_numbering.total_dofs, n_dofs, "DOF count mismatch");
+    let to_free = |full: &[f64]| -> DVector<f64> {
+        DVector::from_iterator(n_free, dof_numbering.free_dofs.iter().map(|&dof| full[dof]))
+    };
+    let mut free_index = vec![None; dof_numbering.total_dofs];
+    for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+        free_index[dof] = Some(i);
+    }
+    let load_to_free = |loads: &[(usize, [f64; 3])]| -> DVector<f64> {
+        let mut force = DVector::zeros(n_free);
+        for &(node, f) in loads {
+            for (component, &dof) in dof_numbering.get_node_dofs(NodeId(node)).iter().enumerate() {
+                if let Some(free) = free_index[dof] {
+                    force[free] += f[component];
+                }
+            }
+        }
+        force
+    };
+
+    // Record cadence (the replay dt). The march's own dt is 5x finer.
+    let spacings: Vec<f64> = records.windows(2).map(|w| w[1].t - w[0].t).collect();
+    let dt_rec = spacings.iter().sum::<f64>() / spacings.len() as f64;
+    let worst_spacing = spacings
+        .iter()
+        .map(|s| (s - dt_rec).abs())
+        .fold(0.0f64, f64::max);
+    println!(
+        "  {} records, dt_rec {dt_rec:.6e} (worst spacing deviation {worst_spacing:.2e})",
+        records.len()
+    );
+    assert!(
+        worst_spacing < 1e-9,
+        "record spacing is not uniform — the fixed-dt replay is invalid"
+    );
+
+    // The basis: trained exactly as phase 2's measurement 1 (interleaved
+    // half), so the projection numbers line up with the campaign record.
+    let train: Vec<DVector<f64>> = records
+        .iter()
+        .step_by(2)
+        .map(|r| to_free(&r.displacement))
+        .collect();
+    let basis_full = pod_basis(&train, 1e-14).unwrap();
+
+    let scale = records
+        .iter()
+        .map(|r| to_free(&r.displacement).norm())
+        .fold(0.0f64, f64::max);
+
+    // The dt-vs-reduction control: RTX_REPLAY_IDENTITY runs the same
+    // replay with the identity basis — the exact subspace, so any death
+    // there is the record-cadence integration (or the plain Newton
+    // without a rescue ladder), NOT the reduction.
+    let bases: Vec<(String, DMatrix<f64>)> = if std::env::var("RTX_REPLAY_IDENTITY").is_ok() {
+        vec![("identity".to_string(), DMatrix::identity(n_free, n_free))]
+    } else {
+        [12usize, 20]
+            .iter()
+            .map(|&m| (format!("m={m}"), basis_full.columns(0, m).into_owned()))
+            .collect()
+    };
+    for (label, v) in bases {
+        let modes = &label;
+
+        // Projection floor over the whole recorded trajectory, same
+        // normalization as the tracking metric.
+        let mut proj_sum = 0.0;
+        for r in &records {
+            let d = to_free(&r.displacement);
+            let err = (&d - &v * (v.transpose() * &d)).norm();
+            proj_sum += err * err;
+        }
+        let proj_rms = (proj_sum / records.len() as f64).sqrt() / scale;
+
+        let model = ReducedNonlinearModel::new_formulated(
+            &mesh,
+            &materials,
+            &dof_numbering,
+            v.clone(),
+            Formulation::TotalLagrangian,
+        )
+        .unwrap();
+        let newmark = ReducedNewmark::new(&model, dt_rec).unwrap();
+
+        // Initial state: the first record, projected.
+        let mut state = rtx_fea::mor::ReducedState {
+            q: v.transpose() * to_free(&records[0].displacement),
+            q_dot: v.transpose() * to_free(&records[0].velocity),
+            q_ddot: v.transpose() * to_free(&records[0].acceleration),
+        };
+
+        let mut sum_sq_release = 0.0f64;
+        let mut n_release = 0usize;
+        let mut sum_sq_cycle = 0.0f64;
+        let mut max_cycle = 0.0f64;
+        let mut n_cycle = 0usize;
+        let mut total_iterations = 0usize;
+        let mut died_at: Option<f64> = None;
+        let start = std::time::Instant::now();
+        for k in 0..records.len() - 1 {
+            let external = load_to_free(&records[k + 1].loads);
+            match newmark.step(&state, &external) {
+                Ok((next, iterations)) => {
+                    state = next;
+                    total_iterations += iterations;
+                }
+                Err(e) => {
+                    died_at = Some(records[k + 1].t);
+                    println!(
+                        "  {modes}: reduced Newton DIED at t = {:.4} (step {k} of {}): {e:?}",
+                        records[k + 1].t,
+                        records.len() - 1
+                    );
+                    break;
+                }
+            }
+            let err = (to_free(&records[k + 1].displacement) - model.expand(&state.q)).norm();
+            if records[k + 1].t < 6.0 {
+                sum_sq_release += err * err;
+                n_release += 1;
+            } else {
+                sum_sq_cycle += err * err;
+                max_cycle = max_cycle.max(err);
+                n_cycle += 1;
+            }
+        }
+        let steps_done = n_release + n_cycle;
+        let per_step = start.elapsed().as_secs_f64() / steps_done.max(1) as f64;
+        println!(
+            "  {modes}: {steps_done} steps, {:.2} Newton iters/step, {:.0} us/step \
+             (full-order structural: 7200 us/pass -> {:.1}x)",
+            total_iterations as f64 / steps_done.max(1) as f64,
+            per_step * 1e6,
+            7.2e-3 / per_step
+        );
+        println!(
+            "    tracking rms: release {:.3e}, cycle {:.3e} (max {:.3e}); projection floor \
+             {:.3e} (of max |d| {scale:.3e} m)",
+            (sum_sq_release / n_release.max(1) as f64).sqrt() / scale,
+            (sum_sq_cycle / n_cycle.max(1) as f64).sqrt() / scale,
+            max_cycle / scale,
+            proj_rms
+        );
+        if let Some(t) = died_at {
+            println!("    DIED at t = {t:.4} — recorded as a phase-4 finding");
+        }
+    }
+}