deep_causality_cfd 0.2.0

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

//! Compressible Newtonian Navier-Stokes regime evaluators.
//!
//! **These are pointwise right-hand-side contributions only.** Each function evaluates the RHS of one
//! conservation law at a single sample point from divergences the caller supplies. No thermodynamic
//! closure is provided or checked here: the equation of state `p(ρ, e)` that ties pressure to the
//! conserved variables is the **caller's responsibility**, and nothing in this module verifies that the
//! `p` passed in is consistent with the `ρ` and `ρE` passed in. `deep_causality_physics` ships
//! `speed_of_sound_ideal_gas_kernel` and the ideal-gas pressure kernels if an ideal-gas closure is
//! wanted; a real-gas or reacting closure is the caller's to supply.
//!
//! Three pointwise RHS kernels for the compressible NS system:
//!
//! ```text
//! ∂ρ/∂t   = − ∇·(ρ u)
//! ∂u/∂t   = − (u·∇)u − (1/ρ) ∇p + (1/ρ) ∇·τ + g
//! ∂(ρE)/∂t = − ∇·(ρ u E) − ∇·(p u) + ∇·(τ·u) − ∇·q + ρ (u·g)
//! ```
//!
//! Conserved variables: `ρ` (density), `ρu` (momentum), `ρE` (total energy
//! per unit volume), with `E = e + 0.5‖u‖²`. The momentum kernel returns the
//! Eulerian acceleration `∂u/∂t` (primitive form) for consistency with the
//! other regime evaluators; the energy kernel returns the conservative-form
//! scalar `∂(ρE)/∂t`.
//!
//! Sign convention follows continuum mechanics: viscous stress positive in
//! tension; heat-flux vector `q` follows Fourier's law `q = −κ∇T`, so the
//! `-∇·q` term in the energy equation is a heat *source*.
//!
//! Caller computes the spatial divergences (`∇·τ`, `∇·q`, `∇·(p u)`,
//! `∇·(τ·u)`, `∇·(ρ u E)`) at the sample point; these kernels do not
//! discretise space.

use deep_causality_algebra::RealField;
use deep_causality_physics::PhysicsError;
use deep_causality_physics::{AccelerationVector, Density, Velocity3, VelocityGradient};
use deep_causality_physics::{
    continuity_rhs_kernel, convective_acceleration_kernel, pressure_gradient_force_kernel,
};

/// Continuity equation RHS: `∂ρ/∂t = − u·∇ρ − ρ ∇·u`.
///
/// Reduces to `0` for incompressible divergence-free flow. Returned as a
/// scalar in kg/(m³·s).
pub fn compressible_ns_continuity_rhs<R>(
    rho: &Density<R>,
    u: &Velocity3<R>,
    grad_rho: &[R; 3],
    div_u: R,
) -> R
where
    R: RealField,
{
    continuity_rhs_kernel(rho, u, grad_rho, div_u)
}

/// Momentum equation RHS in primitive velocity form:
/// `∂u/∂t = − (u·∇)u − (1/ρ) ∇p + (1/ρ) ∇·τ + g`.
///
/// Reduces to `incompressible_ns_rhs` when the viscous-stress
/// divergence equals `ρ ν ∇²u` (constant `μ`, divergence-free flow).
///
/// - `div_tau`      — divergence of the viscous stress tensor (Pa/m)
/// - errors when `ρ = 0` (inherited from `pressure_gradient_force_kernel`)
pub fn compressible_ns_momentum_rhs<R>(
    u: &Velocity3<R>,
    grad_u: &VelocityGradient<R>,
    grad_p: &[R; 3],
    div_tau: &[R; 3],
    rho: &Density<R>,
    body_force_per_mass: &AccelerationVector<R>,
) -> Result<AccelerationVector<R>, PhysicsError>
where
    R: RealField,
{
    // ρ = 0 is rejected by pressure_gradient_force_kernel below; we propagate
    // its error rather than re-checking here. The inv_rho factor on div_tau
    // is only reached after that call succeeds, so ρ > 0 is guaranteed.
    let conv = convective_acceleration_kernel(u, grad_u).into_inner();
    let press = pressure_gradient_force_kernel(rho, grad_p)?.into_inner();
    let inv_rho = R::one() / rho.value();
    let g = body_force_per_mass.value();

    Ok(AccelerationVector::new_unchecked([
        -conv[0] + press[0] + inv_rho * div_tau[0] + g[0],
        -conv[1] + press[1] + inv_rho * div_tau[1] + g[1],
        -conv[2] + press[2] + inv_rho * div_tau[2] + g[2],
    ]))
}

/// Total-energy equation RHS in conservative form:
///
/// `∂(ρE)/∂t = − ∇·(ρ u E) − ∇·(p u) + ∇·(τ·u) − ∇·q + ρ (u·g)`
///
/// All four divergences are supplied by the caller at the sample point.
///
/// **`div_q` sign convention, which this function cannot check.** It takes `∇·q` where `q` is the
/// Fourier heat flux `q = −k∇T`, pointing down-gradient. The term enters as `−∇·q`, so `−∇·q > 0` is
/// net heat *deposited* at the point. The convention lives in this docstring rather than in the type:
/// `div_q` is a bare scalar, so passing `∇·(k∇T)` (the opposite sign) compiles and silently flips the
/// energy source. A `HeatFlux`/`HeatFluxDivergence` newtype would carry the convention in the type
/// instead; that is a follow-up, and until then the caller owns the sign.
pub fn compressible_ns_energy_rhs<R>(
    rho: &Density<R>,
    u: &Velocity3<R>,
    div_rho_u_e: R,
    div_p_u: R,
    div_tau_dot_u: R,
    div_q: R,
    body_force_per_mass: &AccelerationVector<R>,
) -> R
where
    R: RealField,
{
    let u_raw = u.value();
    let g = body_force_per_mass.value();
    let u_dot_g = u_raw[0] * g[0] + u_raw[1] * g[1] + u_raw[2] * g[2];
    -div_rho_u_e - div_p_u + div_tau_dot_u - div_q + rho.value() * u_dot_g
}