deep_causality_cfd 0.2.0

Counterfactual Computational fluid dynamics solvers and the Flow DSL.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

/*
 * SPDX-License-Identifier: MIT
 * Copyright (c) 2023 - 2026. The DeepCausality Authors and Contributors. All Rights Reserved.
 */

//! The Encke↔Cowell **integrator regime switch** (Gap-3 B4) — the `grmhd/select_metric` pattern applied
//! to the trajectory integrator. The indicator `ε = a_aero / a_grav` (the g-load) is computed **from
//! state**; it is compared to **config** thresholds to choose the integrator:
//!
//! * `ε < ε_switch` — the **perturbed-conformal** (Encke-like) split: the exact KS matrix-exponential
//!   core + a small between-step aero kick, branch-cheap. The perturbative regime (coast, blackout onset).
//! * `ε ≥ ε_switch` — **direct** (Cowell) integration, accurate where aero dominates (peak dynamic
//!   pressure). Branch-cheapness is given up *observably*, by an explicit coupling decision.
//!
//! The switch carries **hysteresis** (a Schmitt trigger — separate enter/exit thresholds) so it does not
//! chatter near `ε_switch`. In the overlap band the two integrators agree (the KS Strang split is
//! 2nd-order against a direct solve), so the handover is seamless.

use deep_causality_algebra::RealField;
use deep_causality_physics::PhysicsError;

/// Which trajectory integrator the regime detector has selected.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum IntegratorRegime {
    /// Encke-like perturbed-conformal split (exact KS core + between-step aero kick) — `ε` small.
    PerturbedConformal,
    /// Cowell direct integration — `ε` large (aero-dominated).
    Direct,
}

/// The regime indicator `ε = a_aero / a_grav = |a_aero| / (GM/r²)` — the g-load, computed from state.
///
/// # Errors
/// [`PhysicsError::Singularity`] on a non-positive radius.
pub fn aero_gravity_ratio<R: RealField>(
    aero_accel: [R; 3],
    radius: R,
    gm: R,
) -> Result<R, PhysicsError> {
    if radius <= R::zero() {
        return Err(PhysicsError::Singularity(
            "aero/gravity ratio needs a positive radius".into(),
        ));
    }
    let a_aero = (aero_accel[0] * aero_accel[0]
        + aero_accel[1] * aero_accel[1]
        + aero_accel[2] * aero_accel[2])
        .sqrt();
    let a_grav = gm / (radius * radius);
    Ok(a_aero / a_grav)
}

/// A hysteresis (Schmitt-trigger) integrator switch: it flips to [`IntegratorRegime::Direct`] when the
/// g-load `ε` rises above `enter_direct`, and back to [`IntegratorRegime::PerturbedConformal`] only when
/// `ε` falls below the lower `exit_direct` — the dead band prevents chatter around `ε_switch`.
#[derive(Clone, Copy, Debug)]
pub struct RegimeSwitch<R> {
    exit_direct: R,
    enter_direct: R,
    regime: IntegratorRegime,
}

impl<R: RealField> RegimeSwitch<R> {
    /// A switch with the lower (`exit_direct`) and upper (`enter_direct`) g-load thresholds and an initial
    /// regime.
    ///
    /// **`exit_direct ≤ enter_direct` is a stated precondition, not an enforced one.** This constructor
    /// is infallible and does not reject an inverted band. Supplying `exit_direct > enter_direct` inverts
    /// the Schmitt trigger, so instead of suppressing chatter the switch flips on every step inside the
    /// band. Nothing reports it. Order the thresholds at the call site.
    pub fn new(exit_direct: R, enter_direct: R, initial: IntegratorRegime) -> Self {
        Self {
            exit_direct,
            enter_direct,
            regime: initial,
        }
    }

    /// Feed the current g-load `ε` and return the (possibly updated) regime, applying hysteresis.
    pub fn select(&mut self, epsilon: R) -> IntegratorRegime {
        self.regime = match self.regime {
            IntegratorRegime::PerturbedConformal => {
                if epsilon > self.enter_direct {
                    IntegratorRegime::Direct
                } else {
                    IntegratorRegime::PerturbedConformal
                }
            }
            IntegratorRegime::Direct => {
                if epsilon < self.exit_direct {
                    IntegratorRegime::PerturbedConformal
                } else {
                    IntegratorRegime::Direct
                }
            }
        };
        self.regime
    }

    /// The current regime.
    pub fn regime(&self) -> IntegratorRegime {
        self.regime
    }
}