Skip to main content

deep_causality_physics/kernels/fluids/
mechanics.rs

1/*
2 * SPDX-License-Identifier: MIT
3 * Copyright (c) 2023 - 2026. The DeepCausality Authors and Contributors. All Rights Reserved.
4 */
5
6use crate::{Density, G, Length, PhysicsError, Pressure, Speed};
7use deep_causality_algebra::RealField;
8use deep_causality_num::FromPrimitive;
9
10/// Calculates hydrostatic pressure: $P = P_0 + \rho g h$.
11///
12/// Returns an error when the computed total pressure is negative — the
13/// [`Pressure<R>`] newtype rejects sub-vacuum configurations (`P < 0`).
14/// This guards against physically unrealisable inputs (e.g. a numerical
15/// experiment with a negative reference pressure `P_0` and small `\rho g h`)
16/// rather than producing a silently invalid pressure value.
17pub fn hydrostatic_pressure_kernel<R>(
18    p0: &Pressure<R>,
19    density: &Density<R>,
20    depth: &Length<R>,
21) -> Result<Pressure<R>, PhysicsError>
22where
23    R: RealField + FromPrimitive,
24{
25    let g = R::from_f64(G)
26        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(G) failed".into()))?;
27    let rho_g_h = density.value() * g * depth.value();
28    let p_total = p0.value() + rho_g_h;
29
30    Pressure::new(p_total)
31}
32
33/// Calculates pressure $P_2$ using Bernoulli's principle.
34///
35/// $P_1 + \frac{1}{2}\rho v_1^2 + \rho g h_1 = P_2 + \frac{1}{2}\rho v_2^2 + \rho g h_2$
36///
37/// Solves for $P_2$:
38/// $P_2 = P_1 + \frac{1}{2}\rho(v_1^2 - v_2^2) + \rho g(h_1 - h_2)$
39pub fn bernoulli_pressure_kernel<R>(
40    p1: &Pressure<R>,
41    v1: &Speed<R>,
42    h1: &Length<R>,
43    v2: &Speed<R>,
44    h2: &Length<R>,
45    density: &Density<R>,
46) -> Result<Pressure<R>, PhysicsError>
47where
48    R: RealField + FromPrimitive,
49{
50    let half = R::from_f64(0.5)
51        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(0.5) failed".into()))?;
52    let g = R::from_f64(G)
53        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(G) failed".into()))?;
54
55    let v1_r = v1.value();
56    let v2_r = v2.value();
57    let h1_r = h1.value();
58    let h2_r = h2.value();
59
60    let rho = density.value();
61    let term_kinetic = half * rho * (v1_r * v1_r - v2_r * v2_r);
62    let term_potential = rho * g * (h1_r - h2_r);
63
64    let p2 = p1.value() + term_kinetic + term_potential;
65
66    Pressure::new(p2)
67}