runmat_analysis_fea/physics/structural/
mod.rs

1pub fn displacement_increment_norm(a: &[f64], b: &[f64]) -> f64 {
2    let mut sum = 0.0_f64;
3    for (av, bv) in a.iter().zip(b.iter()) {
4        let d = av - bv;
5        sum += d * d;
6    }
7    sum.sqrt()
8}
9
10#[derive(Debug, Clone, Copy, PartialEq)]
11pub struct TransientAdaptivityInput {
12    pub step_dt: f64,
13    pub residual_norm: f64,
14    pub residual_target: f64,
15    pub min_dt: f64,
16    pub max_dt: f64,
17    pub converged: bool,
18    pub retries: usize,
19    pub adapt_nonconverged_shrink: f64,
20    pub adapt_growth_exponent: f64,
21    pub adapt_min_scale: f64,
22    pub adapt_max_scale: f64,
23    pub adapt_retry_growth_cap: f64,
24    pub thermo_growth_limit: f64,
25    pub thermo_nonconverged_shrink: f64,
26}
27
28pub fn recommend_next_time_step(input: TransientAdaptivityInput) -> f64 {
29    if !input.converged {
30        return (input.step_dt
31            * input.adapt_nonconverged_shrink.clamp(0.2, 1.0)
32            * input.thermo_nonconverged_shrink)
33            .clamp(input.min_dt, input.max_dt);
34    }
35
36    let target = input.residual_target.max(1.0e-12);
37    let ratio = (target / input.residual_norm.max(1.0e-12)).clamp(0.25, 4.0);
38    let mut factor = ratio.powf(input.adapt_growth_exponent.clamp(0.1, 1.0));
39    factor = factor.clamp(
40        input.adapt_min_scale.clamp(0.2, 1.0),
41        input.adapt_max_scale.clamp(1.0, 2.0),
42    );
43    if input.retries > 0 {
44        factor = factor.min(input.adapt_retry_growth_cap.clamp(1.0, 1.5));
45    } else if input.residual_norm <= target * 0.1 {
46        factor = factor.max((1.0 + (input.adapt_max_scale - 1.0) * 0.6).clamp(1.0, 1.5));
47    }
48    factor = factor.min(input.thermo_growth_limit.max(0.5));
49    (input.step_dt * factor).clamp(input.min_dt, input.max_dt)
50}
51
52pub fn nonlinear_iteration_damping(
53    line_search_enabled: bool,
54    refresh_tangent: bool,
55    thermo_severity: f64,
56) -> f64 {
57    let mut damping = if line_search_enabled { 0.62 } else { 0.72 };
58    if refresh_tangent {
59        damping *= 0.85;
60    }
61    damping * (1.0 - 0.08 * thermo_severity).clamp(0.65, 1.0)
62}
63
64pub fn nonlinear_trial_residual(residual: f64, trial_scale: f64) -> f64 {
65    residual * (0.85 * trial_scale + 0.1)
66}
67
68pub fn nonlinear_spike_threshold(max_newton_iters: usize, complexity_scale: f64) -> usize {
69    ((max_newton_iters as f64) * 0.7 / complexity_scale.sqrt()).ceil() as usize
70}
71
72#[cfg(test)]
73mod tests {
74    use super::*;
75
76    #[test]
77    fn displacement_increment_norm_matches_expected_delta() {
78        let a = vec![1.0, 2.0, 3.0];
79        let b = vec![1.0, 0.0, 3.0];
80        let norm = displacement_increment_norm(&a, &b);
81        assert!((norm - 2.0).abs() < 1.0e-12);
82    }
83
84    #[test]
85    fn transient_adaptivity_recommendation_stays_bounded() {
86        let next = recommend_next_time_step(TransientAdaptivityInput {
87            step_dt: 1.0e-3,
88            residual_norm: 1.0e-7,
89            residual_target: 1.0e-6,
90            min_dt: 1.0e-6,
91            max_dt: 2.0e-2,
92            converged: true,
93            retries: 0,
94            adapt_nonconverged_shrink: 0.75,
95            adapt_growth_exponent: 0.35,
96            adapt_min_scale: 0.8,
97            adapt_max_scale: 1.25,
98            adapt_retry_growth_cap: 1.05,
99            thermo_growth_limit: 1.0,
100            thermo_nonconverged_shrink: 1.0,
101        });
102        assert!((1.0e-6..=2.0e-2).contains(&next));
103    }
104
105    #[test]
106    fn nonlinear_policy_helpers_produce_finite_values() {
107        let damping = nonlinear_iteration_damping(true, true, 0.8);
108        let trial = nonlinear_trial_residual(0.1, 0.5);
109        let spike = nonlinear_spike_threshold(24, 4.0);
110        assert!(damping.is_finite() && damping > 0.0);
111        assert!(trial.is_finite() && trial > 0.0);
112        assert!(spike >= 1);
113    }
114}
runmat_analysis_fea/physics/structural/mod.rs

runmat_analysis_fea/physics/structural/
mod.rs
