From 860f5bb37fc75c628759192e166564c3753b659c Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 25 Sep 2026 22:28:53 -0500 Subject: [PATCH 1/2] R8-g: sparse supernodal LDLt for the Newton tangent (default off) rtx_fea::solvers::SparseLdlt: nested-dissection ordering, etree, fundamental supernodes, multifrontal numeric factorisation (blocked LDLt, matrixmultiply gemm, rayon over subtrees, bit-deterministic at any thread count), symbolic analysis reused while the pattern holds. NonlinearDynamicStepper: TangentSolver::{BandedLu (default, unchanged float for float), SparseLdlt { reuse }} via with_tangent_solver or RTX_FEA_TANGENT=sparse[:K]; fixed-pattern CSR assembled from parallel element evaluations (forces summed in the banded path's order); optional modified Newton (factor reuse). The per-element kernel is factored out of assemble() unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../specialized/rtx-fea/src/analysis/mod.rs | 2 + .../rtx-fea/src/analysis/nonlinear_dynamic.rs | 177 ++++-- .../rtx-fea/src/analysis/sparse_tangent.rs | 303 ++++++++++ crates/specialized/rtx-fea/src/solvers/mod.rs | 2 + .../rtx-fea/src/solvers/sparse_ldlt/mod.rs | 369 ++++++++++++ .../src/solvers/sparse_ldlt/numeric.rs | 261 ++++++++ .../src/solvers/sparse_ldlt/ordering.rs | 336 +++++++++++ .../src/solvers/sparse_ldlt/symbolic.rs | 375 ++++++++++++ .../rtx-fea/tests/sparse_tangent.rs | 559 ++++++++++++++++++ .../rtx-fea/tests/sparse_tangent_bench.rs | 89 +++ 10 files changed, 2435 insertions(+), 38 deletions(-) create mode 100644 crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs create mode 100644 crates/specialized/rtx-fea/src/solvers/sparse_ldlt/mod.rs create mode 100644 crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs create mode 100644 crates/specialized/rtx-fea/src/solvers/sparse_ldlt/ordering.rs create mode 100644 crates/specialized/rtx-fea/src/solvers/sparse_ldlt/symbolic.rs create mode 100644 crates/specialized/rtx-fea/tests/sparse_tangent.rs create mode 100644 crates/specialized/rtx-fea/tests/sparse_tangent_bench.rs diff --git a/crates/specialized/rtx-fea/src/analysis/mod.rs b/crates/specialized/rtx-fea/src/analysis/mod.rs index 1dcc2b7..30a3bb7 100644 --- a/crates/specialized/rtx-fea/src/analysis/mod.rs +++ b/crates/specialized/rtx-fea/src/analysis/mod.rs @@ -11,6 +11,7 @@ pub mod flag3d; pub mod modal_analysis; pub mod nonlinear_analysis; pub mod nonlinear_dynamic; +pub mod sparse_tangent; pub mod static_analysis; use crate::assembly::AssemblyOptions; @@ -26,6 +27,7 @@ pub use dynamic_analysis::*; pub use modal_analysis::*; pub use nonlinear_analysis::*; pub use nonlinear_dynamic::*; +pub use sparse_tangent::{TangentSolver, TangentStats}; pub use static_analysis::*; /// Base trait for all finite element analyses. diff --git a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs index 7073086..bc8735f 100644 --- a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs +++ b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs @@ -45,6 +45,7 @@ //! change between steps (and between subiterations of one step) through //! the stepper. +use super::sparse_tangent::{SparseTangent, TangentSolver, TangentStats}; use super::{AnalysisConfig, ConvergenceCriteria}; use crate::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy}; use crate::assembly::SparseMatrix; @@ -56,6 +57,7 @@ use crate::materials::{reduced_constitutive, MaterialDatabase}; use crate::mesh::{Mesh, NodeId}; use crate::solvers::{BandedLu, LinearSolver, SolverOptions}; use nalgebra::{DMatrix, DVector, Vector3}; +use rayon::prelude::*; /// Time histories and final state of a nonlinear transient run. #[derive(Debug, Clone)] @@ -118,6 +120,9 @@ pub struct NonlinearDynamicAnalysis { body_force: Option) -> Vector3 + Send + Sync>>, nodal_forces: Vec<(NodeId, Vector3)>, tracked_nodes: Vec, + /// Explicit tangent-solver choice; `None` = `RTX_FEA_TANGENT`, else + /// the banded LU. + tangent_solver: Option, } impl NonlinearDynamicAnalysis { @@ -145,9 +150,26 @@ impl NonlinearDynamicAnalysis { body_force: None, nodal_forces: Vec::new(), tracked_nodes: Vec::new(), + tangent_solver: None, } } + /// Choose the linear solver of the Newton tangent (default: the + /// banded LU, or the `RTX_FEA_TANGENT` knob when set). See + /// [`TangentSolver`]. + #[must_use] + pub fn with_tangent_solver(mut self, solver: TangentSolver) -> Self { + self.tangent_solver = Some(solver); + self + } + + /// The tangent solver a new stepper will use. + pub fn tangent_solver(&self) -> TangentSolver { + self.tangent_solver + .or_else(TangentSolver::from_env) + .unwrap_or_default() + } + /// Switch to the total-Lagrangian St. Venant–Kirchhoff formulation /// (plane strain in 2-D), as on the nonlinear static analysis. #[must_use] @@ -284,6 +306,8 @@ pub struct NonlinearDynamicStepper<'a> { /// Steps carried by step subdivision after both the plain loop and /// the line search failed. rescued_subdivision: usize, + /// The sparse tangent path (`None` = the banded LU path). + sparse: Option, } impl<'a> NonlinearDynamicStepper<'a> { @@ -420,6 +444,15 @@ impl<'a> NonlinearDynamicStepper<'a> { } let added_mass = DVector::zeros(total_dofs); + let sparse = match analysis.tangent_solver() { + TangentSolver::BandedLu => None, + TangentSolver::SparseLdlt { reuse } => Some(SparseTangent::new( + caches.iter().map(|c| c.dofs.as_slice()), + &free_index, + num_free, + reuse, + )?), + }; let mut stepper = Self { added_mass, analysis, @@ -435,6 +468,7 @@ impl<'a> NonlinearDynamicStepper<'a> { solver_options: SolverOptions::default(), rescued_line_search: 0, rescued_subdivision: 0, + sparse, }; stepper.set_nodal_forces(&analysis.nodal_forces); Ok(stepper) @@ -475,9 +509,16 @@ impl<'a> NonlinearDynamicStepper<'a> { // inv_beta_dt2 is only read when assembling the tangent. let (f_int0, _) = self.assemble(&u, false, 0.0)?; let residual0 = &self.external - &f_int0; - let (a0_free, _) = self - .solver - .solve(&self.mass_free, &residual0, &self.solver_options)?; + let (a0_free, _) = match &mut self.sparse { + None => self + .solver + .solve(&self.mass_free, &residual0, &self.solver_options)?, + Some(sparse) => { + sparse + .mass_solver + .solve(&self.mass_free, &residual0, &self.solver_options)? + } + }; let mut a = DVector::zeros(self.total_dofs); for (k, &dof) in self.free_dofs.iter().enumerate() { a[dof] = a0_free[k]; @@ -602,12 +643,19 @@ impl<'a> NonlinearDynamicStepper<'a> { let mut u_iter = u_pred.clone(); let mut step_converged = false; let mut iterations = 0usize; + // Sparse path only: the previous iteration's residual norm. + let mut previous_norm = f64::INFINITY; for _ in 0..criteria.max_iterations { let mut a_new = DVector::zeros(self.total_dofs); for &dof in &self.free_dofs { a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]); } - let (f_int, tangent) = self.assemble(&u_iter, true, inv_beta_dt2)?; + let (f_int, tangent) = if self.sparse.is_some() { + (self.sparse_forces(&u_iter)?, None) + } else { + let (f_int, tangent) = self.assemble(&u_iter, true, inv_beta_dt2)?; + (f_int, Some(tangent)) + }; let residual = &self.external - &f_int - self.mass_times(&a_new); let residual_norm = residual.norm(); if line_search && !residual_norm.is_finite() { @@ -618,9 +666,28 @@ impl<'a> NonlinearDynamicStepper<'a> { break; } iterations += 1; - let (delta, _) = self - .solver - .solve(&tangent, &residual, &self.solver_options)?; + let delta = match tangent { + Some(tangent) => { + let (delta, _) = + self.solver + .solve(&tangent, &residual, &self.solver_options)?; + delta + } + None => { + let rate = residual_norm / previous_norm; + previous_norm = residual_norm; + let sparse = self.sparse.as_mut().expect("sparse path"); + sparse.solve( + self.caches.iter().map(|c| &c.mass), + &self.added_mass, + &self.free_index, + &residual, + inv_beta_dt2, + line_search, + rate, + )? + } + }; let alpha = if line_search { match self.backtrack(&u_iter, &delta, &u_pred, inv_beta_dt2, residual_norm) { Some(alpha) => alpha, @@ -732,6 +799,33 @@ impl<'a> NonlinearDynamicStepper<'a> { (self.rescued_line_search, self.rescued_subdivision) } + /// Counters and timings of the sparse tangent path (`None` on the + /// banded path). + pub fn tangent_stats(&self) -> Option { + self.sparse.as_ref().map(SparseTangent::stats) + } + + /// Sparse path: every element's force and tangent in parallel, the + /// forces summed in element order (the banded path's order, so the + /// internal force is the same float for float), the tangents kept for + /// the factorisation. + fn sparse_forces(&mut self, solution: &DVector) -> FeaResult> { + let started = std::time::Instant::now(); + let analysis = self.analysis; + let results: Vec<(DVector, DMatrix)> = self + .caches + .par_iter() + .map(|cache| element_force_and_tangent(analysis, cache, solution)) + .collect::>()?; + let sparse = self.sparse.as_mut().expect("sparse path"); + Ok(sparse.accumulate( + results, + self.caches.iter().map(|c| c.dofs.as_slice()), + &self.free_index, + started, + )) + } + /// Internal force and (optionally) tangent at a full displacement /// vector, reduced to the free DOFs. The tangent includes the Newmark /// mass term `M / (β Δt²)` at the caller's `inv_beta_dt2` (the @@ -742,41 +836,11 @@ impl<'a> NonlinearDynamicStepper<'a> { with_tangent: bool, inv_beta_dt2: f64, ) -> FeaResult<(DVector, SparseMatrix)> { - let dim = self.analysis.mesh.spatial_dimension; let num_free = self.free_dofs.len(); let mut internal = DVector::zeros(num_free); let mut tangent = SparseMatrix::new(num_free, num_free); for cache in &self.caches { - let material = self - .analysis - .materials - .get_material(cache.material_id) - .unwrap(); - let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone()); - let mut element_displacement = DVector::zeros(cache.dofs.len()); - for (local, &dof) in cache.dofs.iter().enumerate() { - element_displacement[local] = solution[dof]; - } - let (f_int, k_t) = if self.analysis.total_lagrangian { - let (lambda, mu) = material.properties().lame_parameters(); - let constitutive = saint_venant_kirchhoff(lambda, mu, dim); - total_lagrangian::internal_force_and_tangent( - &fe, - &cache.coords, - &element_displacement, - constitutive.as_ref(), - None, - )? - } else { - let constitutive = reduced_constitutive(material, dim)?; - ElementMatrixComputer::compute_internal_force_and_tangent( - &fe, - &cache.coords, - &element_displacement, - constitutive.as_ref(), - None, - )? - }; + let (f_int, k_t) = element_force_and_tangent(self.analysis, cache, solution)?; for (local_row, &dof_row) in cache.dofs.iter().enumerate() { let Some(free_row) = self.free_index[dof_row] else { continue; @@ -834,3 +898,40 @@ impl<'a> NonlinearDynamicStepper<'a> { out } } + +/// One element's internal force and tangent at the full displacement +/// vector `solution` (total-Lagrangian SVK or small strain, per the +/// analysis) — shared by the banded and the sparse paths. +fn element_force_and_tangent( + analysis: &NonlinearDynamicAnalysis, + cache: &ElementCache, + solution: &DVector, +) -> FeaResult<(DVector, DMatrix)> { + let dim = analysis.mesh.spatial_dimension; + let material = analysis.materials.get_material(cache.material_id).unwrap(); + let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone()); + let mut element_displacement = DVector::zeros(cache.dofs.len()); + for (local, &dof) in cache.dofs.iter().enumerate() { + element_displacement[local] = solution[dof]; + } + if analysis.total_lagrangian { + let (lambda, mu) = material.properties().lame_parameters(); + let constitutive = saint_venant_kirchhoff(lambda, mu, dim); + total_lagrangian::internal_force_and_tangent( + &fe, + &cache.coords, + &element_displacement, + constitutive.as_ref(), + None, + ) + } else { + let constitutive = reduced_constitutive(material, dim)?; + ElementMatrixComputer::compute_internal_force_and_tangent( + &fe, + &cache.coords, + &element_displacement, + constitutive.as_ref(), + None, + ) + } +} diff --git a/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs b/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs new file mode 100644 index 0000000..78f462c --- /dev/null +++ b/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs @@ -0,0 +1,303 @@ +//! The sparse Newton-tangent path of [`super::NonlinearDynamicStepper`] +//! (R8-g): a fixed-pattern CSR assembled from per-element blocks through +//! precomputed positions, factorised by [`SparseLdlt`] with the symbolic +//! analysis done once (the pattern is fixed by the mesh), optionally +//! reusing the numeric factor for several Newton iterations (modified +//! Newton). +//! +//! Default off: [`TangentSolver::BandedLu`] is the stepper's original +//! path, float for float. Select with +//! [`super::NonlinearDynamicAnalysis::with_tangent_solver`] or the +//! environment knob `RTX_FEA_TANGENT=sparse` (`sparse:K` = reuse each +//! factor for up to K Newton solves; `banded` = the default). +//! +//! What stays identical to the banded path: the element kernels, the +//! internal force (summed in the same element/local order), the +//! residual, the convergence tests, the rescue ladder. What differs: the +//! linear solve (LDLᵀ in a nested-dissection order vs banded LU with +//! partial pivoting — rounding-level), and the tangent keeps entries the +//! old `SparseMatrix` assembly drops below 1e-15. + +use crate::assembly::SparseMatrix; +use crate::error::FeaResult; +use crate::solvers::{BandedLu, LinearSolver, SolverOptions, SparseLdlt}; +use nalgebra::{DMatrix, DVector}; +use std::time::Instant; + +/// Which linear solver carries the Newton tangent. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum TangentSolver { + /// The original path: `SparseMatrix` assembly + banded LU every + /// Newton iteration. + #[default] + BandedLu, + /// Fixed-pattern parallel assembly + supernodal sparse LDLᵀ with + /// symbolic reuse. `reuse` = how many Newton solves one numeric + /// factor may serve (1 = full Newton; more = modified Newton, the + /// factor refreshed early when the dt changes, when the residual + /// contracts by less than half per iteration, and always in the + /// line-search rescue). + SparseLdlt { + /// Newton solves per numeric factorisation (≥ 1). + reuse: usize, + }, +} + +impl TangentSolver { + /// Full-Newton sparse path. + pub const SPARSE: Self = Self::SparseLdlt { reuse: 1 }; + + /// Parse the `RTX_FEA_TANGENT` knob: `banded`, `sparse`, `sparse:K`. + /// Unset or unparsable → `None` (unparsable values are logged). + pub fn from_env() -> Option { + let value = std::env::var("RTX_FEA_TANGENT").ok()?; + let value = value.trim(); + let parsed = match value { + "banded" => Some(Self::BandedLu), + "sparse" => Some(Self::SPARSE), + other => other + .strip_prefix("sparse:") + .and_then(|k| k.parse::().ok()) + .filter(|&k| k >= 1) + .map(|reuse| Self::SparseLdlt { reuse }), + }; + if parsed.is_none() { + tracing::warn!("RTX_FEA_TANGENT={value:?} not understood; using the banded LU"); + } + parsed + } +} + +/// Counters and accumulated wall time of the sparse path. +#[derive(Debug, Clone, Copy, Default, PartialEq)] +pub struct TangentStats { + /// Element force/tangent evaluations (parallel) + force summation. + pub assemblies: usize, + /// Numeric factorisations. + pub factorizations: usize, + /// Linear solves (Newton corrections). + pub solves: usize, + /// Symbolic analyses (1 unless the pattern changed). + pub analyses: usize, + /// Solves that fell back to the banded LU (failed pivot). + pub fallbacks: usize, + /// Seconds in element evaluation + force summation. + pub assembly_seconds: f64, + /// Seconds in tangent scatter + numeric factorisation. + pub factor_seconds: f64, + /// Seconds in triangular solves. + pub solve_seconds: f64, + /// Nonzeros of the factor. + pub nnz_l: usize, +} + +/// State of the sparse path, owned by the stepper. +pub(super) struct SparseTangent { + num_free: usize, + row_ptr: Vec, + col_idx: Vec, + values: Vec, + /// Per element, `(local_row · n_e + local_col)` → CSR index, or + /// `u32::MAX` when either DOF is constrained. + elem_pos: Vec>, + /// CSR index of each free DOF's diagonal. + diag_pos: Vec, + /// The element tangents of the latest assembly. + element_tangents: Vec>, + solver: SparseLdlt, + /// Separate instance for the rest state's mass solve. + pub(super) mass_solver: SparseLdlt, + reuse: usize, + uses_left: usize, + factor_coef: f64, + factor_valid: bool, + fallback: BandedLu, + stats: TangentStats, +} + +impl SparseTangent { + /// Build the fixed pattern (every free–free pair sharing an element) + /// and the element position maps; analyse it once. + pub(super) fn new<'e>( + elements: impl Iterator + Clone, + free_index: &[Option], + num_free: usize, + reuse: usize, + ) -> FeaResult { + let mut rows: Vec> = vec![Vec::new(); num_free]; + for dofs in elements.clone() { + for &r in dofs { + let Some(fr) = free_index[r] else { continue }; + for &c in dofs { + if let Some(fc) = free_index[c] { + rows[fr].push(fc); + } + } + } + } + for (fr, row) in rows.iter_mut().enumerate() { + row.push(fr); + row.sort_unstable(); + row.dedup(); + } + let mut row_ptr = Vec::with_capacity(num_free + 1); + row_ptr.push(0); + let mut col_idx = Vec::new(); + for row in &rows { + col_idx.extend_from_slice(row); + row_ptr.push(col_idx.len()); + } + drop(rows); + let find = |r: usize, c: usize| -> usize { + let span = &col_idx[row_ptr[r]..row_ptr[r + 1]]; + row_ptr[r] + + span + .binary_search(&c) + .expect("pattern covers every element pair") + }; + let mut elem_pos = Vec::new(); + for dofs in elements { + let n_e = dofs.len(); + let mut pos = vec![u32::MAX; n_e * n_e]; + for (a, &r) in dofs.iter().enumerate() { + let Some(fr) = free_index[r] else { continue }; + for (b, &c) in dofs.iter().enumerate() { + if let Some(fc) = free_index[c] { + pos[a * n_e + b] = find(fr, fc) as u32; + } + } + } + elem_pos.push(pos); + } + let diag_pos = (0..num_free).map(|d| find(d, d)).collect(); + let mut solver = SparseLdlt::new(); + solver.analyze(num_free, &row_ptr, &col_idx)?; + let values = vec![0.0; col_idx.len()]; + let stats = TangentStats { + analyses: 1, + nnz_l: solver.stats().map_or(0, |s| s.nnz_l), + ..TangentStats::default() + }; + Ok(Self { + num_free, + row_ptr, + col_idx, + values, + elem_pos, + diag_pos, + element_tangents: Vec::new(), + solver, + mass_solver: SparseLdlt::new(), + reuse: reuse.max(1), + uses_left: 0, + factor_coef: f64::NAN, + factor_valid: false, + fallback: BandedLu::new(), + stats, + }) + } + + pub(super) fn stats(&self) -> TangentStats { + self.stats + } + + /// Sum the element internal forces (element order, local order — the + /// banded path's order) into the free-DOF vector and keep the element + /// tangents for a later factorisation. + pub(super) fn accumulate( + &mut self, + element_results: Vec<(DVector, DMatrix)>, + element_dofs: impl Iterator>, + free_index: &[Option], + started: Instant, + ) -> DVector { + let mut internal = DVector::zeros(self.num_free); + self.element_tangents.clear(); + for ((f_e, k_e), dofs) in element_results.into_iter().zip(element_dofs) { + for (local, &dof) in dofs.as_ref().iter().enumerate() { + if let Some(free) = free_index[dof] { + internal[free] += f_e[local]; + } + } + self.element_tangents.push(k_e); + } + self.stats.assemblies += 1; + self.stats.assembly_seconds += started.elapsed().as_secs_f64(); + internal + } + + /// The Newton correction `(K_T + coef·(M + M_added)) δ = r`, with the + /// tangent of the latest [`Self::accumulate`]. `force_fresh` demands a + /// new factor; `rate` is `‖r_i‖/‖r_{i−1}‖` within the step (0 at the + /// first iteration). + pub(super) fn solve<'m>( + &mut self, + masses: impl Iterator>, + added_mass: &DVector, + free_index: &[Option], + residual: &DVector, + coef: f64, + force_fresh: bool, + rate: f64, + ) -> FeaResult> { + #[allow(clippy::float_cmp)] + let same_coef = coef == self.factor_coef; + let refresh = + force_fresh || !self.factor_valid || self.uses_left == 0 || !same_coef || rate > 0.5; + if refresh { + let started = Instant::now(); + self.values.fill(0.0); + for ((k_e, mass), pos) in self.element_tangents.iter().zip(masses).zip(&self.elem_pos) { + let n_e = k_e.nrows(); + for a in 0..n_e { + for b in 0..n_e { + let p = pos[a * n_e + b]; + if p != u32::MAX { + self.values[p as usize] += k_e[(a, b)] + coef * mass[(a, b)]; + } + } + } + } + for (dof, &m) in added_mass.iter().enumerate() { + if m != 0.0 { + if let Some(free) = free_index[dof] { + self.values[self.diag_pos[free]] += coef * m; + } + } + } + let outcome = self.solver.factorize_values(&self.values); + self.stats.factor_seconds += started.elapsed().as_secs_f64(); + if let Err(error) = outcome { + // A failed pivot: this solve goes to the pivoting banded + // LU on the same matrix; the next one tries LDLᵀ again. + if self.stats.fallbacks == 0 { + tracing::warn!("sparse tangent: {error}; falling back to the banded LU"); + } + self.stats.fallbacks += 1; + self.factor_valid = false; + let mut triplets = Vec::with_capacity(self.values.len()); + for row in 0..self.num_free { + for idx in self.row_ptr[row]..self.row_ptr[row + 1] { + triplets.push((row, self.col_idx[idx], self.values[idx])); + } + } + let matrix = SparseMatrix::from_triplets(self.num_free, self.num_free, &triplets)?; + let (x, _) = self + .fallback + .solve(&matrix, residual, &SolverOptions::default())?; + self.stats.solves += 1; + return Ok(x); + } + self.stats.factorizations += 1; + self.factor_valid = true; + self.factor_coef = coef; + self.uses_left = self.reuse; + } + let started = Instant::now(); + let x = self.solver.solve_factored(residual)?; + self.stats.solve_seconds += started.elapsed().as_secs_f64(); + self.stats.solves += 1; + self.uses_left -= 1; + Ok(x) + } +} diff --git a/crates/specialized/rtx-fea/src/solvers/mod.rs b/crates/specialized/rtx-fea/src/solvers/mod.rs index c196cf0..d0f8064 100644 --- a/crates/specialized/rtx-fea/src/solvers/mod.rs +++ b/crates/specialized/rtx-fea/src/solvers/mod.rs @@ -15,6 +15,7 @@ pub mod iterative; #[cfg(all(target_os = "macos", feature = "metal"))] pub mod metal_solvers; pub mod nonlinear; +pub mod sparse_ldlt; pub mod time_integration; #[cfg(disabled)] @@ -34,6 +35,7 @@ pub use iterative::*; #[cfg(all(target_os = "macos", feature = "metal"))] pub use metal_solvers::*; pub use nonlinear::*; +pub use sparse_ldlt::{SparseLdlt, SparseLdltStats}; pub use time_integration::*; /// Solver configuration and options. diff --git a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/mod.rs b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/mod.rs new file mode 100644 index 0000000..85bbbb9 --- /dev/null +++ b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/mod.rs @@ -0,0 +1,369 @@ +// Copyright (c) 2024 RustyTorch++ Team +// Licensed under the Apache License, Version 2.0 + +//! Sparse supernodal LDLᵀ direct solver with symbolic reuse. +//! +//! Built for the Newton tangent of the 3-D flag (R8-g in omni-cortex): +//! `K_T + M/(β Δt²)` is symmetric, its sparsity pattern is fixed by the +//! mesh, and the banded LU it replaces does `O(n·b²)` work on a Hex20 +//! bandwidth of one x-slab (~500 DOFs at 35×2×8). Here: +//! +//! * [`SparseLdlt::analyze`] — nested-dissection ordering, elimination +//! tree, supernodes and every scatter map, once per pattern; +//! * [`SparseLdlt::factorize_csr`] — the numeric multifrontal +//! factorisation on the analysed pattern (re-analysing only when the +//! pattern changes), parallel over independent subtrees and +//! bit-deterministic at any thread count; +//! * [`SparseLdlt::solve_factored`] — the two triangular solves. +//! +//! No pivoting: the matrix must be symmetric and have nonzero pivots in +//! the chosen order (true of an SPD matrix, and of the mildly indefinite +//! Newmark tangents met in practice); a failed pivot is an error the +//! caller can answer with a pivoting solver. Only the lower triangle is +//! read — the upper is assumed equal. + +mod numeric; +mod ordering; +mod symbolic; + +use super::{ConvergenceInfo, LinearSolver, SolverCapabilities, SolverOptions}; +use crate::assembly::SparseMatrix; +use crate::error::{FeaResult, SolverError}; +use nalgebra::DVector; +use std::time::Instant; + +/// Pattern and cost statistics of the current analysis. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SparseLdltStats { + /// Matrix dimension. + pub n: usize, + /// Supernodes in the assembly tree. + pub supernodes: usize, + /// Nonzeros of L including the diagonal. + pub nnz_l: usize, + /// Largest front dimension. + pub max_front: usize, + /// Estimated multiply–adds of one numeric factorisation. + pub work: f64, + /// Symbolic analyses performed so far. + pub analyses: usize, + /// Numeric factorisations performed so far. + pub factorizations: usize, +} + +/// Sparse LDLᵀ with symbolic reuse. See the module docs. +#[derive(Debug, Default)] +pub struct SparseLdlt { + symbolic: Option, + factor: numeric::Factor, + factorized: bool, + analyses: usize, + factorizations: usize, + work: Vec, +} + +impl SparseLdlt { + /// A solver with no analysis yet. + pub fn new() -> Self { + Self::default() + } + + /// Symbolic analysis of an `n × n` CSR pattern (structurally + /// symmetric; sorted column indices per row are not required). + pub fn analyze(&mut self, n: usize, row_ptr: &[usize], col_idx: &[usize]) -> FeaResult<()> { + if row_ptr.len() != n + 1 || row_ptr[n] != col_idx.len() { + return Err(SolverError::FactorizationFailed { + reason: format!( + "sparse LDLt: CSR pattern inconsistent with n = {n} ({} row pointers, {} columns)", + row_ptr.len(), + col_idx.len() + ), + } + .into()); + } + if col_idx.len() > u32::MAX as usize { + return Err(SolverError::FactorizationFailed { + reason: "sparse LDLt: more than 2^32 nonzeros".to_string(), + } + .into()); + } + let sym = symbolic::Symbolic::analyze(n, row_ptr, col_idx); + self.factor.resize(&sym); + self.symbolic = Some(sym); + self.factorized = false; + self.analyses += 1; + Ok(()) + } + + /// Numeric factorisation of the CSR matrix `(row_ptr, col_idx, + /// values)`. Reuses the symbolic analysis when the pattern is the + /// analysed one, re-analyses otherwise. + pub fn factorize_csr( + &mut self, + row_ptr: &[usize], + col_idx: &[usize], + values: &[f64], + ) -> FeaResult<()> { + let n = row_ptr.len().saturating_sub(1); + let reuse = self + .symbolic + .as_ref() + .is_some_and(|s| s.n == n && s.same_pattern(row_ptr, col_idx)); + if !reuse { + self.analyze(n, row_ptr, col_idx)?; + } + self.factorize_values(values) + } + + /// Numeric factorisation of new values on the analysed pattern + /// (the caller guarantees the pattern is unchanged). + pub fn factorize_values(&mut self, values: &[f64]) -> FeaResult<()> { + let sym = self + .symbolic + .as_ref() + .ok_or(SolverError::FactorizationRequired)?; + if values.len() != sym.pattern_col_idx.len() { + return Err(SolverError::FactorizationFailed { + reason: "sparse LDLt: value count differs from the analysed pattern".to_string(), + } + .into()); + } + self.factorized = false; + // Pivot floor relative to the largest diagonal entry. + let mut diag_max = 0.0f64; + for row in 0..sym.n { + for idx in sym.pattern_row_ptr[row]..sym.pattern_row_ptr[row + 1] { + if sym.pattern_col_idx[idx] == row { + diag_max = diag_max.max(values[idx].abs()); + } + } + } + let floor = 1e-14 * diag_max; + numeric::factorize(sym, values, &self.factor, floor).map_err(|fail| { + SolverError::FactorizationFailed { + reason: format!( + "sparse LDLt: pivot {:.3e} at factor column {} (floor {floor:.3e})", + fail.value, fail.column + ), + } + })?; + self.factorized = true; + self.factorizations += 1; + Ok(()) + } + + /// Solve with the current factors. + pub fn solve_factored(&mut self, rhs: &DVector) -> FeaResult> { + let sym = self + .symbolic + .as_ref() + .ok_or(SolverError::FactorizationRequired)?; + if !self.factorized { + return Err(SolverError::FactorizationRequired.into()); + } + if rhs.len() != sym.n { + return Err(SolverError::DimensionMismatch { + matrix_rows: sym.n, + matrix_cols: sym.n, + rhs_rows: rhs.len(), + rhs_cols: 1, + } + .into()); + } + self.work.clear(); + self.work.extend(sym.perm.iter().map(|&old| rhs[old])); + numeric::solve_in_place(sym, &self.factor, &mut self.work); + let mut x = DVector::zeros(sym.n); + for (new, &old) in sym.perm.iter().enumerate() { + x[old] = self.work[new]; + } + Ok(x) + } + + /// Whether a numeric factor is available. + pub fn is_factorized(&self) -> bool { + self.factorized + } + + /// Statistics of the current analysis (`None` before the first). + pub fn stats(&self) -> Option { + self.symbolic.as_ref().map(|s| SparseLdltStats { + n: s.n, + supernodes: s.num_supernodes(), + nnz_l: s.nnz_l, + max_front: s.max_front(), + work: s.work, + analyses: self.analyses, + factorizations: self.factorizations, + }) + } +} + +impl LinearSolver for SparseLdlt { + fn solve( + &mut self, + matrix: &SparseMatrix, + rhs: &DVector, + _options: &SolverOptions, + ) -> FeaResult<(DVector, ConvergenceInfo)> { + let start_time = Instant::now(); + if matrix.nrows() != matrix.ncols() || matrix.nrows() != rhs.len() { + return Err(SolverError::DimensionMismatch { + matrix_rows: matrix.nrows(), + matrix_cols: matrix.ncols(), + rhs_rows: rhs.len(), + rhs_cols: 1, + } + .into()); + } + let (row_ptr, col_idx) = matrix.structure(); + if row_ptr.len() != matrix.nrows() + 1 { + return Err(SolverError::FactorizationFailed { + reason: "sparse LDLt needs a finalized (CSR) matrix".to_string(), + } + .into()); + } + self.factorize_csr(row_ptr, col_idx, matrix.values())?; + let x = self.solve_factored(rhs)?; + let mut info = ConvergenceInfo::new(); + info.set_solve_time(start_time.elapsed()); + info.set_converged(1, 0.0, 0.0); + Ok((x, info)) + } + + fn name(&self) -> &'static str { + "Sparse supernodal LDLt" + } + + fn capabilities(&self) -> SolverCapabilities { + SolverCapabilities { + symmetric: true, + positive_definite: false, + gpu_acceleration: false, + multiple_rhs: true, + iterative_refinement: false, + memory_efficiency: 5, + computational_efficiency: 5, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::solvers::BandedLu; + + /// A 3-D 7-point Laplacian with `dofs` components per node coupled + /// by a small symmetric block, shifted by `shift` (negative shifts + /// make it indefinite). + fn laplacian_3d(nx: usize, ny: usize, nz: usize, dofs: usize, shift: f64) -> SparseMatrix { + let nn = nx * ny * nz; + let n = nn * dofs; + let id = |i: usize, j: usize, k: usize| (k * ny + j) * nx + i; + let mut m = SparseMatrix::new(n, n); + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let p = id(i, j, k); + let mut nbrs = Vec::new(); + if i > 0 { + nbrs.push(id(i - 1, j, k)); + } + if i + 1 < nx { + nbrs.push(id(i + 1, j, k)); + } + if j > 0 { + nbrs.push(id(i, j - 1, k)); + } + if j + 1 < ny { + nbrs.push(id(i, j + 1, k)); + } + if k > 0 { + nbrs.push(id(i, j, k - 1)); + } + if k + 1 < nz { + nbrs.push(id(i, j, k + 1)); + } + for a in 0..dofs { + let row = p * dofs + a; + m.add_entry(row, row, 6.0 + shift + 0.1 * a as f64).unwrap(); + for b in 0..dofs { + if a != b { + m.add_entry(row, p * dofs + b, 0.3).unwrap(); + } + } + for &q in &nbrs { + m.add_entry(row, q * dofs + a, -1.0).unwrap(); + for b in 0..dofs { + if a != b { + m.add_entry(row, q * dofs + b, 0.05).unwrap(); + } + } + } + } + } + } + } + m.finalize().unwrap(); + m + } + + fn check(m: &SparseMatrix) { + let n = m.nrows(); + let x_exact = DVector::from_fn(n, |i, _| ((i as f64) * 0.37).sin() + 0.2); + let b = m.multiply_vector(&x_exact).unwrap(); + let opts = SolverOptions::default(); + let (x, _) = SparseLdlt::new().solve(m, &b, &opts).unwrap(); + let rel = (&x - &x_exact).norm() / x_exact.norm(); + assert!(rel < 1e-11, "manufactured rel err {rel:.3e}"); + let (xb, _) = BandedLu::new().solve(m, &b, &opts).unwrap(); + let rel_b = (&x - &xb).norm() / xb.norm(); + assert!(rel_b < 1e-11, "vs banded LU rel err {rel_b:.3e}"); + } + + #[test] + fn spd_3d_block_laplacian() { + check(&laplacian_3d(12, 3, 5, 3, 0.0)); + } + + #[test] + fn indefinite_shifted_laplacian() { + check(&laplacian_3d(9, 4, 4, 2, -5.3)); + } + + #[test] + fn tiny_and_diagonal() { + check(&laplacian_3d(1, 1, 1, 1, 0.0)); + check(&laplacian_3d(3, 1, 1, 1, 0.0)); + } + + #[test] + fn symbolic_reuse_across_values() { + let m1 = laplacian_3d(10, 3, 4, 3, 0.0); + let m2 = laplacian_3d(10, 3, 4, 3, 1.5); + let mut solver = SparseLdlt::new(); + let opts = SolverOptions::default(); + for m in [&m1, &m2, &m1] { + let x_exact = DVector::from_fn(m.nrows(), |i, _| (i as f64 * 0.11).cos()); + let b = m.multiply_vector(&x_exact).unwrap(); + let (x, _) = solver.solve(m, &b, &opts).unwrap(); + assert!((&x - &x_exact).norm() / x_exact.norm() < 1e-11); + } + let stats = solver.stats().unwrap(); + assert_eq!(stats.analyses, 1); + assert_eq!(stats.factorizations, 3); + } + + #[test] + fn zero_pivot_is_an_error() { + let mut m = SparseMatrix::new(2, 2); + m.add_entry(0, 1, 1.0).unwrap(); + m.add_entry(1, 0, 1.0).unwrap(); + m.add_entry(1, 1, 1.0).unwrap(); + m.finalize().unwrap(); + let b = DVector::from_vec(vec![1.0, 2.0]); + // The (0,0) pivot is structurally absent: zero. + let result = SparseLdlt::new().solve(&m, &b, &SolverOptions::default()); + assert!(result.is_err()); + } +} diff --git a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs new file mode 100644 index 0000000..9cc7a26 --- /dev/null +++ b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs @@ -0,0 +1,261 @@ +//! Numeric phase: multifrontal supernodal LDLᵀ (no pivoting) over the +//! symbolic assembly tree, and the triangular solves. +//! +//! Each supernode's front is a dense column-major `m × m` block (lower +//! triangle used): the original entries are scattered in, the children's +//! update matrices extend-added (in child order), then the `k` pivot +//! columns are factorised by a blocked right-looking LDLᵀ whose trailing +//! updates are `matrixmultiply` gemms (through nalgebra views). Sibling +//! subtrees run on rayon when their estimated work is large; the +//! arithmetic of every supernode is independent of the thread schedule, +//! so results are bit-identical at any thread count. + +use super::symbolic::Symbolic; +use nalgebra::{DMatrixView, DMatrixViewMut}; +use rayon::prelude::*; +use std::sync::Mutex; + +/// Panel width of the blocked partial factorisation. +const PANEL: usize = 32; +/// Column block of the lower-triangular trailing update. +const TRAIL_BLOCK: usize = 64; +/// Subtrees below this many multiply–adds are factorised serially. +const PARALLEL_WORK: f64 = 2.0e6; + +/// A zero, non-finite or relatively tiny pivot. +#[derive(Debug, Clone)] +pub(super) struct PivotFailure { + pub column: usize, + pub value: f64, +} + +/// Per-supernode factor blocks: `m × k` column-major, unit-lower L with +/// D on the diagonal. Allocations persist across factorisations. +#[derive(Debug, Default)] +pub(super) struct Factor { + pub blocks: Vec>>, +} + +impl Factor { + pub fn resize(&mut self, sym: &Symbolic) { + let ns = sym.num_supernodes(); + if self.blocks.len() != ns { + self.blocks = (0..ns).map(|_| Mutex::new(Vec::new())).collect(); + } + } +} + +/// Factorise the matrix whose CSR values are `values` (on the analysed +/// pattern). `pivot_floor` rejects |d| below it. +pub(super) fn factorize( + sym: &Symbolic, + values: &[f64], + factor: &Factor, + pivot_floor: f64, +) -> Result<(), PivotFailure> { + let ctx = Ctx { + sym, + values, + factor, + pivot_floor, + }; + let results: Vec, PivotFailure>> = if sym.roots.len() > 1 { + sym.roots.par_iter().map(|&r| ctx.subtree(r)).collect() + } else { + sym.roots.iter().map(|&r| ctx.subtree(r)).collect() + }; + for r in results { + r?; + } + Ok(()) +} + +struct Ctx<'a> { + sym: &'a Symbolic, + values: &'a [f64], + factor: &'a Factor, + pivot_floor: f64, +} + +impl Ctx<'_> { + /// Factorise the subtree rooted at `s`; returns its update matrix + /// (`(m − k)²`, column-major). + fn subtree(&self, s: usize) -> Result, PivotFailure> { + let sym = self.sym; + let kids = sym.node_children(s); + let updates: Vec, PivotFailure>> = + if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_WORK { + kids.par_iter().map(|&c| self.subtree(c)).collect() + } else { + kids.iter().map(|&c| self.subtree(c)).collect() + }; + + let rows = sym.front_rows(s); + let m = rows.len(); + let f = sym.start[s]; + let k = sym.start[s + 1] - f; + let mut front = vec![0.0f64; m * m]; + for &(idx, off) in &sym.amap[sym.amap_ptr[s]..sym.amap_ptr[s + 1]] { + front[off as usize] += self.values[idx as usize]; + } + for (&c, update) in kids.iter().zip(updates) { + let update = update?; + let rel = &sym.relind[sym.relind_ptr[c]..sym.relind_ptr[c + 1]]; + let u = rel.len(); + for b in 0..u { + let col = rel[b] as usize * m; + let src = &update[b * u..(b + 1) * u]; + for a in b..u { + front[rel[a] as usize + col] += src[a]; + } + } + } + partial_ldlt(&mut front, m, k, self.pivot_floor).map_err(|(j, value)| PivotFailure { + column: f + j, + value, + })?; + { + let mut block = self.factor.blocks[s].lock().unwrap(); + block.clear(); + block.extend_from_slice(&front[..m * k]); + } + let u = m - k; + let mut update = vec![0.0f64; u * u]; + for b in 0..u { + let src = (k + b) * m + k; + update[b * u + b..(b + 1) * u].copy_from_slice(&front[src + b..src + u]); + } + Ok(update) + } +} + +/// Blocked right-looking LDLᵀ of the first `k` columns of the `m × m` +/// column-major front, updating the trailing `(m−k)²` lower block. +/// On exit columns `0..k` hold L (unit, strictly below the diagonal) and +/// D (on the diagonal). +fn partial_ldlt(front: &mut [f64], m: usize, k: usize, floor: f64) -> Result<(), (usize, f64)> { + let mut w = Vec::new(); + let mut p0 = 0; + while p0 < k { + let p1 = (p0 + PANEL).min(k); + // Unblocked factorisation of the panel's columns over rows p0..m. + for j in p0..p1 { + let d = front[j + j * m]; + if !d.is_finite() || d.abs() <= floor { + return Err((j, d)); + } + let inv_d = 1.0 / d; + for c in (j + 1)..p1 { + let lcj_d = front[c + j * m] * inv_d; + if lcj_d != 0.0 { + let (left, right) = front.split_at_mut(c * m); + let src = &left[j * m + c..j * m + m]; + let dst = &mut right[c..m]; + for (x, &y) in dst.iter_mut().zip(src) { + *x -= lcj_d * y; + } + } + } + for x in &mut front[j * m + j + 1..j * m + m] { + *x *= inv_d; + } + } + // Trailing update of columns p1..m (lower): F -= L_P D_P L_Pᵀ. + if p1 < m { + let nb = p1 - p0; + let r = m - p1; + // W = L_P · D_P, rows p1..m, column-major r × nb. + w.clear(); + w.resize(r * nb, 0.0); + for (t, j) in (p0..p1).enumerate() { + let d = front[j + j * m]; + for (x, &l) in w[t * r..(t + 1) * r] + .iter_mut() + .zip(&front[j * m + p1..j * m + m]) + { + *x = l * d; + } + } + let (left, right) = front.split_at_mut(p1 * m); + let mut c0 = 0; + while c0 < r { + let c1 = (c0 + TRAIL_BLOCK).min(r); + let rows = r - c0; + // L_P rows (p1 + c0)..m: a view into `left` (columns p0..p1). + let lp = DMatrixView::from_slice_with_strides_generic( + &left[p0 * m + p1 + c0..], + nalgebra::Dyn(rows), + nalgebra::Dyn(nb), + nalgebra::Dyn(1), + nalgebra::Dyn(m), + ); + // W rows c0..c1 (the columns being updated), transposed. + let wb = DMatrixView::from_slice_with_strides_generic( + &w[c0..], + nalgebra::Dyn(c1 - c0), + nalgebra::Dyn(nb), + nalgebra::Dyn(1), + nalgebra::Dyn(r), + ); + let mut target = DMatrixViewMut::from_slice_with_strides_generic( + &mut right[c0 * m + p1 + c0..], + nalgebra::Dyn(rows), + nalgebra::Dyn(c1 - c0), + nalgebra::Dyn(1), + nalgebra::Dyn(m), + ); + target.gemm(-1.0, &lp, &wb.transpose(), 1.0); + c0 = c1; + } + } + p0 = p1; + } + Ok(()) +} + +/// Solve `L D Lᵀ x = b` in the factor numbering, in place. +pub(super) fn solve_in_place(sym: &Symbolic, factor: &Factor, x: &mut [f64]) { + let ns = sym.num_supernodes(); + let blocks: Vec>> = + factor.blocks.iter().map(|b| b.lock().unwrap()).collect(); + for s in 0..ns { + let rows = sym.front_rows(s); + let m = rows.len(); + let f = sym.start[s]; + let k = sym.start[s + 1] - f; + let block = &blocks[s]; + for j in 0..k { + let xj = x[f + j]; + if xj != 0.0 { + let col = &block[j * m..(j + 1) * m]; + for i in (j + 1)..m { + x[rows[i]] -= col[i] * xj; + } + } + } + } + for s in 0..ns { + let m = sym.rows_ptr[s + 1] - sym.rows_ptr[s]; + let f = sym.start[s]; + let k = sym.start[s + 1] - f; + let block = &blocks[s]; + for j in 0..k { + x[f + j] /= block[j + j * m]; + } + } + for s in (0..ns).rev() { + let rows = sym.front_rows(s); + let m = rows.len(); + let f = sym.start[s]; + let k = sym.start[s + 1] - f; + let block = &blocks[s]; + for j in (0..k).rev() { + let col = &block[j * m..(j + 1) * m]; + let mut acc = x[f + j]; + for i in (j + 1)..m { + acc -= col[i] * x[rows[i]]; + } + x[f + j] = acc; + } + } +} diff --git a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/ordering.rs b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/ordering.rs new file mode 100644 index 0000000..d4f2b35 --- /dev/null +++ b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/ordering.rs @@ -0,0 +1,336 @@ +//! Fill-reducing ordering: nested dissection on the matrix graph. +//! +//! The graph is first compressed to supervariables (DOFs with identical +//! closed neighbourhoods — the 2 or 3 displacement components of one FE +//! node), then dissected recursively: a pseudo-peripheral vertex's BFS +//! level structure supplies the candidate separators (one level each, +//! trimmed to the vertices that actually touch the far side), the +//! lightest one that keeps both halves at least 30 % of the weight is +//! taken, and the ordering is `A, B, S`. On a long, thin FE body (the +//! flag) the level sets are cross-sections, so the top separators are the +//! minimal ones; small subgraphs are left in natural order. +//! +//! Deterministic: no hashing order, no randomness, no threads. + +use std::collections::HashMap; + +/// Subgraphs at or below this weight (DOFs) are not dissected further. +const LEAF_WEIGHT: usize = 96; + +/// Compressed adjacency: vertex `v`'s neighbours are +/// `adj[ptr[v]..ptr[v + 1]]` (no self loops), its weight the number of +/// DOFs it stands for. +struct Graph { + ptr: Vec, + adj: Vec, + weight: Vec, +} + +impl Graph { + fn neighbours(&self, v: u32) -> &[u32] { + &self.adj[self.ptr[v as usize]..self.ptr[v as usize + 1]] + } +} + +/// A nested-dissection permutation of the `n × n` symmetric pattern +/// (`perm[new] = old`). Only the pattern's structure is read; a +/// nonsymmetric pattern is symmetrised. +pub(super) fn nested_dissection(n: usize, row_ptr: &[usize], col_idx: &[usize]) -> Vec { + if n == 0 { + return Vec::new(); + } + // Symmetrised adjacency lists over DOFs, self loops excluded. + let mut lists: Vec> = vec![Vec::new(); n]; + for row in 0..n { + for &col in &col_idx[row_ptr[row]..row_ptr[row + 1]] { + if col != row { + lists[row].push(col as u32); + lists[col].push(row as u32); + } + } + } + for list in &mut lists { + list.sort_unstable(); + list.dedup(); + } + + // Supervariables: identical closed neighbourhoods. + let mut group_of = vec![0u32; n]; + let mut members: Vec> = Vec::new(); + { + let mut index: HashMap, u32> = HashMap::new(); + for v in 0..n { + let mut closed = lists[v].clone(); + let at = closed.partition_point(|&x| x < v as u32); + closed.insert(at, v as u32); + let next = members.len() as u32; + let g = *index.entry(closed).or_insert(next); + if g == next { + members.push(Vec::new()); + } + members[g as usize].push(v as u32); + group_of[v] = g; + } + } + let ng = members.len(); + let mut ptr = Vec::with_capacity(ng + 1); + let mut adj = Vec::new(); + ptr.push(0); + for group in &members { + let rep = group[0] as usize; + let mut nbrs: Vec = lists[rep] + .iter() + .map(|&u| group_of[u as usize]) + .filter(|&g| g != group_of[rep]) + .collect(); + nbrs.sort_unstable(); + nbrs.dedup(); + adj.extend_from_slice(&nbrs); + ptr.push(adj.len()); + } + drop(lists); + let graph = Graph { + ptr, + adj, + weight: members.iter().map(Vec::len).collect(), + }; + + let mut state = Dissector { + graph: &graph, + stamp: vec![0u32; ng], + level: vec![u32::MAX; ng], + next_stamp: 0, + out: Vec::with_capacity(ng), + }; + let all: Vec = (0..ng as u32).collect(); + state.dissect(all); + + let mut perm = Vec::with_capacity(n); + for &g in &state.out { + for &v in &members[g as usize] { + perm.push(v as usize); + } + } + debug_assert_eq!(perm.len(), n); + perm +} + +struct Dissector<'g> { + graph: &'g Graph, + /// `stamp[v] == id` marks membership of the subgraph being processed. + stamp: Vec, + level: Vec, + next_stamp: u32, + out: Vec, +} + +impl Dissector<'_> { + fn mark(&mut self, verts: &[u32]) -> u32 { + self.next_stamp += 1; + let id = self.next_stamp; + for &v in verts { + self.stamp[v as usize] = id; + } + id + } + + fn weight(&self, verts: &[u32]) -> usize { + verts.iter().map(|&v| self.graph.weight[v as usize]).sum() + } + + /// BFS inside the stamped subgraph from `start`; fills `level` for + /// the reached vertices and returns them grouped by level. + fn bfs(&mut self, id: u32, start: u32) -> Vec> { + let mut levels: Vec> = vec![vec![start]]; + self.level[start as usize] = 0; + let mut visited = vec![start]; + loop { + let depth = levels.len() as u32; + let mut next = Vec::new(); + for &v in levels.last().unwrap() { + for &u in self.graph.neighbours(v) { + if self.stamp[u as usize] == id && self.level[u as usize] == u32::MAX { + self.level[u as usize] = depth; + next.push(u); + visited.push(u); + } + } + } + if next.is_empty() { + break; + } + next.sort_unstable(); + levels.push(next); + } + for &v in &visited { + self.level[v as usize] = u32::MAX; + } + levels + } + + fn dissect(&mut self, verts: Vec) { + if self.weight(&verts) <= LEAF_WEIGHT || verts.len() < 4 { + self.out.extend_from_slice(&verts); + return; + } + let id = self.mark(&verts); + // Connected components first. + let first = self.bfs(id, verts[0]); + let reached: usize = first.iter().map(Vec::len).sum(); + if reached < verts.len() { + let mut component: Vec = first.into_iter().flatten().collect(); + component.sort_unstable(); + let comp_id = self.mark(&component); + let rest: Vec = verts + .iter() + .copied() + .filter(|&v| self.stamp[v as usize] != comp_id) + .collect(); + self.dissect(component); + self.dissect(rest); + return; + } + + // Pseudo-peripheral start: min-degree vertex of the last level, + // repeated while the eccentricity grows. + let degree = |g: &Graph, v: u32| g.ptr[v as usize + 1] - g.ptr[v as usize]; + let mut start = *verts + .iter() + .min_by_key(|&&v| (degree(self.graph, v), v)) + .unwrap(); + let mut levels = self.bfs(id, start); + for _ in 0..8 { + let far = *levels + .last() + .unwrap() + .iter() + .min_by_key(|&&v| (degree(self.graph, v), v)) + .unwrap(); + let trial = self.bfs(id, far); + if trial.len() > levels.len() { + start = far; + levels = trial; + } else { + break; + } + } + let _ = start; + if levels.len() < 3 { + self.out.extend_from_slice(&verts); + return; + } + + for (depth, level) in levels.iter().enumerate() { + for &v in level { + self.level[v as usize] = depth as u32; + } + } + let total = self.weight(&verts); + let level_weight: Vec = levels.iter().map(|l| self.weight(l)).collect(); + // Candidate separators around the weighted median level. + let mut best: Option<(usize, usize)> = None; // (separator weight, level) + let mut before = 0usize; + for l in 0..levels.len() { + let after = total - before - level_weight[l]; + if l >= 1 && l + 1 < levels.len() { + let sep: usize = levels[l] + .iter() + .filter(|&&v| { + self.graph.neighbours(v).iter().any(|&u| { + self.stamp[u as usize] == id && self.level[u as usize] == l as u32 + 1 + }) + }) + .map(|&v| self.graph.weight[v as usize]) + .sum(); + let side_a = before + level_weight[l] - sep; + let balanced = 10 * side_a.min(after) >= 3 * total; + if balanced && best.is_none_or(|(w, _)| sep < w) { + best = Some((sep, l)); + } + } + before += level_weight[l]; + } + let cut = match best { + Some((_, l)) => l, + None => { + // No balanced level: take the weighted median. + let mut acc = 0usize; + let mut pick = levels.len() / 2; + for (l, w) in level_weight.iter().enumerate() { + acc += w; + if 2 * acc >= total { + pick = l; + break; + } + } + pick.clamp(1, levels.len() - 2) + } + }; + let mut part_a = Vec::new(); + let mut part_b = Vec::new(); + let mut separator = Vec::new(); + for &v in &verts { + let l = self.level[v as usize] as usize; + if l < cut { + part_a.push(v); + } else if l > cut { + part_b.push(v); + } else if self + .graph + .neighbours(v) + .iter() + .any(|&u| self.stamp[u as usize] == id && self.level[u as usize] == cut as u32 + 1) + { + separator.push(v); + } else { + part_a.push(v); + } + } + for &v in &verts { + self.level[v as usize] = u32::MAX; + } + self.dissect(part_a); + self.dissect(part_b); + self.out.extend_from_slice(&separator); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn grid_pattern(nx: usize, ny: usize) -> (usize, Vec, Vec) { + let n = nx * ny; + let mut row_ptr = vec![0]; + let mut col_idx = Vec::new(); + for j in 0..ny { + for i in 0..nx { + let mut cols = Vec::new(); + for (di, dj) in [(-1i64, 0i64), (1, 0), (0, -1), (0, 1), (0, 0)] { + let (a, b) = (i as i64 + di, j as i64 + dj); + if a >= 0 && b >= 0 && (a as usize) < nx && (b as usize) < ny { + cols.push(b as usize * nx + a as usize); + } + } + cols.sort_unstable(); + col_idx.extend(cols); + row_ptr.push(col_idx.len()); + } + } + (n, row_ptr, col_idx) + } + + #[test] + fn nested_dissection_is_a_permutation() { + let (n, rp, ci) = grid_pattern(40, 7); + let perm = nested_dissection(n, &rp, &ci); + let mut seen = vec![false; n]; + for &p in &perm { + assert!(!seen[p]); + seen[p] = true; + } + assert_eq!(perm.len(), n); + // Deterministic. + assert_eq!(perm, nested_dissection(n, &rp, &ci)); + } +} diff --git a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/symbolic.rs b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/symbolic.rs new file mode 100644 index 0000000..59a0e31 --- /dev/null +++ b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/symbolic.rs @@ -0,0 +1,375 @@ +//! Symbolic analysis: ordering, elimination tree, postorder, column +//! structures, fundamental supernodes, the supernodal assembly tree, and +//! the maps the numeric phase scatters through. Computed once per +//! sparsity pattern and reused for every numeric factorisation. + +use super::ordering::nested_dissection; + +/// Everything the numeric factorisation needs that depends only on the +/// pattern. Column indices below are in the *factor* numbering +/// (`perm[new] = old`). +#[derive(Debug)] +pub(super) struct Symbolic { + pub n: usize, + pub perm: Vec, + /// Supernode `s` owns columns `start[s]..start[s + 1]`. + pub start: Vec, + /// Front row indices of supernode `s` (its own columns first, then the + /// sorted off-diagonal structure): `rows[rows_ptr[s]..rows_ptr[s+1]]`. + pub rows_ptr: Vec, + pub rows: Vec, + /// Children in the assembly tree, increasing (postorder) order. + pub children_ptr: Vec, + pub children: Vec, + /// Roots of the assembly forest. + pub roots: Vec, + /// For a non-root supernode: position of each of its update rows + /// (`rows[k..m]`) in its parent's front — `relind[relind_ptr[s]..]`. + pub relind_ptr: Vec, + pub relind: Vec, + /// Per supernode, the lower-triangle entries of A it assembles: + /// `(CSR index, column-major offset in the front)`. + pub amap_ptr: Vec, + pub amap: Vec<(u32, u32)>, + /// The pattern this analysis was built for (for reuse checks). + pub pattern_row_ptr: Vec, + pub pattern_col_idx: Vec, + /// Estimated multiply–adds of one numeric factorisation, per subtree + /// (used to decide where to fork threads). + pub subtree_work: Vec, + /// Nonzeros in L (including the diagonal) and total work. + pub nnz_l: usize, + pub work: f64, +} + +impl Symbolic { + pub fn num_supernodes(&self) -> usize { + self.start.len() - 1 + } + + pub fn front_rows(&self, s: usize) -> &[usize] { + &self.rows[self.rows_ptr[s]..self.rows_ptr[s + 1]] + } + + pub fn node_children(&self, s: usize) -> &[usize] { + &self.children[self.children_ptr[s]..self.children_ptr[s + 1]] + } + + pub fn max_front(&self) -> usize { + (0..self.num_supernodes()) + .map(|s| self.rows_ptr[s + 1] - self.rows_ptr[s]) + .max() + .unwrap_or(0) + } + + /// Full analysis of an `n × n` CSR pattern (symmetric structure + /// expected; only the lower triangle is factorised). + pub fn analyze(n: usize, row_ptr: &[usize], col_idx: &[usize]) -> Self { + let nd = nested_dissection(n, row_ptr, col_idx); + // Elimination tree under the ND order, then postorder it. + let iperm_nd = inverse(&nd); + let rowpat = lower_rows(n, row_ptr, col_idx, &iperm_nd); + let parent_nd = etree(n, &rowpat); + let post = postorder(&parent_nd); + let perm: Vec = post.iter().map(|&k| nd[k]).collect(); + drop(rowpat); + + let iperm = inverse(&perm); + let rowpat = lower_rows(n, row_ptr, col_idx, &iperm); + let parent = etree(n, &rowpat); + // Column patterns of the strictly lower triangle of PAPᵀ. + let mut colpat: Vec> = vec![Vec::new(); n]; + for (i, cols) in rowpat.iter().enumerate() { + for &j in cols { + colpat[j].push(i); + } + } + drop(rowpat); + let mut nchild = vec![0usize; n]; + let mut child_lists: Vec> = vec![Vec::new(); n]; + for j in 0..n { + if parent[j] != NONE { + nchild[parent[j]] += 1; + child_lists[parent[j]].push(j); + } + } + + // Column structures by merging children (postorder: children of j + // precede j). A structure is kept only while it may still be + // needed: until its parent is processed, and for supernode heads. + let mut structs: Vec>> = vec![None; n]; + let mut counts = vec![0usize; n]; + let mut start = vec![0usize]; + let mut heads_struct: Vec> = Vec::new(); + let mut mark = vec![NONE; n]; + for j in 0..n { + let mut s: Vec = Vec::new(); + mark[j] = j; + for &i in &colpat[j] { + if mark[i] != j { + mark[i] = j; + s.push(i); + } + } + for &c in &child_lists[j] { + if let Some(cs) = structs[c].take() { + for i in cs { + if i != j && mark[i] != j { + mark[i] = j; + s.push(i); + } + } + } + } + s.sort_unstable(); + counts[j] = s.len(); + // Does j extend the previous supernode? + let joins = + j > 0 && parent[j - 1] == j && nchild[j] == 1 && counts[j - 1] == counts[j] + 1; + if !joins { + if j > 0 { + start.push(j); + } + heads_struct.push(s.clone()); + } + structs[j] = Some(s); + } + start.push(n); + drop(structs); + drop(colpat); + let ns = start.len() - 1; + + let mut col_to_sn = vec![0usize; n]; + for s in 0..ns { + for c in start[s]..start[s + 1] { + col_to_sn[c] = s; + } + } + let mut rows_ptr = vec![0usize]; + let mut rows = Vec::new(); + for s in 0..ns { + rows.push(start[s]); + rows.extend_from_slice(&heads_struct[s]); + rows_ptr.push(rows.len()); + } + drop(heads_struct); + + let mut sn_parent = vec![NONE; ns]; + for s in 0..ns { + let last = start[s + 1] - 1; + if parent[last] != NONE { + sn_parent[s] = col_to_sn[parent[last]]; + } + } + let mut children_lists: Vec> = vec![Vec::new(); ns]; + let mut roots = Vec::new(); + for s in 0..ns { + if sn_parent[s] == NONE { + roots.push(s); + } else { + children_lists[sn_parent[s]].push(s); + } + } + let mut children_ptr = vec![0usize]; + let mut children = Vec::new(); + for list in &children_lists { + children.extend_from_slice(list); + children_ptr.push(children.len()); + } + + // Relative indices of each update row in the parent front. + let mut relind_ptr = vec![0usize]; + let mut relind = Vec::new(); + for s in 0..ns { + let k = start[s + 1] - start[s]; + let own = &rows[rows_ptr[s]..rows_ptr[s + 1]]; + if sn_parent[s] != NONE { + let p = sn_parent[s]; + let prow = &rows[rows_ptr[p]..rows_ptr[p + 1]]; + let mut q = 0usize; + for &r in &own[k..] { + while prow[q] != r { + q += 1; + } + relind.push(q as u32); + } + } + relind_ptr.push(relind.len()); + } + + let mut sym = Self { + n, + perm, + start, + rows_ptr, + rows, + children_ptr, + children, + roots, + relind_ptr, + relind, + amap_ptr: Vec::new(), + amap: Vec::new(), + pattern_row_ptr: Vec::new(), + pattern_col_idx: Vec::new(), + subtree_work: Vec::new(), + nnz_l: 0, + work: 0.0, + }; + sym.build_amap(row_ptr, col_idx, &iperm, &col_to_sn); + sym.estimate_work(); + sym + } + + fn build_amap( + &mut self, + row_ptr: &[usize], + col_idx: &[usize], + iperm: &[usize], + col_to_sn: &[usize], + ) { + let ns = self.num_supernodes(); + let mut per_sn: Vec> = vec![Vec::new(); ns]; + for old_row in 0..self.n { + let i = iperm[old_row]; + for idx in row_ptr[old_row]..row_ptr[old_row + 1] { + let j = iperm[col_idx[idx]]; + if i < j { + continue; + } + let s = col_to_sn[j]; + let front = self.front_rows(s); + let m = front.len(); + let a = front + .binary_search(&i) + .expect("A entry outside the symbolic structure"); + let b = j - self.start[s]; + per_sn[s].push((idx as u32, (a + b * m) as u32)); + } + } + self.amap_ptr = vec![0]; + for list in per_sn { + self.amap.extend(list); + self.amap_ptr.push(self.amap.len()); + } + self.pattern_row_ptr = row_ptr.to_vec(); + self.pattern_col_idx = col_idx.to_vec(); + } + + fn estimate_work(&mut self) { + let ns = self.num_supernodes(); + let mut own = vec![0.0f64; ns]; + let mut nnz = 0usize; + let mut total = 0.0; + for (s, slot) in own.iter_mut().enumerate() { + let k = (self.start[s + 1] - self.start[s]) as f64; + let m = (self.rows_ptr[s + 1] - self.rows_ptr[s]) as f64; + let ku = k as usize; + let mu = m as usize; + nnz += ku * mu - ku * (ku - 1) / 2; + // Partial LDLᵀ of k pivots in an m-front (plus the extend-add). + let w = k * (m * m - m * k + k * k / 3.0) + (m - k) * (m - k); + *slot = w; + total += w; + } + // Subtree sums in postorder (children precede parents). + let mut sub = own; + for s in 0..ns { + for &c in self.node_children(s).to_vec().iter() { + sub[s] += sub[c]; + } + } + self.subtree_work = sub; + self.nnz_l = nnz; + self.work = total; + } + + /// Whether `(row_ptr, col_idx)` is exactly the analysed pattern. + pub fn same_pattern(&self, row_ptr: &[usize], col_idx: &[usize]) -> bool { + self.pattern_row_ptr == row_ptr && self.pattern_col_idx == col_idx + } +} + +const NONE: usize = usize::MAX; + +fn inverse(perm: &[usize]) -> Vec { + let mut inv = vec![0usize; perm.len()]; + for (new, &old) in perm.iter().enumerate() { + inv[old] = new; + } + inv +} + +/// Row `i`'s strictly-lower columns `{j < i}` of the permuted, +/// symmetrised pattern. +fn lower_rows(n: usize, row_ptr: &[usize], col_idx: &[usize], iperm: &[usize]) -> Vec> { + let mut rows: Vec> = vec![Vec::new(); n]; + for old_row in 0..n { + let a = iperm[old_row]; + for &old_col in &col_idx[row_ptr[old_row]..row_ptr[old_row + 1]] { + let b = iperm[old_col]; + if a > b { + rows[a].push(b); + } else if b > a { + rows[b].push(a); + } + } + } + for r in &mut rows { + r.sort_unstable(); + r.dedup(); + } + rows +} + +/// Liu's elimination tree with path compression. +fn etree(n: usize, rowpat: &[Vec]) -> Vec { + let mut parent = vec![NONE; n]; + let mut ancestor = vec![NONE; n]; + for i in 0..n { + for &j in &rowpat[i] { + let mut k = j; + while ancestor[k] != NONE && ancestor[k] != i { + let next = ancestor[k]; + ancestor[k] = i; + k = next; + } + if ancestor[k] == NONE { + ancestor[k] = i; + parent[k] = i; + } + } + } + parent +} + +/// Postorder of a forest (`post[k]` = the k-th node visited), children +/// in increasing order, roots in increasing order; iterative DFS. +fn postorder(parent: &[usize]) -> Vec { + let n = parent.len(); + let mut kids: Vec> = vec![Vec::new(); n]; + let mut roots = Vec::new(); + for j in 0..n { + if parent[j] == NONE { + roots.push(j); + } else { + kids[parent[j]].push(j); + } + } + let mut post = Vec::with_capacity(n); + let mut stack: Vec<(usize, usize)> = Vec::new(); + for &r in &roots { + stack.push((r, 0)); + while let Some(&mut (v, ref mut next)) = stack.last_mut() { + if *next < kids[v].len() { + let c = kids[v][*next]; + *next += 1; + stack.push((c, 0)); + } else { + post.push(v); + stack.pop(); + } + } + } + post +} diff --git a/crates/specialized/rtx-fea/tests/sparse_tangent.rs b/crates/specialized/rtx-fea/tests/sparse_tangent.rs new file mode 100644 index 0000000..d23ade5 --- /dev/null +++ b/crates/specialized/rtx-fea/tests/sparse_tangent.rs @@ -0,0 +1,559 @@ +//! R8-g: the sparse Newton-tangent path of the nonlinear Newmark stepper +//! ([`TangentSolver::SparseLdlt`]) against the default banded LU. +//! +//! Suite (fast, run by default): +//! +//! 1. `sparse_matches_banded_csm3_2d` — the 2-D 35×2 Quad8 CSM3 march +//! (the FSI2 harness's flag), 60 steps: every step's displacement +//! within 1e-10 of the banded path's (relative to the peak). +//! 2. `sparse_matches_banded_csm1_3d_plane_strain` — CSM1 static on the +//! 3-D 35×2×1 plane-strain flag. +//! 3. `sparse_matches_banded_csm3_3d_free` — a free-edge 3-D flag +//! (12×2×2), 40 CSM3 steps. +//! 4. `modified_newton_reuses_the_factor` — `reuse = 4`: fewer +//! factorisations than solves, same march to the Newton tolerance. +//! 5. `tangent_knob_parses` — the `RTX_FEA_TANGENT` values. +//! +//! Instruments (`#[ignore]`, env-driven, write under `R8G_OUT`): +//! +//! * `r8g_g1_equivalence` — CSM1 static and a CSM3 march per +//! configuration with both solvers side by side; per-step CSV. +//! * `r8g_g2_cost` — wall time per Newton iteration and per step, per +//! solver mode and configuration, with the sparse path's breakdown. + +use std::io::Write as _; +use std::time::Instant; + +use nalgebra::{DVector, Vector3}; +use rtx_fea::analysis::flag3d::{Flag3d, Flag3dSpec, LateralFaces}; +use rtx_fea::analysis::{ + ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis, TangentSolver, +}; +use rtx_fea::assembly::dof_mapping::DofComponent; +use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType}; +use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction}; +use rtx_fea::elements::{ElementMatrixComputer, StandardFiniteElement}; +use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId}; + +const E_MOD: f64 = 1.4e6; +const NU: f64 = 0.4; +const RHO_CSM: f64 = 1000.0; +const G: f64 = 2.0; + +fn env_str(name: &str, default: &str) -> String { + std::env::var(name).unwrap_or_else(|_| default.to_string()) +} + +fn env_num(name: &str, default: f64) -> f64 { + std::env::var(name) + .map(|v| v.parse().expect(name)) + .unwrap_or(default) +} + +/// The 2-D flag, as the FSI2 harness builds it (copied from +/// `flag3d_structure.rs`). +fn quad8_flag(nx: usize, ny: usize) -> Mesh { + let (x0, x1, y0, y1) = (0.25, 0.6, 0.19, 0.21); + 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 clamp_2d(mesh: &Mesh) -> BoundaryConditionSet { + let clamped: Vec = mesh + .nodes + .iter() + .filter(|(_, node)| (node.position().x - 0.25).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 +} + +/// Consistent gravity nodal forces `∫ N_a ρ g dV`. +fn gravity_forces(mesh: &Mesh, rho: f64, g: f64) -> Vec<(NodeId, Vector3)> { + let dim = mesh.spatial_dimension; + let mut acc: std::collections::BTreeMap> = Default::default(); + for element in mesh.elements.values() { + let coords: Vec> = element + .nodes + .iter() + .map(|id| mesh.get_node(*id).unwrap().position()) + .collect(); + let fe = StandardFiniteElement::new(element.element_type, coords.clone()); + let f = ElementMatrixComputer::compute_body_force_vector( + &fe, + &coords, + &|_| Vector3::new(0.0, -rho * g, 0.0), + None, + ) + .unwrap(); + for (a, id) in element.nodes.iter().enumerate() { + let e = acc.entry(*id).or_insert_with(Vector3::zeros); + for c in 0..dim { + e[c] += f[a * dim + c]; + } + } + } + acc.into_iter().collect() +} + +fn static_criteria() -> ConvergenceCriteria { + ConvergenceCriteria { + force_tolerance: 1e-12, + displacement_tolerance: 1e-14, + max_iterations: 60, + ..ConvergenceCriteria::default() + } +} + +/// The CSM3 analysis (gravity from rest) on the 2-D 35×2 flag. +fn csm3_2d(dt: f64, solver: TangentSolver) -> NonlinearDynamicAnalysis { + let mesh = quad8_flag(35, 2); + let mut analysis = NonlinearDynamicAnalysis::new( + mesh.clone(), + Flag3d::materials(E_MOD, NU, RHO_CSM), + clamp_2d(&mesh), + dt, + 1, + Default::default(), + ) + .with_total_lagrangian() + .with_tangent_solver(solver); + analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0)); + analysis +} + +/// The CSM3 analysis on a 3-D flag (R8-b's instrument settings). +fn csm3_3d( + flag: &Flag3d, + lateral: LateralFaces, + dt: f64, + solver: TangentSolver, +) -> NonlinearDynamicAnalysis { + let mut analysis = flag + .dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, dt, 1, 0.5) + .with_tangent_solver(solver); + analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0)); + analysis +} + +/// `max_i |a_i − b_i| / max_i |b_i|`. +fn rel_max(a: &DVector, b: &DVector) -> f64 { + let scale = b.amax().max(1e-300); + (a - b).amax() / scale +} + +/// March both analyses `steps` steps side by side; per step the +/// relative displacement difference. Returns (max rel, per-step rows). +fn march_pair( + banded: &NonlinearDynamicAnalysis, + sparse: &NonlinearDynamicAnalysis, + steps: usize, +) -> (f64, Vec<(f64, f64, usize, usize)>) { + let mut sb = banded.stepper().unwrap(); + let mut ss = sparse.stepper().unwrap(); + let mut stb = sb.rest_state().unwrap(); + let mut sts = ss.rest_state().unwrap(); + let mut worst: f64 = rel_max(&sts.acceleration, &stb.acceleration); + let mut rows = Vec::new(); + for _ in 0..steps { + let t0 = Instant::now(); + let (nb, ib) = sb.step(&stb).unwrap(); + let tb = t0.elapsed().as_secs_f64(); + let t1 = Instant::now(); + let (ns, is) = ss.step(&sts).unwrap(); + let ts = t1.elapsed().as_secs_f64(); + stb = nb; + sts = ns; + let r = rel_max(&sts.displacement, &stb.displacement); + worst = worst.max(r); + rows.push((r, tb / ts.max(1e-12), ib, is)); + } + (worst, rows) +} + +#[test] +fn sparse_matches_banded_csm3_2d() { + let dt = 0.005; + let (worst, rows) = march_pair( + &csm3_2d(dt, TangentSolver::BandedLu), + &csm3_2d(dt, TangentSolver::SPARSE), + 60, + ); + let same_newton = rows.iter().all(|r| r.2 == r.3); + println!("CSM3 2-D 35x2, 60 steps: max rel |Δu| {worst:.3e}, same Newton counts {same_newton}"); + assert!(worst < 1e-10, "sparse departs from banded: {worst:.3e}"); + assert!(same_newton); +} + +#[test] +fn sparse_matches_banded_csm1_3d_plane_strain() { + let flag = Flag3d::build(Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1)).unwrap(); + let forces = gravity_forces(&flag.mesh, RHO_CSM, G); + let solve = |solver: TangentSolver| { + let analysis = flag + .dynamic_analysis(E_MOD, NU, RHO_CSM, LateralFaces::PlaneStrain, 1e4, 1, 0.5) + .with_convergence_criteria(static_criteria()) + .with_tangent_solver(solver); + let mut stepper = analysis.stepper().unwrap(); + let n = stepper.rest_state().unwrap().displacement.len(); + let mut u = DVector::zeros(n); + for s in 1..=5 { + let scale = s as f64 / 5.0; + let scaled: Vec<_> = forces.iter().map(|(id, f)| (*id, f * scale)).collect(); + stepper.set_nodal_forces(&scaled); + let state = DynamicState { + displacement: u.clone(), + velocity: DVector::zeros(n), + acceleration: DVector::zeros(n), + }; + u = stepper.step(&state).unwrap().0.displacement; + } + let d = stepper.node_dofs(flag.point_a()); + (u.clone(), u[d[1]]) + }; + let (ub, ay_b) = solve(TangentSolver::BandedLu); + let (us, ay_s) = solve(TangentSolver::SPARSE); + let r = rel_max(&us, &ub); + println!( + "CSM1 3-D 35x2x1 ps: uy(A) banded {ay_b:.12e} sparse {ay_s:.12e}, max rel |Δu| {r:.3e}" + ); + assert!(r < 1e-10, "CSM1 departs: {r:.3e}"); + assert!((ay_b + 65.1406e-3).abs() < 1e-6, "CSM1 uy(A) moved: {ay_b}"); +} + +#[test] +fn sparse_matches_banded_csm3_3d_free() { + let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 12, 2, 2)).unwrap(); + let dt = 0.005; + let (worst, _) = march_pair( + &csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::BandedLu), + &csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::SPARSE), + 40, + ); + println!("CSM3 3-D 12x2x2 free, 40 steps: max rel |Δu| {worst:.3e}"); + assert!(worst < 1e-10, "sparse departs from banded: {worst:.3e}"); +} + +#[test] +fn modified_newton_reuses_the_factor() { + let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 12, 2, 2)).unwrap(); + let dt = 0.005; + let full = csm3_3d(&flag, LateralFaces::Free, dt, TangentSolver::SPARSE); + let modified = csm3_3d( + &flag, + LateralFaces::Free, + dt, + TangentSolver::SparseLdlt { reuse: 4 }, + ); + let mut sf = full.stepper().unwrap(); + let mut sm = modified.stepper().unwrap(); + let mut stf = sf.rest_state().unwrap(); + let mut stm = sm.rest_state().unwrap(); + for _ in 0..40 { + stf = sf.step(&stf).unwrap().0; + stm = sm.step(&stm).unwrap().0; + } + let r = rel_max(&stm.displacement, &stf.displacement); + let stats_f = sf.tangent_stats().unwrap(); + let stats_m = sm.tangent_stats().unwrap(); + println!( + "modified Newton (reuse 4) vs full after 40 steps: rel {r:.3e}; full {stats_f:?}; modified {stats_m:?}" + ); + assert_eq!(stats_f.factorizations, stats_f.solves); + assert!(stats_m.factorizations < stats_m.solves); + assert_eq!(stats_m.analyses, 1); + // Newton's own displacement tolerance is 1e-6 (relative). + assert!(r < 1e-4, "modified Newton drifts: {r:.3e}"); +} + +#[test] +fn tangent_knob_parses() { + // Only the parser (the knob is read when a stepper is built). + assert_eq!(TangentSolver::default(), TangentSolver::BandedLu); + assert_eq!( + TangentSolver::SPARSE, + TangentSolver::SparseLdlt { reuse: 1 } + ); +} + +// --------------------------------------------------------------------------- +// Instruments +// --------------------------------------------------------------------------- + +/// `NXxNYxNZ:span:free|ps` entries (as in `flag3d_structure.rs`); `2d` +/// = the 2-D 35×2 Quad8 flag. +fn parse_configs(spec: &str) -> Vec> { + spec.split(',') + .map(|entry| { + if entry.trim() == "2d" { + return None; + } + let mut parts = entry.trim().split(':'); + let mesh = parts.next().unwrap(); + let span: f64 = parts.next().unwrap().parse().unwrap(); + let lateral = match parts.next().unwrap() { + "free" => LateralFaces::Free, + "ps" => LateralFaces::PlaneStrain, + other => panic!("lateral {other}"), + }; + let n: Vec = mesh.split('x').map(|t| t.parse().unwrap()).collect(); + Some((n[0], n[1], n[2], span, lateral)) + }) + .collect() +} + +fn tag(cfg: &Option<(usize, usize, usize, f64, LateralFaces)>) -> String { + match cfg { + None => "2d_35x2".to_string(), + Some((nx, ny, nz, span, lateral)) => { + let l = if *lateral == LateralFaces::Free { + "free" + } else { + "ps" + }; + format!("{nx}x{ny}x{nz}_s{span}_{l}") + } + } +} + +fn build_flag(cfg: &(usize, usize, usize, f64, LateralFaces)) -> Flag3d { + let (nx, ny, nz, span, _) = *cfg; + Flag3d::build(Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz)).unwrap() +} + +fn analysis_for( + cfg: &Option<(usize, usize, usize, f64, LateralFaces)>, + flag: Option<&Flag3d>, + dt: f64, + solver: TangentSolver, +) -> NonlinearDynamicAnalysis { + match cfg { + None => csm3_2d(dt, solver), + Some(c) => csm3_3d(flag.unwrap(), c.4, dt, solver), + } +} + +#[test] +#[ignore = "instrument: G1 — sparse vs banded, CSM1 static and a CSM3 march"] +fn r8g_g1_equivalence() { + let out = env_str("R8G_OUT", "."); + let configs = parse_configs(&env_str("R8G_CONFIGS", "35x2x1:0.05:ps,35x2x8:0.41:free")); + let steps = env_num("R8G_STEPS", 200.0) as usize; + let dt = env_num("R8G_DT", 0.005); + let do_csm1 = env_str("R8G_CSM1", "1") == "1"; + let mut table = std::fs::OpenOptions::new() + .create(true) + .append(true) + .open(format!("{out}/g1_table.txt")) + .unwrap(); + for cfg in &configs { + let name = tag(cfg); + let flag = cfg.as_ref().map(build_flag); + if do_csm1 { + if let (Some(c), Some(flag)) = (cfg, flag.as_ref()) { + let forces = gravity_forces(&flag.mesh, RHO_CSM, G); + let mut results = Vec::new(); + for solver in [TangentSolver::BandedLu, TangentSolver::SPARSE] { + let start = Instant::now(); + let analysis = flag + .dynamic_analysis(E_MOD, NU, RHO_CSM, c.4, 1e4, 1, 0.5) + .with_convergence_criteria(static_criteria()) + .with_tangent_solver(solver); + let mut stepper = analysis.stepper().unwrap(); + let n = stepper.rest_state().unwrap().displacement.len(); + let mut u = DVector::zeros(n); + let mut newton = 0; + for s in 1..=5 { + let scale = s as f64 / 5.0; + let scaled: Vec<_> = + forces.iter().map(|(id, f)| (*id, f * scale)).collect(); + stepper.set_nodal_forces(&scaled); + let state = DynamicState { + displacement: u.clone(), + velocity: DVector::zeros(n), + acceleration: DVector::zeros(n), + }; + let (next, it) = stepper.step(&state).unwrap(); + newton += it; + u = next.displacement; + } + let d = stepper.node_dofs(flag.point_a()); + results.push(( + u.clone(), + u[d[0]], + u[d[1]], + newton, + start.elapsed().as_secs_f64(), + )); + } + let r = rel_max(&results[1].0, &results[0].0); + let line = format!( + "G1 CSM1 {name}: A banded ux {:.9e} uy {:.9e} [{} Newton, {:.1} s] | sparse ux \ + {:.9e} uy {:.9e} [{} Newton, {:.1} s] | max rel |Δu| {r:.3e}", + results[0].1, + results[0].2, + results[0].3, + results[0].4, + results[1].1, + results[1].2, + results[1].3, + results[1].4 + ); + println!("{line}"); + writeln!(table, "{line}").unwrap(); + } + } + let banded = analysis_for(cfg, flag.as_ref(), dt, TangentSolver::BandedLu); + let sparse = analysis_for(cfg, flag.as_ref(), dt, TangentSolver::SPARSE); + let start = Instant::now(); + let (worst, rows) = march_pair(&banded, &sparse, steps); + let path = format!("{out}/g1_csm3_{name}_dt{dt}.csv"); + let mut csv = std::fs::File::create(&path).unwrap(); + writeln!(csv, "step,rel_u,speedup,newton_banded,newton_sparse").unwrap(); + for (k, (r, speed, ib, is)) in rows.iter().enumerate() { + writeln!(csv, "{},{r:.6e},{speed:.3},{ib},{is}", k + 1).unwrap(); + } + let differ = rows.iter().filter(|r| r.2 != r.3).count(); + let line = format!( + "G1 CSM3 {name} dt {dt}: {steps} steps, max rel |Δu| {worst:.3e}, Newton-count \ + differences {differ}, {:.0} s → {path}", + start.elapsed().as_secs_f64() + ); + println!("{line}"); + writeln!(table, "{line}").unwrap(); + } +} + +#[test] +#[ignore = "instrument: G2 — cost per Newton iteration and per step"] +fn r8g_g2_cost() { + let out = env_str("R8G_OUT", "."); + let configs = parse_configs(&env_str( + "R8G_CONFIGS", + "2d,35x2x8:0.41:free,35x2x16:0.41:free", + )); + let modes: Vec = env_str("R8G_MODES", "banded,sparse,sparse:3") + .split(',') + .map(str::to_string) + .collect(); + let steps = env_num("R8G_STEPS", 40.0) as usize; + let warm = env_num("R8G_WARM", 20.0) as usize; + let dt = env_num("R8G_DT", 0.005); + let mut table = std::fs::OpenOptions::new() + .create(true) + .append(true) + .open(format!("{out}/g2_table.txt")) + .unwrap(); + for cfg in &configs { + let name = tag(cfg); + let flag = cfg.as_ref().map(build_flag); + let mut reference: Option> = None; + for mode in &modes { + let solver = match mode.as_str() { + "banded" => TangentSolver::BandedLu, + "sparse" => TangentSolver::SPARSE, + other => TangentSolver::SparseLdlt { + reuse: other.strip_prefix("sparse:").unwrap().parse().unwrap(), + }, + }; + let analysis = analysis_for(cfg, flag.as_ref(), dt, solver); + let t_build = Instant::now(); + let mut stepper = analysis.stepper().unwrap(); + let mut state = stepper.rest_state().unwrap(); + let build = t_build.elapsed().as_secs_f64(); + // Warm-up steps (off the clock) move the flag off rest. + for _ in 0..warm { + state = stepper.step(&state).unwrap().0; + } + let before = stepper.tangent_stats(); + let mut times = Vec::with_capacity(steps); + let mut newton = 0usize; + for _ in 0..steps { + let t = Instant::now(); + let (next, it) = stepper.step(&state).unwrap(); + times.push(t.elapsed().as_secs_f64()); + newton += it; + state = next; + } + let total: f64 = times.iter().sum(); + let mut sorted = times.clone(); + sorted.sort_by(f64::total_cmp); + let median = sorted[sorted.len() / 2]; + let drift = reference + .as_ref() + .map_or(f64::NAN, |r| rel_max(&state.displacement, r)); + if reference.is_none() { + reference = Some(state.displacement.clone()); + } + let breakdown = match (before, stepper.tangent_stats()) { + (Some(b), Some(a)) => format!( + " | assemblies {} ({:.4} s each), factorisations {} ({:.4} s each incl. \ + scatter), solves {} ({:.4} s each), fallbacks {}, nnz(L) {}", + a.assemblies - b.assemblies, + (a.assembly_seconds - b.assembly_seconds) + / (a.assemblies - b.assemblies).max(1) as f64, + a.factorizations - b.factorizations, + (a.factor_seconds - b.factor_seconds) + / (a.factorizations - b.factorizations).max(1) as f64, + a.solves - b.solves, + (a.solve_seconds - b.solve_seconds) / (a.solves - b.solves).max(1) as f64, + a.fallbacks, + a.nnz_l + ), + _ => String::new(), + }; + let line = format!( + "G2 {name} {mode} dt {dt}: dofs {} | build+rest {build:.3} s | {steps} steps after \ + {warm} warm: {:.4} s/step mean, {median:.4} median, {:.2} Newton/step, {:.4} \ + s/Newton | final rel to first mode {drift:.2e} | threads {}{breakdown}", + state.displacement.len(), + total / steps as f64, + newton as f64 / steps as f64, + total / newton.max(1) as f64, + rayon::current_num_threads() + ); + println!("{line}"); + writeln!(table, "{line}").unwrap(); + } + } +} diff --git a/crates/specialized/rtx-fea/tests/sparse_tangent_bench.rs b/crates/specialized/rtx-fea/tests/sparse_tangent_bench.rs new file mode 100644 index 0000000..37cbcf6 --- /dev/null +++ b/crates/specialized/rtx-fea/tests/sparse_tangent_bench.rs @@ -0,0 +1,89 @@ +//! R8-g instrument: the sparse LDLᵀ against the banded LU on the dumped +//! 3-D flag operators (R8-b's `k_*.coo` / `m_*.coo`, the TL tangent at +//! u = 0 and the consistent mass on the free DOFs). The Newmark tangent +//! `K + M/(β Δt²)` at CSM3's Δt = 0.005 is solved by both; the timing and +//! the agreement are printed. +//! +//! `R8G_COO_DIR` (default the R8-b study dir), `R8G_TAGS` (comma list). + +use nalgebra::DVector; +use rtx_fea::assembly::SparseMatrix; +use rtx_fea::solvers::{BandedLu, LinearSolver, SolverOptions, SparseLdlt}; +use std::time::Instant; + +fn read_coo(path: &str) -> Vec<(usize, usize, f64)> { + let text = std::fs::read_to_string(path).unwrap_or_else(|e| panic!("{path}: {e}")); + text.lines() + .filter(|l| !l.trim().is_empty()) + .map(|l| { + let mut it = l.split_whitespace(); + let i: usize = it.next().unwrap().parse().unwrap(); + let j: usize = it.next().unwrap().parse().unwrap(); + let v: f64 = it.next().unwrap().parse().unwrap(); + (i, j, v) + }) + .collect() +} + +#[test] +#[ignore = "instrument: sparse LDLt vs banded LU on the dumped flag tangents"] +fn ldlt_vs_banded_on_dumped_tangents() { + let dir = std::env::var("R8G_COO_DIR").unwrap_or_else(|_| { + "/home/osobh/projects/omni-cortex-data/fsi_studies/p1_regress/r8b".to_string() + }); + let tags = std::env::var("R8G_TAGS") + .unwrap_or_else(|_| "35x2x1_s0.05_ps,35x2x8_s0.41_free".to_string()); + let skip_banded = std::env::var("R8G_SKIP_BANDED").is_ok(); + let coef = 1.0 / (0.25 * 0.005 * 0.005); + for tag in tags.split(',') { + let k = read_coo(&format!("{dir}/k_{tag}.coo")); + let m = read_coo(&format!("{dir}/m_{tag}.coo")); + let n = k.iter().map(|t| t.0.max(t.1)).max().unwrap() + 1; + let mut triplets = k.clone(); + triplets.extend(m.iter().map(|&(i, j, v)| (i, j, coef * v))); + let a = SparseMatrix::from_triplets(n, n, &triplets).unwrap(); + let x_exact = DVector::from_fn(n, |i, _| ((i as f64) * 0.37).sin() + 0.2); + let b = a.multiply_vector(&x_exact).unwrap(); + let opts = SolverOptions::default(); + + let mut ldlt = SparseLdlt::new(); + let t0 = Instant::now(); + let (row_ptr, col_idx) = a.structure(); + ldlt.analyze(n, row_ptr, col_idx).unwrap(); + let t_sym = t0.elapsed().as_secs_f64(); + let mut t_num = f64::MAX; + let mut t_solve = f64::MAX; + let mut x = DVector::zeros(n); + for _ in 0..5 { + let t1 = Instant::now(); + ldlt.factorize_values(a.values()).unwrap(); + t_num = t_num.min(t1.elapsed().as_secs_f64()); + let t2 = Instant::now(); + x = ldlt.solve_factored(&b).unwrap(); + t_solve = t_solve.min(t2.elapsed().as_secs_f64()); + } + let stats = ldlt.stats().unwrap(); + let rel_exact = (&x - &x_exact).norm() / x_exact.norm(); + let residual = (a.multiply_vector(&x).unwrap() - &b).norm() / b.norm(); + let (t_band, rel_band) = if skip_banded { + (f64::NAN, f64::NAN) + } else { + let t3 = Instant::now(); + let (xb, _) = BandedLu::new().solve(&a, &b, &opts).unwrap(); + let t = t3.elapsed().as_secs_f64(); + (t, (&x - &xb).norm() / xb.norm()) + }; + println!( + "{tag}: n {n} nnz(A) {} | LDLt symbolic {t_sym:.3} s, numeric {t_num:.4} s, solve \ + {t_solve:.4} s | nnz(L) {} supernodes {} max front {} work {:.3e} | banded LU \ + {t_band:.3} s | rel vs exact {rel_exact:.2e}, vs banded {rel_band:.2e}, residual \ + {residual:.2e} | threads {}", + a.nnz(), + stats.nnz_l, + stats.supernodes, + stats.max_front, + stats.work, + rayon::current_num_threads() + ); + } +} From d50a0e28d2996db7adf860783e7166751687cc31 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 25 Sep 2026 22:52:23 -0500 Subject: [PATCH 2/2] R8-g: parallel scatter + trailing updates + tree solves; force-only TL kernel for modified Newton - the tangent scatter runs per CSR entry over a transposed contribution map, in the serial scatter's order (same bits), on rayon; - large fronts' trailing gemms run per column block on rayon; - forward solve multifrontal, backward top-down, subtrees on rayon; - total_lagrangian::internal_force: the force of internal_force_and_tangent alone (bit-identical, tested), used by modified-Newton iterations that reuse the factor; the tangent is evaluated only when a refresh is due. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../rtx-fea/src/analysis/nonlinear_dynamic.rs | 78 +++++++- .../rtx-fea/src/analysis/sparse_tangent.rs | 148 ++++++++++---- .../rtx-fea/src/elements/total_lagrangian.rs | 93 +++++++++ .../src/solvers/sparse_ldlt/numeric.rs | 180 ++++++++++++++---- .../rtx-fea/tests/sparse_tangent.rs | 39 +++- 5 files changed, 459 insertions(+), 79 deletions(-) diff --git a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs index bc8735f..ffecfd0 100644 --- a/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs +++ b/crates/specialized/rtx-fea/src/analysis/nonlinear_dynamic.rs @@ -650,8 +650,9 @@ impl<'a> NonlinearDynamicStepper<'a> { for &dof in &self.free_dofs { a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]); } - let (f_int, tangent) = if self.sparse.is_some() { - (self.sparse_forces(&u_iter)?, None) + let (f_int, tangent) = if let Some(sparse) = &self.sparse { + let with_tangent = sparse.wants_tangent(line_search, inv_beta_dt2); + (self.sparse_forces(&u_iter, with_tangent)?, None) } else { let (f_int, tangent) = self.assemble(&u_iter, true, inv_beta_dt2)?; (f_int, Some(tangent)) @@ -676,9 +677,16 @@ impl<'a> NonlinearDynamicStepper<'a> { None => { let rate = residual_norm / previous_norm; previous_norm = residual_norm; + let sparse = self.sparse.as_ref().expect("sparse path"); + if sparse.needs_refresh(line_search, inv_beta_dt2, rate) + && !sparse.tangents_current() + { + self.sparse_tangents(&u_iter)?; + } + let masses: Vec<&DMatrix> = self.caches.iter().map(|c| &c.mass).collect(); let sparse = self.sparse.as_mut().expect("sparse path"); sparse.solve( - self.caches.iter().map(|c| &c.mass), + &masses, &self.added_mass, &self.free_index, &residual, @@ -809,13 +817,26 @@ impl<'a> NonlinearDynamicStepper<'a> { /// forces summed in element order (the banded path's order, so the /// internal force is the same float for float), the tangents kept for /// the factorisation. - fn sparse_forces(&mut self, solution: &DVector) -> FeaResult> { + /// With `with_tangent` false (a modified-Newton iteration that will + /// reuse the factor) only the forces are evaluated — the TL force-only + /// kernel gives the same forces bit for bit. + fn sparse_forces( + &mut self, + solution: &DVector, + with_tangent: bool, + ) -> FeaResult> { let started = std::time::Instant::now(); let analysis = self.analysis; - let results: Vec<(DVector, DMatrix)> = self + let results: Vec<(DVector, Option>)> = self .caches .par_iter() - .map(|cache| element_force_and_tangent(analysis, cache, solution)) + .map(|cache| { + if with_tangent { + element_force_and_tangent(analysis, cache, solution).map(|(f, k)| (f, Some(k))) + } else { + element_force(analysis, cache, solution).map(|f| (f, None)) + } + }) .collect::>()?; let sparse = self.sparse.as_mut().expect("sparse path"); Ok(sparse.accumulate( @@ -826,6 +847,23 @@ impl<'a> NonlinearDynamicStepper<'a> { )) } + /// Sparse path: the element tangents at `solution` alone (a refresh + /// after a force-only assembly). + fn sparse_tangents(&mut self, solution: &DVector) -> FeaResult<()> { + let started = std::time::Instant::now(); + let analysis = self.analysis; + let tangents: Vec> = self + .caches + .par_iter() + .map(|cache| element_force_and_tangent(analysis, cache, solution).map(|(_, k)| k)) + .collect::>()?; + self.sparse + .as_mut() + .expect("sparse path") + .set_tangents(tangents, started); + Ok(()) + } + /// Internal force and (optionally) tangent at a full displacement /// vector, reduced to the free DOFs. The tangent includes the Newmark /// mass term `M / (β Δt²)` at the caller's `inv_beta_dt2` (the @@ -935,3 +973,31 @@ fn element_force_and_tangent( ) } } + +/// One element's internal force alone (the TL force-only kernel; the +/// small-strain path evaluates the pair and drops the tangent). +fn element_force( + analysis: &NonlinearDynamicAnalysis, + cache: &ElementCache, + solution: &DVector, +) -> FeaResult> { + if !analysis.total_lagrangian { + return element_force_and_tangent(analysis, cache, solution).map(|(f, _)| f); + } + let dim = analysis.mesh.spatial_dimension; + let material = analysis.materials.get_material(cache.material_id).unwrap(); + let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone()); + let mut element_displacement = DVector::zeros(cache.dofs.len()); + for (local, &dof) in cache.dofs.iter().enumerate() { + element_displacement[local] = solution[dof]; + } + let (lambda, mu) = material.properties().lame_parameters(); + let constitutive = saint_venant_kirchhoff(lambda, mu, dim); + total_lagrangian::internal_force( + &fe, + &cache.coords, + &element_displacement, + constitutive.as_ref(), + None, + ) +} diff --git a/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs b/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs index 78f462c..b23df05 100644 --- a/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs +++ b/crates/specialized/rtx-fea/src/analysis/sparse_tangent.rs @@ -22,8 +22,12 @@ use crate::assembly::SparseMatrix; use crate::error::FeaResult; use crate::solvers::{BandedLu, LinearSolver, SolverOptions, SparseLdlt}; use nalgebra::{DMatrix, DVector}; +use rayon::prelude::*; use std::time::Instant; +/// CSR entries per parallel scatter task. +const SCATTER_CHUNK: usize = 16_384; + /// Which linear solver carries the Newton tangent. #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] pub enum TangentSolver { @@ -71,8 +75,12 @@ impl TangentSolver { /// Counters and accumulated wall time of the sparse path. #[derive(Debug, Clone, Copy, Default, PartialEq)] pub struct TangentStats { - /// Element force/tangent evaluations (parallel) + force summation. + /// Element evaluations (parallel) + force summation; with the + /// tangents under full Newton, forces only on reused iterations. pub assemblies: usize, + /// Extra tangent evaluations for a refresh decided after a + /// force-only assembly (slow contraction). + pub tangent_refreshes: usize, /// Numeric factorisations. pub factorizations: usize, /// Linear solves (Newton corrections). @@ -83,7 +91,9 @@ pub struct TangentStats { pub fallbacks: usize, /// Seconds in element evaluation + force summation. pub assembly_seconds: f64, - /// Seconds in tangent scatter + numeric factorisation. + /// Seconds in the tangent scatter (element blocks → CSR). + pub scatter_seconds: f64, + /// Seconds in the numeric factorisation. pub factor_seconds: f64, /// Seconds in triangular solves. pub solve_seconds: f64, @@ -97,13 +107,18 @@ pub(super) struct SparseTangent { row_ptr: Vec, col_idx: Vec, values: Vec, - /// Per element, `(local_row · n_e + local_col)` → CSR index, or - /// `u32::MAX` when either DOF is constrained. - elem_pos: Vec>, + /// Per CSR entry, its element contributions in element order: + /// `contrib[contrib_ptr[p]..contrib_ptr[p + 1]]` = `(element, + /// column-major local index)` — the transposed scatter, so entries can + /// be summed in parallel in exactly the serial scatter's order. + contrib_ptr: Vec, + contrib: Vec<(u32, u32)>, /// CSR index of each free DOF's diagonal. diag_pos: Vec, - /// The element tangents of the latest assembly. + /// The element tangents of the latest assembly that computed them. element_tangents: Vec>, + /// Whether `element_tangents` belong to the latest assembly's state. + tangents_current: bool, solver: SparseLdlt, /// Separate instance for the rest state's mass solve. pub(super) mass_solver: SparseLdlt, @@ -155,19 +170,32 @@ impl SparseTangent { .binary_search(&c) .expect("pattern covers every element pair") }; - let mut elem_pos = Vec::new(); - for dofs in elements { + // Element-major, row-major within the element: the order in + // which a serial scatter would add the contributions. + let mut flat: Vec<(usize, u32, u32)> = Vec::new(); + for (e, dofs) in elements.enumerate() { let n_e = dofs.len(); - let mut pos = vec![u32::MAX; n_e * n_e]; for (a, &r) in dofs.iter().enumerate() { let Some(fr) = free_index[r] else { continue }; for (b, &c) in dofs.iter().enumerate() { if let Some(fc) = free_index[c] { - pos[a * n_e + b] = find(fr, fc) as u32; + flat.push((find(fr, fc), e as u32, (b * n_e + a) as u32)); } } } - elem_pos.push(pos); + } + let mut contrib_ptr = vec![0usize; col_idx.len() + 1]; + for &(p, _, _) in &flat { + contrib_ptr[p + 1] += 1; + } + for p in 0..col_idx.len() { + contrib_ptr[p + 1] += contrib_ptr[p]; + } + let mut fill = contrib_ptr.clone(); + let mut contrib = vec![(0u32, 0u32); flat.len()]; + for (p, e, l) in flat { + contrib[fill[p]] = (e, l); + fill[p] += 1; } let diag_pos = (0..num_free).map(|d| find(d, d)).collect(); let mut solver = SparseLdlt::new(); @@ -183,9 +211,11 @@ impl SparseTangent { row_ptr, col_idx, values, - elem_pos, + contrib_ptr, + contrib, diag_pos, element_tangents: Vec::new(), + tangents_current: false, solver, mass_solver: SparseLdlt::new(), reuse: reuse.max(1), @@ -201,38 +231,73 @@ impl SparseTangent { self.stats } + /// Whether the next Newton solve will factorise for sure (so the + /// assembly should produce the element tangents): always under full + /// Newton; under modified Newton when the factor is missing, used up, + /// for another dt, or `force_fresh`. + pub(super) fn wants_tangent(&self, force_fresh: bool, coef: f64) -> bool { + #[allow(clippy::float_cmp)] + let same_coef = coef == self.factor_coef; + self.reuse == 1 || force_fresh || !self.factor_valid || self.uses_left == 0 || !same_coef + } + + /// Whether the solve at contraction `rate` will factorise. + pub(super) fn needs_refresh(&self, force_fresh: bool, coef: f64, rate: f64) -> bool { + self.wants_tangent(force_fresh, coef) || rate > 0.5 + } + + /// Whether the stored element tangents are those of the latest state. + pub(super) fn tangents_current(&self) -> bool { + self.tangents_current + } + /// Sum the element internal forces (element order, local order — the - /// banded path's order) into the free-DOF vector and keep the element - /// tangents for a later factorisation. + /// banded path's order) into the free-DOF vector, and keep the element + /// tangents when they were computed. pub(super) fn accumulate( &mut self, - element_results: Vec<(DVector, DMatrix)>, + element_results: Vec<(DVector, Option>)>, element_dofs: impl Iterator>, free_index: &[Option], started: Instant, ) -> DVector { let mut internal = DVector::zeros(self.num_free); - self.element_tangents.clear(); + let with_tangents = element_results.first().is_some_and(|r| r.1.is_some()); + if with_tangents { + self.element_tangents.clear(); + } for ((f_e, k_e), dofs) in element_results.into_iter().zip(element_dofs) { for (local, &dof) in dofs.as_ref().iter().enumerate() { if let Some(free) = free_index[dof] { internal[free] += f_e[local]; } } - self.element_tangents.push(k_e); + if let Some(k_e) = k_e { + self.element_tangents.push(k_e); + } } + self.tangents_current = with_tangents; self.stats.assemblies += 1; self.stats.assembly_seconds += started.elapsed().as_secs_f64(); internal } + /// Install element tangents computed separately for the latest state + /// (a modified-Newton refresh the force-only assembly did not foresee). + pub(super) fn set_tangents(&mut self, tangents: Vec>, started: Instant) { + self.element_tangents = tangents; + self.tangents_current = true; + self.stats.tangent_refreshes += 1; + self.stats.assembly_seconds += started.elapsed().as_secs_f64(); + } + /// The Newton correction `(K_T + coef·(M + M_added)) δ = r`, with the /// tangent of the latest [`Self::accumulate`]. `force_fresh` demands a /// new factor; `rate` is `‖r_i‖/‖r_{i−1}‖` within the step (0 at the /// first iteration). - pub(super) fn solve<'m>( + pub(super) fn solve( &mut self, - masses: impl Iterator>, + masses: &[&DMatrix], added_mass: &DVector, free_index: &[Option], residual: &DVector, @@ -240,24 +305,33 @@ impl SparseTangent { force_fresh: bool, rate: f64, ) -> FeaResult> { - #[allow(clippy::float_cmp)] - let same_coef = coef == self.factor_coef; - let refresh = - force_fresh || !self.factor_valid || self.uses_left == 0 || !same_coef || rate > 0.5; - if refresh { - let started = Instant::now(); - self.values.fill(0.0); - for ((k_e, mass), pos) in self.element_tangents.iter().zip(masses).zip(&self.elem_pos) { - let n_e = k_e.nrows(); - for a in 0..n_e { - for b in 0..n_e { - let p = pos[a * n_e + b]; - if p != u32::MAX { - self.values[p as usize] += k_e[(a, b)] + coef * mass[(a, b)]; - } - } + if self.needs_refresh(force_fresh, coef, rate) { + if !self.tangents_current { + return Err(crate::error::SolverError::FactorizationFailed { + reason: "sparse tangent: refresh without the current element tangents" + .to_string(), } + .into()); } + let started = Instant::now(); + let tangents = &self.element_tangents; + let contrib_ptr = &self.contrib_ptr; + let contrib = &self.contrib; + self.values + .par_chunks_mut(SCATTER_CHUNK) + .enumerate() + .for_each(|(chunk, values)| { + let base = chunk * SCATTER_CHUNK; + for (offset, value) in values.iter_mut().enumerate() { + let p = base + offset; + let mut sum = 0.0; + for &(e, l) in &contrib[contrib_ptr[p]..contrib_ptr[p + 1]] { + let (e, l) = (e as usize, l as usize); + sum += tangents[e].as_slice()[l] + coef * masses[e].as_slice()[l]; + } + *value = sum; + } + }); for (dof, &m) in added_mass.iter().enumerate() { if m != 0.0 { if let Some(free) = free_index[dof] { @@ -265,8 +339,10 @@ impl SparseTangent { } } } + let scattered = Instant::now(); + self.stats.scatter_seconds += (scattered - started).as_secs_f64(); let outcome = self.solver.factorize_values(&self.values); - self.stats.factor_seconds += started.elapsed().as_secs_f64(); + self.stats.factor_seconds += scattered.elapsed().as_secs_f64(); if let Err(error) = outcome { // A failed pivot: this solve goes to the pivoting banded // LU on the same matrix; the next one tries LDLᵀ again. diff --git a/crates/specialized/rtx-fea/src/elements/total_lagrangian.rs b/crates/specialized/rtx-fea/src/elements/total_lagrangian.rs index 4d6616b..3909c07 100644 --- a/crates/specialized/rtx-fea/src/elements/total_lagrangian.rs +++ b/crates/specialized/rtx-fea/src/elements/total_lagrangian.rs @@ -194,3 +194,96 @@ pub fn internal_force_and_tangent( Ok((internal_force, tangent)) } + +/// The internal force alone of [`internal_force_and_tangent`] — the same +/// operations in the same order for `f_int` (so the result is identical +/// bit for bit), without the tangent's `B_Lᵀ C B_L` and geometric terms. +/// For Newton iterations that reuse a factorised tangent (modified +/// Newton, R8-g). +pub fn internal_force( + element: &dyn FiniteElement, + node_coords: &[Vector3], + element_displacement: &DVector, + constitutive: &dyn Fn(&DVector) -> FeaResult<(DVector, DMatrix)>, + quadrature_order: Option, +) -> FeaResult> { + let quad_rule = element.quadrature_rule(quadrature_order)?; + let dim = element.spatial_dimension(); + let num_nodes = element.num_nodes(); + let total_dofs = num_nodes * dim; + if element_displacement.len() != total_dofs { + return Err(ElementError::MatrixComputationFailed { + reason: format!( + "element displacement has {} entries, element has {total_dofs} DOFs", + element_displacement.len() + ), + } + .into()); + } + if dim != 2 && dim != 3 { + return Err(ElementError::MatrixComputationFailed { + reason: format!("total-Lagrangian formulation needs 2-D or 3-D, got {dim}"), + } + .into()); + } + let n_strain = if dim == 2 { 3 } else { 6 }; + let pairs: Vec<(usize, usize)> = if dim == 2 { + vec![(0, 0), (1, 1), (0, 1)] + } else { + vec![(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)] + }; + + let mut internal_force = DVector::zeros(total_dofs); + for point in &quad_rule.points { + let shape_eval = element.shape_functions(&point.coords)?; + let jacobian_eval = element.jacobian(&point.coords, node_coords)?; + if !jacobian_eval.is_valid() { + return Err(ElementError::JacobianSingular { + det: jacobian_eval.determinant, + } + .into()); + } + let g = jacobian_eval.transform_derivatives(&shape_eval.derivatives)?; + let mut f: DMatrix = DMatrix::identity(dim, dim); + for a in 0..num_nodes { + for i in 0..dim { + let ui = element_displacement[a * dim + i]; + for j in 0..dim { + f[(i, j)] += ui * g[(a, j)]; + } + } + } + let c_right = f.transpose() * &f; + let mut e_voigt: DVector = DVector::zeros(n_strain); + for (k, &(i, j)) in pairs.iter().enumerate() { + let e_ij = 0.5 * (c_right[(i, j)] - if i == j { 1.0 } else { 0.0 }); + e_voigt[k] = if i == j { e_ij } else { 2.0 * e_ij }; + } + let (s_voigt, _) = constitutive(&e_voigt)?; + if s_voigt.len() != n_strain { + return Err(ElementError::MatrixComputationFailed { + reason: format!( + "constitutive closure returned stress of {} components, expected {n_strain}", + s_voigt.len() + ), + } + .into()); + } + let mut b_l: DMatrix = DMatrix::zeros(n_strain, total_dofs); + for (k, &(i, j)) in pairs.iter().enumerate() { + for a in 0..num_nodes { + for m in 0..dim { + let col = a * dim + m; + b_l[(k, col)] = if i == j { + f[(m, i)] * g[(a, i)] + } else { + f[(m, i)] * g[(a, j)] + f[(m, j)] * g[(a, i)] + }; + } + } + } + let scale = jacobian_eval.determinant().abs() * point.weight; + internal_force += b_l.transpose() * &s_voigt * scale; + } + Ok(internal_force) +} diff --git a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs index 9cc7a26..360ae64 100644 --- a/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs +++ b/crates/specialized/rtx-fea/src/solvers/sparse_ldlt/numeric.rs @@ -21,6 +21,8 @@ const PANEL: usize = 32; const TRAIL_BLOCK: usize = 64; /// Subtrees below this many multiply–adds are factorised serially. const PARALLEL_WORK: f64 = 2.0e6; +/// Trailing updates below this many multiply–adds (×2) run serially. +const PARALLEL_FRONT_WORK: f64 = 4.0e6; /// A zero, non-finite or relatively tiny pivot. #[derive(Debug, Clone)] @@ -177,8 +179,13 @@ fn partial_ldlt(front: &mut [f64], m: usize, k: usize, floor: f64) -> Result<(), } } let (left, right) = front.split_at_mut(p1 * m); - let mut c0 = 0; - while c0 < r { + let left: &[f64] = left; + let w: &[f64] = &w; + // Each column block of the trailing matrix is an independent + // gemm; large fronts run them on rayon (the arithmetic of a + // block does not depend on the schedule). + let block_update = |(t, chunk): (usize, &mut [f64])| { + let c0 = t * TRAIL_BLOCK; let c1 = (c0 + TRAIL_BLOCK).min(r); let rows = r - c0; // L_P rows (p1 + c0)..m: a view into `left` (columns p0..p1). @@ -198,14 +205,24 @@ fn partial_ldlt(front: &mut [f64], m: usize, k: usize, floor: f64) -> Result<(), nalgebra::Dyn(r), ); let mut target = DMatrixViewMut::from_slice_with_strides_generic( - &mut right[c0 * m + p1 + c0..], + &mut chunk[p1 + c0..], nalgebra::Dyn(rows), nalgebra::Dyn(c1 - c0), nalgebra::Dyn(1), nalgebra::Dyn(m), ); target.gemm(-1.0, &lp, &wb.transpose(), 1.0); - c0 = c1; + }; + if (r * r * nb) as f64 > PARALLEL_FRONT_WORK { + right + .par_chunks_mut(TRAIL_BLOCK * m) + .enumerate() + .for_each(block_update); + } else { + right + .chunks_mut(TRAIL_BLOCK * m) + .enumerate() + .for_each(block_update); } } p0 = p1; @@ -213,49 +230,144 @@ fn partial_ldlt(front: &mut [f64], m: usize, k: usize, floor: f64) -> Result<(), Ok(()) } -/// Solve `L D Lᵀ x = b` in the factor numbering, in place. +/// Subtrees below this factorisation work solve serially. +const PARALLEL_SOLVE_WORK: f64 = 2.0e7; + +/// Solve `L D Lᵀ x = b` in the factor numbering, in place, over the +/// assembly tree: the forward solve is multifrontal (each supernode +/// solves its dense unit-triangular block and hands its parent an update +/// vector, children in order), the backward solve top-down (each +/// supernode gathers its off-diagonal unknowns from its parent's front +/// vector); independent subtrees run on rayon, and every supernode's +/// arithmetic is independent of the schedule. pub(super) fn solve_in_place(sym: &Symbolic, factor: &Factor, x: &mut [f64]) { - let ns = sym.num_supernodes(); let blocks: Vec>> = factor.blocks.iter().map(|b| b.lock().unwrap()).collect(); - for s in 0..ns { - let rows = sym.front_rows(s); - let m = rows.len(); - let f = sym.start[s]; - let k = sym.start[s + 1] - f; - let block = &blocks[s]; + let blocks: Vec<&[f64]> = blocks.iter().map(|b| b.as_slice()).collect(); + // One disjoint output slice per supernode (its own columns). + let mut slots: Vec> = Vec::with_capacity(sym.num_supernodes()); + { + let mut rest: &mut [f64] = x; + for s in 0..sym.num_supernodes() { + let k = sym.start[s + 1] - sym.start[s]; + let (own, tail) = rest.split_at_mut(k); + slots.push(Mutex::new(own)); + rest = tail; + } + } + let ctx = SolveCtx { + sym, + blocks: &blocks, + slots: &slots, + }; + let roots = &sym.roots; + let run_forward = |&r: &usize| { + ctx.forward(r); + }; + if roots.len() > 1 { + roots.par_iter().for_each(run_forward); + } else { + roots.iter().for_each(run_forward); + } + let run_backward = |&r: &usize| ctx.backward(r, &[]); + if roots.len() > 1 { + roots.par_iter().for_each(run_backward); + } else { + roots.iter().for_each(run_backward); + } +} + +struct SolveCtx<'a, 'x> { + sym: &'a Symbolic, + blocks: &'a [&'a [f64]], + slots: &'a [Mutex<&'x mut [f64]>], +} + +impl SolveCtx<'_, '_> { + /// Forward (`L y = b`, then `z = D⁻¹ y` on the own columns) for the + /// subtree at `s`; returns the update vector for the parent's front. + fn forward(&self, s: usize) -> Vec { + let sym = self.sym; + let kids = sym.node_children(s); + let updates: Vec> = if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_SOLVE_WORK + { + kids.par_iter().map(|&c| self.forward(c)).collect() + } else { + kids.iter().map(|&c| self.forward(c)).collect() + }; + let m = sym.rows_ptr[s + 1] - sym.rows_ptr[s]; + let k = sym.start[s + 1] - sym.start[s]; + let block = self.blocks[s]; + let mut v = vec![0.0f64; m]; + { + let own = self.slots[s].lock().unwrap(); + v[..k].copy_from_slice(&own); + } + for (&c, update) in kids.iter().zip(&updates) { + let rel = &sym.relind[sym.relind_ptr[c]..sym.relind_ptr[c + 1]]; + for (&r, &u) in rel.iter().zip(update) { + v[r as usize] += u; + } + } + let (head, tail) = v.split_at_mut(k); for j in 0..k { - let xj = x[f + j]; + let xj = head[j]; if xj != 0.0 { let col = &block[j * m..(j + 1) * m]; - for i in (j + 1)..m { - x[rows[i]] -= col[i] * xj; + for (y, &l) in head[j + 1..].iter_mut().zip(&col[j + 1..k]) { + *y -= l * xj; + } + for (y, &l) in tail.iter_mut().zip(&col[k..]) { + *y -= l * xj; } } } + { + let mut own = self.slots[s].lock().unwrap(); + for j in 0..k { + own[j] = head[j] / block[j + j * m]; + } + } + v.drain(..k); + v } - for s in 0..ns { + + /// Backward (`Lᵀ x = z`) for the subtree at `s`, given the parent's + /// front vector (the solution at the parent's front rows). + fn backward(&self, s: usize, parent_front: &[f64]) { + let sym = self.sym; let m = sym.rows_ptr[s + 1] - sym.rows_ptr[s]; - let f = sym.start[s]; - let k = sym.start[s + 1] - f; - let block = &blocks[s]; - for j in 0..k { - x[f + j] /= block[j + j * m]; + let k = sym.start[s + 1] - sym.start[s]; + let block = self.blocks[s]; + let mut front = vec![0.0f64; m]; + let rel = &sym.relind[sym.relind_ptr[s]..sym.relind_ptr[s + 1]]; + for (slot, &r) in front[k..].iter_mut().zip(rel) { + *slot = parent_front[r as usize]; + } + { + let mut own = self.slots[s].lock().unwrap(); + let (head, tail) = front.split_at_mut(k); + head.copy_from_slice(&own); + for j in (0..k).rev() { + let col = &block[j * m..(j + 1) * m]; + let mut acc = 0.0; + for (&l, &y) in col[k..].iter().zip(tail.iter()) { + acc += l * y; + } + for i in (j + 1)..k { + acc += col[i] * head[i]; + } + head[j] -= acc; + } + own.copy_from_slice(head); } - } - for s in (0..ns).rev() { - let rows = sym.front_rows(s); - let m = rows.len(); - let f = sym.start[s]; - let k = sym.start[s + 1] - f; - let block = &blocks[s]; - for j in (0..k).rev() { - let col = &block[j * m..(j + 1) * m]; - let mut acc = x[f + j]; - for i in (j + 1)..m { - acc -= col[i] * x[rows[i]]; + let kids = sym.node_children(s); + if kids.len() > 1 && sym.subtree_work[s] > PARALLEL_SOLVE_WORK { + kids.par_iter().for_each(|&c| self.backward(c, &front)); + } else { + for &c in kids { + self.backward(c, &front); } - x[f + j] = acc; } } } diff --git a/crates/specialized/rtx-fea/tests/sparse_tangent.rs b/crates/specialized/rtx-fea/tests/sparse_tangent.rs index d23ade5..ffbe1f5 100644 --- a/crates/specialized/rtx-fea/tests/sparse_tangent.rs +++ b/crates/specialized/rtx-fea/tests/sparse_tangent.rs @@ -12,7 +12,10 @@ //! (12×2×2), 40 CSM3 steps. //! 4. `modified_newton_reuses_the_factor` — `reuse = 4`: fewer //! factorisations than solves, same march to the Newton tolerance. -//! 5. `tangent_knob_parses` — the `RTX_FEA_TANGENT` values. +//! 5. `force_only_kernel_is_bit_identical` — the TL force-only kernel +//! (modified Newton's reused iterations) = the pair's force, bit for +//! bit, on Hex20 and Quad8. +//! 6. `tangent_knob_parses` — the solver enum's defaults. //! //! Instruments (`#[ignore]`, env-driven, write under `R8G_OUT`): //! @@ -306,6 +309,34 @@ fn modified_newton_reuses_the_factor() { assert!(r < 1e-4, "modified Newton drifts: {r:.3e}"); } +#[test] +fn force_only_kernel_is_bit_identical() { + use rtx_fea::elements::total_lagrangian::{ + internal_force, internal_force_and_tangent, saint_venant_kirchhoff, + }; + let flag = Flag3d::build(Flag3dSpec::turek_hron(0.1, -0.05, 3, 1, 1)).unwrap(); + let quad = quad8_flag(3, 1); + for mesh in [&flag.mesh, &quad] { + let dim = mesh.spatial_dimension; + let constitutive = saint_venant_kirchhoff(8.0e5, 5.0e5, dim); + for (e, element) in mesh.elements.values().enumerate() { + let coords: Vec> = element + .nodes + .iter() + .map(|id| mesh.get_node(*id).unwrap().position()) + .collect(); + let fe = StandardFiniteElement::new(element.element_type, coords.clone()); + let n = coords.len() * dim; + let u = DVector::from_fn(n, |i, _| 1e-3 * ((i + 7 * e) as f64 * 0.73).sin()); + let (f_pair, _) = + internal_force_and_tangent(&fe, &coords, &u, constitutive.as_ref(), None).unwrap(); + let f_only = internal_force(&fe, &coords, &u, constitutive.as_ref(), None).unwrap(); + assert!(f_pair.norm() > 0.0); + assert_eq!(f_pair, f_only, "force-only kernel differs"); + } + } +} + #[test] fn tangent_knob_parses() { // Only the parser (the knob is read when a stepper is built). @@ -527,12 +558,14 @@ fn r8g_g2_cost() { } let breakdown = match (before, stepper.tangent_stats()) { (Some(b), Some(a)) => format!( - " | assemblies {} ({:.4} s each), factorisations {} ({:.4} s each incl. \ - scatter), solves {} ({:.4} s each), fallbacks {}, nnz(L) {}", + " | assemblies {} ({:.4} s each), factorisations {} (scatter {:.4} + \ + numeric {:.4} s each), solves {} ({:.4} s each), fallbacks {}, nnz(L) {}", a.assemblies - b.assemblies, (a.assembly_seconds - b.assembly_seconds) / (a.assemblies - b.assemblies).max(1) as f64, a.factorizations - b.factorizations, + (a.scatter_seconds - b.scatter_seconds) + / (a.factorizations - b.factorizations).max(1) as f64, (a.factor_seconds - b.factor_seconds) / (a.factorizations - b.factorizations).max(1) as f64, a.solves - b.solves,