deep_causality_physics 0.7.0

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

//! Compressible-flow thermodynamic kernels.
//!
//! Includes ideal-gas speed of sound, specific and total enthalpy, isentropic
//! stagnation (total) pressure / temperature, and the local entropy-production
//! rate for a Newtonian fluid with heat conduction.
//!
//! Conventions: SI units throughout. `γ` is the ratio of specific heats
//! (dimensionless); `R_specific` is the specific gas constant (J/(kg·K));
//! `c_p` is the specific heat at constant pressure (J/(kg·K)).

use crate::PhysicsError;
use crate::Speed;
use crate::Temperature;
use crate::{Pressure, SpecificEnthalpy, Velocity3, VelocityGradient, ViscousStress};
use deep_causality_algebra::RealField;
use deep_causality_num::FromPrimitive;

/// Ideal-gas speed of sound `a = √(γ · R_s · T)`.
///
/// Errors when `γ · R_s · T ≤ 0` or when [`Speed::new`] rejects the result.
pub fn speed_of_sound_ideal_gas_kernel<R>(
    gamma: R,
    r_specific: R,
    temperature: &Temperature<R>,
) -> Result<Speed<R>, PhysicsError>
where
    R: RealField,
{
    let arg = gamma * r_specific * temperature.value();
    if arg <= R::zero() {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "speed_of_sound_ideal_gas_kernel: γ·R_s·T must be positive".into(),
        ));
    }
    Speed::new(arg.sqrt())
}

/// Specific enthalpy `h = c_p · T` (J/kg) for an ideal gas with constant `c_p`.
pub fn specific_enthalpy_kernel<R>(cp: R, temperature: &Temperature<R>) -> SpecificEnthalpy<R>
where
    R: RealField,
{
    SpecificEnthalpy::new_unchecked(cp * temperature.value())
}

/// Total (stagnation) enthalpy `h_0 = h + 0.5 · ‖u‖²` (J/kg).
pub fn total_enthalpy_kernel<R>(
    h: &SpecificEnthalpy<R>,
    u: &Velocity3<R>,
) -> Result<SpecificEnthalpy<R>, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let half = R::from_f64(0.5)
        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(0.5) failed".into()))?;
    let u_raw = u.value();
    let speed_sq = u_raw[0] * u_raw[0] + u_raw[1] * u_raw[1] + u_raw[2] * u_raw[2];
    Ok(SpecificEnthalpy::new_unchecked(h.value() + half * speed_sq))
}

/// Isentropic stagnation pressure
/// `p_0 = p · (1 + (γ−1)/2 · M²)^(γ/(γ−1))`.
///
/// Errors when `γ ≤ 1` (the exponent diverges) or when [`Pressure::new`]
/// rejects the result.
pub fn total_pressure_isentropic_kernel<R>(
    p: &Pressure<R>,
    mach: R,
    gamma: R,
) -> Result<Pressure<R>, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let one = R::one();
    if gamma <= one {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "total_pressure_isentropic_kernel: γ must be > 1".into(),
        ));
    }
    let half = R::from_f64(0.5)
        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(0.5) failed".into()))?;
    let base = one + (gamma - one) * half * mach * mach;
    if base <= R::zero() {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "total_pressure_isentropic_kernel: base of exponent must be positive".into(),
        ));
    }
    let exponent = gamma / (gamma - one);
    Pressure::new(p.value() * base.powf(exponent))
}

/// Isentropic stagnation temperature
/// `T_0 = T · (1 + (γ−1)/2 · M²)`.
pub fn total_temperature_isentropic_kernel<R>(
    t: &Temperature<R>,
    mach: R,
    gamma: R,
) -> Result<Temperature<R>, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let one = R::one();
    if gamma <= one {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "total_temperature_isentropic_kernel: γ must be > 1".into(),
        ));
    }
    let half = R::from_f64(0.5)
        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(0.5) failed".into()))?;
    let factor = one + (gamma - one) * half * mach * mach;
    Temperature::new(t.value() * factor)
}

/// Shared guard + stagnation factor `1 + (γ−1)/2 · M²` for the isentropic
/// ratio kernels. Inclusion-style checks (`!(x >= lo && ...)`) reject NaN.
fn isentropic_stagnation_factor<R>(mach: R, gamma: R, ctx: &str) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    if !(mach >= R::zero() && mach.is_finite()) {
        return Err(PhysicsError::PhysicalInvariantBroken(format!(
            "{ctx}: Mach number must be finite and >= 0"
        )));
    }
    if !(gamma > R::one() && gamma.is_finite()) {
        return Err(PhysicsError::PhysicalInvariantBroken(format!(
            "{ctx}: γ must be finite and > 1"
        )));
    }
    let half = R::from_f64(0.5)
        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(0.5) failed".into()))?;
    Ok(R::one() + (gamma - R::one()) * half * mach * mach)
}

/// Isentropic stagnation-to-static pressure ratio
/// `p₀/p = (1 + (γ−1)/2 · M²)^(γ/(γ−1))`.
///
/// Dimensionless; valid for a calorically perfect gas along an isentrope.
///
/// # Errors
/// [`PhysicsError::PhysicalInvariantBroken`] unless `M` is finite and `≥ 0`
/// and `γ` is finite and `> 1` (NaN inputs are rejected).
///
/// # References
/// * Anderson, "Modern Compressible Flow: With Historical Perspective,"
///   McGraw-Hill — isentropic flow relations (Ch. 3).
pub fn isentropic_pressure_ratio_kernel<R>(mach: R, gamma: R) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let factor = isentropic_stagnation_factor(mach, gamma, "isentropic_pressure_ratio_kernel")?;
    Ok(factor.powf(gamma / (gamma - R::one())))
}

/// Isentropic stagnation-to-static temperature ratio
/// `T₀/T = 1 + (γ−1)/2 · M²`.
///
/// Dimensionless; valid for a calorically perfect gas (adiabatic — no
/// isentropy needed for the temperature ratio itself).
///
/// # Errors
/// [`PhysicsError::PhysicalInvariantBroken`] unless `M` is finite and `≥ 0`
/// and `γ` is finite and `> 1` (NaN inputs are rejected).
///
/// # References
/// * Anderson, "Modern Compressible Flow: With Historical Perspective,"
///   McGraw-Hill — isentropic flow relations (Ch. 3).
pub fn isentropic_temperature_ratio_kernel<R>(mach: R, gamma: R) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    isentropic_stagnation_factor(mach, gamma, "isentropic_temperature_ratio_kernel")
}

/// Isentropic stagnation-to-static density ratio
/// `ρ₀/ρ = (1 + (γ−1)/2 · M²)^(1/(γ−1))`.
///
/// Dimensionless; valid for a calorically perfect gas along an isentrope.
///
/// # Errors
/// [`PhysicsError::PhysicalInvariantBroken`] unless `M` is finite and `≥ 0`
/// and `γ` is finite and `> 1` (NaN inputs are rejected).
///
/// # References
/// * Anderson, "Modern Compressible Flow: With Historical Perspective,"
///   McGraw-Hill — isentropic flow relations (Ch. 3).
pub fn isentropic_density_ratio_kernel<R>(mach: R, gamma: R) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let factor = isentropic_stagnation_factor(mach, gamma, "isentropic_density_ratio_kernel")?;
    Ok(factor.powf(R::one() / (gamma - R::one())))
}

/// Area–Mach relation for quasi-one-dimensional isentropic duct flow:
/// `A/A* = (1/M) · [ 2/(γ+1) · (1 + (γ−1)/2 · M²) ]^((γ+1)/(2(γ−1)))`,
/// where `A*` is the sonic-throat area.
///
/// Dimensionless. `A/A* = 1` exactly at `M = 1`; the relation is
/// double-valued in `M` for `A/A* > 1` (one subsonic, one supersonic branch).
///
/// # Errors
/// [`PhysicsError::PhysicalInvariantBroken`] unless `M` is finite and `> 0`
/// (the relation divides by `M`, and `A/A*` diverges at rest) and `γ` is
/// finite and `> 1` (NaN inputs are rejected).
///
/// # References
/// * Anderson, "Modern Compressible Flow: With Historical Perspective,"
///   McGraw-Hill — area–Mach relation, quasi-one-dimensional flow (Ch. 5).
pub fn area_mach_ratio_kernel<R>(mach: R, gamma: R) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    if !(mach > R::zero() && mach.is_finite()) {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "area_mach_ratio_kernel: Mach number must be finite and > 0".into(),
        ));
    }
    let factor = isentropic_stagnation_factor(mach, gamma, "area_mach_ratio_kernel")?;
    let two = R::from_f64(2.0)
        .ok_or_else(|| PhysicsError::NumericalInstability("R::from_f64(2.0) failed".into()))?;
    let gp1 = gamma + R::one();
    let exponent = gp1 / (two * (gamma - R::one()));
    Ok((two / gp1 * factor).powf(exponent) / mach)
}

/// Local entropy-production rate density (W/(m³·K)) for a Newtonian fluid:
/// `σ = Φ/T + κ · ‖∇T‖² / T²`
/// where `Φ = τ : ∇u` is viscous dissipation and `κ` is thermal conductivity.
///
/// For a Newtonian fluid both terms are non-negative (Clausius–Duhem
/// inequality). Errors when `T ≤ 0` (division by zero or negative absolute
/// temperature).
///
/// The `tau` argument is typed as [`ViscousStress<R>`], so the `σ ≥ 0`
/// second-law guarantee is preserved at the type level — the full Cauchy
/// stress `σ = −p I + τ` cannot be passed here without an explicit
/// conversion.
pub fn entropy_production_rate_kernel<R>(
    temperature: &Temperature<R>,
    tau: &ViscousStress<R>,
    grad_u: &VelocityGradient<R>,
    thermal_conductivity: R,
    grad_temperature: &[R; 3],
) -> Result<R, PhysicsError>
where
    R: RealField,
{
    let t = temperature.value();
    if t <= R::zero() {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "entropy_production_rate_kernel: temperature must be > 0".into(),
        ));
    }
    // Fourier's law fixes the sign of κ ≥ 0 for any physically realisable
    // material; a negative thermal conductivity would have heat flowing
    // from cold to hot and would let the kernel return σ < 0, breaking
    // the second-law (non-negativity) contract advertised above.
    if thermal_conductivity < R::zero() {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "entropy_production_rate_kernel: thermal_conductivity must be ≥ 0".into(),
        ));
    }

    // Φ = τ : ∇u
    let tv = tau.value();
    let gv = grad_u.value();
    let phi = tv[0][0] * gv[0][0]
        + tv[0][1] * gv[0][1]
        + tv[0][2] * gv[0][2]
        + tv[1][0] * gv[1][0]
        + tv[1][1] * gv[1][1]
        + tv[1][2] * gv[1][2]
        + tv[2][0] * gv[2][0]
        + tv[2][1] * gv[2][1]
        + tv[2][2] * gv[2][2];

    let grad_t_norm_sq = grad_temperature[0] * grad_temperature[0]
        + grad_temperature[1] * grad_temperature[1]
        + grad_temperature[2] * grad_temperature[2];

    Ok(phi / t + thermal_conductivity * grad_t_norm_sq / (t * t))
}