carapace_core/analysis/
convergence.rs

1use nalgebra::DVector;
2
3use super::AnalysisError;
4
5/// Convergence criterion for `Algorithm::Newton`/`KrylovNewton`. Closed enum
6///. Unused by `Algorithm::Linear` (which never iterates).
7///
8/// Every variant is evaluated against the **accepted** state of an
9/// iteration (`IterationContext`): the residual reassembled *after* the
10/// iteration's correction (and any line-search rescaling) was applied, not
11/// the residual that correction was solved against. A force test can
12/// therefore never pass on a stale pre-correction residual.
13#[derive(Debug, Clone, Copy, PartialEq)]
14pub enum ConvergenceTest {
15    /// Converged once the Euclidean norm of the post-update unbalanced
16    /// force residual drops below `tol`.
17    NormUnbalance { tol: f64, max_iter: usize },
18    /// Converged once the displacement correction's norm drops below `tol`.
19    /// Does not certify force equilibrium on its own (a stiff or stalled
20    /// iteration can take a tiny step at a large residual).
21    NormDispIncr { tol: f64, max_iter: usize },
22    /// Converged once the (absolute) incremental work `|du . residual|`
23    /// (accepted correction against the accepted residual) drops below
24    /// `tol`. Pairs force with translation and moment with rotation
25    /// consistently, but its value changes with the work unit and can be
26    /// small through cancellation, so it does not certify force
27    /// equilibrium on its own either.
28    EnergyIncr { tol: f64, max_iter: usize },
29    /// Every configured criterion must pass:
30    ///
31    /// - **Force equilibrium** (always): on each reduced equation `i`,
32    ///   `|r_i| <= abs_i + relative_tol * F_ref_i`, where `abs_i` is
33    ///   `force_tol` for translational equations and `moment_tol` for
34    ///   rotational ones, and `F_ref_i` is the largest of the committed and
35    ///   predicted external force and the committed internal force on that
36    ///   equation, frozen at the start of the step. A per-equation maximum
37    ///   is used rather than a norm so many quiet DOFs can't dilute one
38    ///   unbalanced one, and the reference is never the initial or current
39    ///   residual (either can be zero) or the current load alone (which
40    ///   can pass through zero during continuation).
41    /// - **Displacement correction** (when `displacement_tol` is set): the
42    ///   largest component of the last correction is at most
43    ///   `displacement_tol`. Under arc length this is the last *undamped*
44    ///   correction, so a tiny accepted backtracking step can't fake it.
45    Combined {
46        force_tol: f64,
47        moment_tol: f64,
48        relative_tol: f64,
49        displacement_tol: Option<f64>,
50        max_iter: usize,
51    },
52}
53
54/// What a convergence test sees at the end of an iteration: the residual
55/// at the accepted (post-update) state, the correction actually applied to
56/// reach it, and — for `ConvergenceTest::Combined` — the per-equation
57/// force tolerance frozen at the start of the step.
58pub(crate) struct IterationContext<'a> {
59    pub(crate) residual: &'a DVector<f64>,
60    pub(crate) du: &'a DVector<f64>,
61    pub(crate) force_tolerance: Option<&'a ForceTolerance>,
62}
63
64impl ConvergenceTest {
65    pub(crate) fn max_iter(&self) -> usize {
66        match self {
67            ConvergenceTest::NormUnbalance { max_iter, .. }
68            | ConvergenceTest::NormDispIncr { max_iter, .. }
69            | ConvergenceTest::EnergyIncr { max_iter, .. }
70            | ConvergenceTest::Combined { max_iter, .. } => *max_iter,
71        }
72    }
73
74    /// Rejects nonfinite or nonpositive tolerances and a zero iteration
75    /// budget. `relative_tol` may be zero; the absolute tolerances may not.
76    /// Every Newton-family step runs this; it is public so wire decoders
77    /// can reject bad input up front.
78    pub fn validate(&self) -> Result<(), AnalysisError> {
79        let positive = |value: f64| value.is_finite() && value > 0.0;
80        let invalid = |field| Err(AnalysisError::InvalidOption { field });
81        match *self {
82            ConvergenceTest::NormUnbalance { tol, max_iter }
83            | ConvergenceTest::NormDispIncr { tol, max_iter }
84            | ConvergenceTest::EnergyIncr { tol, max_iter } => {
85                if !positive(tol) {
86                    return invalid("convergence.tol");
87                }
88                if max_iter == 0 {
89                    return invalid("convergence.maxIter");
90                }
91            }
92            ConvergenceTest::Combined {
93                force_tol,
94                moment_tol,
95                relative_tol,
96                displacement_tol,
97                max_iter,
98            } => {
99                if !positive(force_tol) {
100                    return invalid("convergence.forceTol");
101                }
102                if !positive(moment_tol) {
103                    return invalid("convergence.momentTol");
104                }
105                if !relative_tol.is_finite() || relative_tol < 0.0 {
106                    return invalid("convergence.relativeTol");
107                }
108                if displacement_tol.is_some_and(|tol| !positive(tol)) {
109                    return invalid("convergence.displacementTol");
110                }
111                if max_iter == 0 {
112                    return invalid("convergence.maxIter");
113                }
114            }
115        }
116        Ok(())
117    }
118
119    /// Whether this test can certify a state reached without any
120    /// correction (the step's predictor) — true for the force-based tests,
121    /// which then treat the correction as zero. Displacement/work tests
122    /// need a correction to measure, so they always take at least one.
123    pub(crate) fn accepts_predictor(&self) -> bool {
124        matches!(
125            self,
126            ConvergenceTest::NormUnbalance { .. } | ConvergenceTest::Combined { .. }
127        )
128    }
129
130    /// Whether the iteration described by `ctx` has converged.
131    pub(crate) fn check(&self, ctx: &IterationContext<'_>) -> bool {
132        if !ctx.residual.iter().all(|value| value.is_finite())
133            || !ctx.du.iter().all(|value| value.is_finite())
134        {
135            return false;
136        }
137        match self {
138            ConvergenceTest::NormUnbalance { tol, .. } => ctx.residual.norm() < *tol,
139            ConvergenceTest::NormDispIncr { tol, .. } => ctx.du.norm() < *tol,
140            ConvergenceTest::EnergyIncr { tol, .. } => ctx.du.dot(ctx.residual).abs() < *tol,
141            ConvergenceTest::Combined {
142                displacement_tol, ..
143            } => {
144                let force_ok = ctx
145                    .force_tolerance
146                    .expect("Combined convergence always receives a force tolerance")
147                    .scaled_max(ctx.residual)
148                    <= 1.0;
149                force_ok && displacement_tol.is_none_or(|tol| ctx.du.amax() <= tol)
150            }
151        }
152    }
153}
154
155/// Per-equation force tolerance for `ConvergenceTest::Combined` (and the
156/// force gate arc-length continuation applies under either force test),
157/// frozen for one step. `thresholds[i]` is the largest admissible
158/// `|r_i|`; dividing a residual by it gives a dimensionless measure where
159/// `<= 1` means equilibrium on that equation.
160#[derive(Debug, Clone)]
161pub(crate) struct ForceTolerance {
162    thresholds: DVector<f64>,
163}
164
165impl ForceTolerance {
166    /// `rotational[i]` marks moment equations. `magnitudes` are the force
167    /// vectors whose componentwise maximum forms `F_ref` (committed and
168    /// predicted external load, committed internal force).
169    pub(crate) fn combined(
170        force_tol: f64,
171        moment_tol: f64,
172        relative_tol: f64,
173        rotational: &[bool],
174        magnitudes: &[&DVector<f64>],
175    ) -> Self {
176        let thresholds = DVector::from_fn(rotational.len(), |i, _| {
177            let reference = magnitudes
178                .iter()
179                .map(|vector| vector[i].abs())
180                .fold(0.0, f64::max);
181            let absolute = if rotational[i] { moment_tol } else { force_tol };
182            absolute + relative_tol * reference
183        });
184        Self { thresholds }
185    }
186
187    /// The tolerance `ConvergenceTest::NormUnbalance` implies for arc
188    /// length's merit function: a uniform per-equation scale, so the
189    /// scaled squared norm is `||r||^2 / tol^2`.
190    pub(crate) fn uniform(n: usize, tol: f64) -> Self {
191        Self {
192            thresholds: DVector::from_element(n, tol),
193        }
194    }
195
196    /// Builds the tolerance a step needs for `test`, or `None` when `test`
197    /// doesn't use one.
198    pub(crate) fn for_test(
199        test: &ConvergenceTest,
200        rotational: &[bool],
201        magnitudes: &[&DVector<f64>],
202    ) -> Option<Self> {
203        match *test {
204            ConvergenceTest::Combined {
205                force_tol,
206                moment_tol,
207                relative_tol,
208                ..
209            } => Some(Self::combined(
210                force_tol,
211                moment_tol,
212                relative_tol,
213                rotational,
214                magnitudes,
215            )),
216            _ => None,
217        }
218    }
219
220    /// `max_i |r_i| / threshold_i`.
221    pub(crate) fn scaled_max(&self, residual: &DVector<f64>) -> f64 {
222        residual
223            .iter()
224            .zip(self.thresholds.iter())
225            .map(|(r, t)| r.abs() / t)
226            .fold(0.0, f64::max)
227    }
228
229    /// `sum_i (r_i / threshold_i)^2`.
230    pub(crate) fn scaled_norm_squared(&self, residual: &DVector<f64>) -> f64 {
231        residual
232            .iter()
233            .zip(self.thresholds.iter())
234            .map(|(r, t)| (r / t).powi(2))
235            .sum()
236    }
237}
238
239#[cfg(test)]
240mod tests {
241    use super::*;
242
243    fn combined(displacement_tol: Option<f64>) -> ConvergenceTest {
244        ConvergenceTest::Combined {
245            force_tol: 1e-3,
246            moment_tol: 1e-1,
247            relative_tol: 1e-6,
248            displacement_tol,
249            max_iter: 10,
250        }
251    }
252
253    #[test]
254    fn combined_uses_moment_tolerance_on_rotational_equations_only() {
255        let load = DVector::from_column_slice(&[1000.0, 0.0]);
256        let tol = ForceTolerance::combined(1e-3, 1e-1, 1e-6, &[false, true], &[&load]);
257        let test = combined(None);
258        let du = DVector::zeros(2);
259        // Translational threshold 1e-3 + 1e-6 * 1000 = 2e-3; rotational 1e-1.
260        let pass = DVector::from_column_slice(&[1.9e-3, 0.09]);
261        let fail_force = DVector::from_column_slice(&[2.1e-3, 0.0]);
262        let fail_moment = DVector::from_column_slice(&[0.0, 0.11]);
263        let ctx = |r| IterationContext {
264            residual: r,
265            du: &du,
266            force_tolerance: Some(&tol),
267        };
268        assert!(test.check(&ctx(&pass)));
269        assert!(!test.check(&ctx(&fail_force)));
270        assert!(!test.check(&ctx(&fail_moment)));
271    }
272
273    #[test]
274    fn small_correction_never_hides_large_residual_under_combined() {
275        let tol = ForceTolerance::combined(1e-6, 1e-6, 0.0, &[false], &[]);
276        let test = combined(Some(1.0));
277        let residual = DVector::from_column_slice(&[10.0]);
278        let du = DVector::from_column_slice(&[1e-12]);
279        assert!(!test.check(&IterationContext {
280            residual: &residual,
281            du: &du,
282            force_tolerance: Some(&tol),
283        }));
284    }
285
286    #[test]
287    fn nonfinite_state_never_converges() {
288        let test = ConvergenceTest::NormDispIncr {
289            tol: 1.0,
290            max_iter: 1,
291        };
292        let residual = DVector::from_column_slice(&[f64::NAN]);
293        let du = DVector::from_column_slice(&[0.0]);
294        assert!(!test.check(&IterationContext {
295            residual: &residual,
296            du: &du,
297            force_tolerance: None,
298        }));
299    }
300
301    #[test]
302    fn validation_rejects_bad_combined_options() {
303        assert!(combined(None).validate().is_ok());
304        let bad = ConvergenceTest::Combined {
305            force_tol: 0.0,
306            moment_tol: 1.0,
307            relative_tol: 0.0,
308            displacement_tol: None,
309            max_iter: 1,
310        };
311        assert_eq!(
312            bad.validate(),
313            Err(AnalysisError::InvalidOption {
314                field: "convergence.forceTol"
315            })
316        );
317    }
318}
carapace_core/analysis/convergence.rs

carapace_core/analysis/
convergence.rs
