semiflow 0.15.0-beta

no_std Rust PDE solver: evolution equations via Chernoff approximation of operator semigroups (heat, diffusion, Schrodinger)
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

//! Private generic helpers for `Diffusion4thChernoff`.
//!
//! Declared as `#[path = "diffusion4_generic.rs"] mod helpers_generic;` inside
//! `diffusion4.rs` — this file is a child of that module, so `super::` works.

pub(super) use diffusion_zeta_common::{validate_a_x_generic, validate_tau_generic};
use num_traits::Float;

use super::{Diffusion4thChernoff, C1, C2, C3, W0, W1, W2};
use crate::{
    diffusion_zeta_common,
    error::SemiflowError,
    float::{from_f64, half, SemiflowFloat},
    grid_fn::GridFn1D,
};

/// γ-A inner-Strang baseline (generic, scalar path).
#[inline]
pub(super) fn gamma_a_baseline_generic<F: SemiflowFloat>(
    dc: &Diffusion4thChernoff<F>,
    tau: F,
    f: &GridFn1D<F>,
    x: F,
) -> Result<F, SemiflowError> {
    let s_half = half::<F>() * tau;
    let two = from_f64::<F>(2.0);
    let three = from_f64::<F>(3.0);
    let w0 = from_f64::<F>(W0);
    let w1 = from_f64::<F>(W1);
    let w2 = from_f64::<F>(W2);

    let x_pre = x + s_half * (dc.a_prime)(x);
    let a_at_pre = (dc.a)(x_pre);
    validate_a_x_generic(a_at_pre, x_pre)?;

    let h0 = two * Float::sqrt(a_at_pre * tau);
    let h0_3 = two * Float::sqrt(three * a_at_pre * tau);

    let center_pos = x_pre + s_half * (dc.a_prime)(x_pre);

    let near_p_raw = x_pre + h0;
    let near_p_pos = near_p_raw + s_half * (dc.a_prime)(near_p_raw);

    let near_neg_raw = x_pre - h0;
    let near_neg_pos = near_neg_raw + s_half * (dc.a_prime)(near_neg_raw);

    let far_p_raw = x_pre + h0_3;
    let far_p_pos = far_p_raw + s_half * (dc.a_prime)(far_p_raw);

    let far_neg_raw = x_pre - h0_3;
    let far_neg_pos = far_neg_raw + s_half * (dc.a_prime)(far_neg_raw);

    let center = w0 * f.sample_generic(center_pos)?;
    let near = w1 * (f.sample_generic(near_p_pos)? + f.sample_generic(near_neg_pos)?);
    let far = w2 * (f.sample_generic(far_p_pos)? + f.sample_generic(far_neg_pos)?);

    Ok(center + near + far)
}

/// Apply 7-point Fornberg FD stencil for `f^(deriv)` at `x` (generic).
#[inline]
pub(super) fn fd7_generic<F: SemiflowFloat>(
    f: &GridFn1D<F>,
    x: F,
    delta: F,
    coeffs: &[f64; 7],
    deriv: i32,
) -> Result<F, SemiflowError> {
    let ks: [f64; 7] = [-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0];
    let mut sum = F::zero();
    for j in 0..7 {
        let c = from_f64::<F>(coeffs[j]);
        let k = from_f64::<F>(ks[j]);
        sum += c * f.sample_generic(x + k * delta)?;
    }
    let denom = SemiflowFloat::libm_powi(delta, deriv);
    Ok(sum / denom)
}

/// ζ⁴ τ²-correction (generic scalar path).
#[inline]
pub(super) fn zeta4_correction_generic<F: SemiflowFloat>(
    dc: &Diffusion4thChernoff<F>,
    tau: F,
    f: &GridFn1D<F>,
    x: F,
) -> Result<F, SemiflowError> {
    let three = from_f64::<F>(3.0);
    let half_v = half::<F>();
    let quarter = half_v * half_v;
    let delta = Float::max(
        three * dc.grid.dx(),
        SemiflowFloat::libm_powf(tau, from_f64::<F>(0.75)),
    );

    let a_x = (dc.a)(x);
    let a_prime_x = (dc.a_prime)(x);
    let app_x = (dc.a_double_prime)(x);

    // Short-circuit: constant-a path.
    if a_prime_x == F::zero() && app_x == F::zero() {
        return Ok(F::zero());
    }

    let f1 = fd7_generic(f, x, delta, &C1, 1)?;
    let f2 = fd7_generic(f, x, delta, &C2, 2)?;
    let f3 = fd7_generic(f, x, delta, &C3, 3)?;

    Ok(tau
        * tau
        * (a_x * a_prime_x * f3 + (a_x * app_x * half_v) * f2 + (a_prime_x * app_x * quarter) * f1))
}

/// Apply ζ⁴ at a single grid node `i` (generic path).
#[inline]
pub(super) fn apply_at_node_generic<F: SemiflowFloat>(
    dc: &Diffusion4thChernoff<F>,
    tau: F,
    f: &GridFn1D<F>,
    i: usize,
) -> Result<F, SemiflowError> {
    let x = dc.grid.x_at(i);
    Ok(gamma_a_baseline_generic(dc, tau, f, x)? + zeta4_correction_generic(dc, tau, f, x)?)
}