outram_foam_basic_lib/limiters.rs
1//! TVD flux limiters — field-agnostic `psi(r)` functions on plain `f64`.
2//!
3//! A **flux limiter** `psi(r)` blends a high-order (linear/central) face flux
4//! with first-order upwind to suppress spurious oscillations near sharp
5//! gradients, where `r` is the ratio of successive solution gradients. `psi = 0`
6//! is first-order upwind; `psi = 1` recovers second-order (linear) differencing.
7//!
8//! This is a **pure-`f64`, mesh-free** API so any finite-volume code (e.g. the
9//! `outram-park-fork-pflotran` solute/energy transport) can build higher-order
10//! TVD advection without depending on this crate's field/mesh types. A separate,
11//! field-tied limiter for rhoCentralFoam reconstruction lives at
12//! [`crate::fv_operators::fvc::Limiter`]; this module is the reusable,
13//! general one, and the two should be consolidated eventually.
14//!
15//! # Provenance (translated from OpenFOAM upstream source)
16//!
17//! Each limiter here is a Rust translation of the corresponding `limiter()`
18//! method in OpenFOAM's
19//! `src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/<name>/<name>.H`,
20//! **Copyright (C) 2011-2022 OpenFOAM Foundation**, GNU General Public License
21//! version 3 or later (this crate is GPL-3.0). Source read from
22//! `github.com/OpenFOAM/OpenFOAM-dev` (master) on 2026-07-22. OpenFOAM® is a
23//! registered trademark of OpenCFD Ltd (ESI Group); this is an independent
24//! translation, not an official OpenFOAM product (see the workspace
25//! `TRADEMARKS.md`).
26//!
27//! The exact upstream expression is quoted in each variant's doc comment. Only
28//! OpenFOAM's **r-based** limiters are ported: OpenFOAM's NVD-based schemes
29//! (`QUICK`, `Gamma`, `SFCD`, `Phi`) use actual cell values rather than a pure
30//! `psi(r)` and are **not** representable here, so they are deliberately omitted
31//! rather than approximated.
32
33/// A TVD flux limiter `psi(r)`, translated from OpenFOAM's r-based
34/// `limitedSchemes`. `psi = 0` is first-order upwind, `psi = 1` is second-order
35/// linear; the TVD variants clip extrema (`psi(r) = 0` for `r <= 0`).
36#[derive(Debug, Clone, Copy, PartialEq)]
37pub enum FluxLimiter {
38 /// First-order `upwind`: `psi = 0`.
39 Upwind,
40 /// Unlimited `linear` (central) differencing: `psi = 1` (2nd order, not TVD).
41 Linear,
42 /// `vanLeer`: upstream `(r + mag(r))/(1 + mag(r))`. Smooth, symmetric.
43 VanLeer,
44 /// `vanAlbada`: upstream `r*(r + 1)/(sqr(r) + 1)` with `r = max(0, r)`. Symmetric.
45 VanAlbada,
46 /// `Minmod`: upstream `max(min(r, 1), 0)`. Most diffusive TVD limiter. Symmetric.
47 Minmod,
48 /// `SuperBee`: upstream `max(max(min(2r, 1), min(r, 2)), 0)`. Most compressive. Symmetric.
49 SuperBee,
50 /// `MUSCL`: upstream `max(min(min(2r, 0.5r + 0.5), 2), 0)`. Symmetric.
51 Muscl,
52 /// `UMIST`: upstream `max(min(min(min(2r, 0.75r + 0.25), 0.25r + 0.75), 2), 0)`.
53 /// Third-order biased (not symmetric).
54 Umist,
55 /// `OSPRE`: upstream `1.5 r (r + 1)/(r (r + 1) + 1)` with `r = max(0, r)`. Symmetric.
56 Ospre,
57 /// `limitedLinear(k)`: upstream `max(min((2/k) r, 1), 0)`, coefficient
58 /// `k` in `[0, 1]` (`k -> 0` approaches unlimited linear, `k = 1` most
59 /// limited). A k-blended bounded scheme — not strictly within the classic
60 /// `psi <= 2r` Sweby envelope for small `k`.
61 LimitedLinear(f64),
62}
63
64impl FluxLimiter {
65 /// The flux-limiter function `psi(r)`. `r` may be any `f64`; a non-finite `r`
66 /// returns `0`. For the TVD limiters `r <= 0` (a local extremum) returns `0`.
67 /// The result is clamped non-negative.
68 pub fn psi(&self, r: f64) -> f64 {
69 if !r.is_finite() {
70 return 0.0;
71 }
72 let val = match *self {
73 FluxLimiter::Upwind => 0.0,
74 FluxLimiter::Linear => 1.0,
75 // vanLeer: (r + |r|)/(1 + |r|)
76 FluxLimiter::VanLeer => (r + r.abs()) / (1.0 + r.abs()),
77 // vanAlbada: r*(r+1)/(r^2+1), r = max(0,r)
78 FluxLimiter::VanAlbada => {
79 let r = r.max(0.0);
80 r * (r + 1.0) / (r * r + 1.0)
81 }
82 // Minmod: max(min(r,1),0)
83 FluxLimiter::Minmod => r.min(1.0).max(0.0),
84 // SuperBee: max(max(min(2r,1), min(r,2)), 0)
85 FluxLimiter::SuperBee => (2.0 * r).min(1.0).max(r.min(2.0)).max(0.0),
86 // MUSCL: max(min(min(2r, 0.5r+0.5), 2), 0)
87 FluxLimiter::Muscl => (2.0 * r).min(0.5 * r + 0.5).min(2.0).max(0.0),
88 // UMIST: max(min(min(min(2r, 0.75r+0.25), 0.25r+0.75), 2), 0)
89 FluxLimiter::Umist => (2.0 * r)
90 .min(0.75 * r + 0.25)
91 .min(0.25 * r + 0.75)
92 .min(2.0)
93 .max(0.0),
94 // OSPRE: 1.5 r (r+1)/(r(r+1)+1), r = max(0,r)
95 FluxLimiter::Ospre => {
96 let r = r.max(0.0);
97 1.5 * r * (r + 1.0) / (r * (r + 1.0) + 1.0)
98 }
99 // limitedLinear(k): max(min((2/k) r, 1), 0), twoByk = 2/max(k, small)
100 FluxLimiter::LimitedLinear(k) => {
101 let two_by_k = 2.0 / k.max(1e-30);
102 (two_by_k * r).min(1.0).max(0.0)
103 }
104 };
105 val.max(0.0)
106 }
107
108 /// True if this is a second-order **TVD** limiter (everything except
109 /// [`Upwind`](Self::Upwind) and [`Linear`](Self::Linear)).
110 pub fn is_tvd(&self) -> bool {
111 !matches!(self, FluxLimiter::Upwind | FluxLimiter::Linear)
112 }
113
114 /// True if the limiter is **symmetric** (`psi(r)/r == psi(1/r)`), the
115 /// property of the second-order symmetric TVD limiters. UMIST is
116 /// third-order-biased (asymmetric); `limitedLinear` is a k-blend and not in
117 /// this set.
118 pub fn is_symmetric(&self) -> bool {
119 matches!(
120 self,
121 FluxLimiter::VanLeer
122 | FluxLimiter::VanAlbada
123 | FluxLimiter::Minmod
124 | FluxLimiter::SuperBee
125 | FluxLimiter::Muscl
126 | FluxLimiter::Ospre
127 )
128 }
129}
130
131#[cfg(test)]
132mod tests {
133 use super::*;
134
135 /// The classic r-based TVD limiters that must sit inside the Sweby envelope.
136 const CLASSIC: &[FluxLimiter] = &[
137 FluxLimiter::VanLeer,
138 FluxLimiter::VanAlbada,
139 FluxLimiter::Minmod,
140 FluxLimiter::SuperBee,
141 FluxLimiter::Muscl,
142 FluxLimiter::Umist,
143 FluxLimiter::Ospre,
144 ];
145
146 #[test]
147 fn consistency_at_r_one() {
148 for l in CLASSIC {
149 assert!((l.psi(1.0) - 1.0).abs() < 1e-12, "{l:?} psi(1) != 1");
150 }
151 assert_eq!(FluxLimiter::Linear.psi(1.0), 1.0);
152 assert_eq!(FluxLimiter::Upwind.psi(1.0), 0.0);
153 }
154
155 #[test]
156 fn extrema_are_clipped() {
157 for l in CLASSIC {
158 for &r in &[-3.0, -1.0, -0.1, 0.0] {
159 assert_eq!(l.psi(r), 0.0, "{l:?} psi({r}) != 0");
160 }
161 }
162 }
163
164 #[test]
165 fn within_sweby_tvd_envelope() {
166 // Sufficient TVD constraints: 0 <= psi <= 2r and psi <= 2 for r > 0.
167 for l in CLASSIC {
168 let mut r = 0.01;
169 while r <= 6.0 {
170 let p = l.psi(r);
171 assert!(p >= 0.0, "{l:?} psi({r}) < 0");
172 assert!(p <= 2.0 * r + 1e-12, "{l:?} psi({r})={p} > 2r");
173 assert!(p <= 2.0 + 1e-12, "{l:?} psi({r})={p} > 2");
174 r += 0.01;
175 }
176 }
177 }
178
179 #[test]
180 fn symmetric_limiters_satisfy_symmetry() {
181 for l in CLASSIC.iter().filter(|l| l.is_symmetric()) {
182 for &r in &[0.2, 0.5, 0.75, 1.3, 2.0, 4.0] {
183 let lhs = l.psi(r) / r;
184 let rhs = l.psi(1.0 / r);
185 assert!(
186 (lhs - rhs).abs() < 1e-10,
187 "{l:?} asymmetric at r={r}: {lhs} vs {rhs}"
188 );
189 }
190 }
191 }
192
193 #[test]
194 fn known_spot_values() {
195 assert!((FluxLimiter::SuperBee.psi(0.5) - 1.0).abs() < 1e-12);
196 assert!((FluxLimiter::SuperBee.psi(2.0) - 2.0).abs() < 1e-12);
197 assert!((FluxLimiter::Minmod.psi(0.5) - 0.5).abs() < 1e-12);
198 assert!((FluxLimiter::Minmod.psi(3.0) - 1.0).abs() < 1e-12);
199 assert!((FluxLimiter::VanLeer.psi(1.0) - 1.0).abs() < 1e-12);
200 }
201
202 #[test]
203 fn limited_linear_matches_upstream_form() {
204 // psi = max(min((2/k) r, 1), 0); bounded in [0,1]; psi(1)=1 for any k<=1.
205 for &k in &[0.25, 0.5, 1.0] {
206 let l = FluxLimiter::LimitedLinear(k);
207 for &r in &[-1.0, 0.0, 0.1, 0.5, 1.0, 3.0] {
208 let p = l.psi(r);
209 assert!(
210 (0.0..=1.0).contains(&p),
211 "limitedLinear({k}) psi({r})={p} out of [0,1]"
212 );
213 }
214 assert!((l.psi(1.0) - 1.0).abs() < 1e-12);
215 assert_eq!(l.psi(-1.0), 0.0);
216 }
217 }
218
219 #[test]
220 fn non_finite_r_is_safe() {
221 for l in CLASSIC {
222 assert_eq!(l.psi(f64::NAN), 0.0);
223 assert_eq!(l.psi(f64::INFINITY), 0.0);
224 }
225 }
226}