vle-thermo 0.16.0

//! Truncated virial equation of state (second coefficient only).
//!
//! The virial EOS represents the compressibility factor as a power series
//! in density (or pressure): `Z = 1 + B·ρ + C·ρ² + …`. Truncating after
//! the second term gives the simplest non-ideal corrector to the
//! ideal-gas law:
//!
//! ```text
//!   Z = 1 + B·P/(R·T)
//! ```
//!
//! where `B(T)` is the temperature-dependent **second virial coefficient**
//! (units **cm³/mol** when paired with R in kJ/(kmol·K) and P in kPa, with
//! a careful 1e3 factor — see `pitzer_b` below).
//!
//! The Pitzer correlation (Ref (5), Abbott 1989) gives `B` in
//! corresponding-states form:
//!
//! ```text
//!   B·Pc/(R·Tc) = B⁰(Tr) + ω · B¹(Tr)
//!   B⁰(Tr) = 0.083 − 0.422 / Tr^1.6
//!   B¹(Tr) = 0.139 − 0.172 / Tr^4.2
//! ```
//!
//! The truncated form is valid at **low to moderate reduced pressures**
//! (roughly Pr < Tr / 2). Beyond that, cubic EOS is the better tool —
//! `engine::eos` is the production path. Virial stays in the toolkit
//! because (a) it gives clean, closed-form fugacity / departure
//! expressions that are useful for K-value initialization, and (b) the
//! Chapter IV validation cases include a sanity check against it.
//!
//! # Source
//! - VB6 `legacy/vb6/clsVirial.cls` (pure) and `clsVirialMulticomp.cls`
//!   (mixture). The Pitzer B⁰/B¹ coefficients here are bit-identical to
//!   those used by the legacy code.

use crate::types::{Component, R_GAS};
use thiserror::Error;

/// Errors raised by the virial layer.
#[derive(Debug, Error, PartialEq)]
pub enum VirialError {
    /// Pure-component or mixture input was inconsistent (zero or
    /// negative critical properties, mismatched composition vector).
    #[error("invalid input to virial layer: {0}")]
    InvalidInput(String),
}

/// Pitzer B⁰(Tr) — non-polar component contribution.
///
/// `B⁰ = 0.083 − 0.422 / Tr^1.6`.
///
/// Returned value is **dimensionless** (it's the simple-fluid part of
/// the reduced second virial coefficient `B·Pc/(R·Tc)`).
pub fn pitzer_b0(tr: f64) -> f64 {
    0.083 - 0.422 / tr.powf(1.6)
}

/// Pitzer B¹(Tr) — acentric-factor correction.
///
/// `B¹ = 0.139 − 0.172 / Tr^4.2`.
///
/// Dimensionless. Multiplied by ω in the full Pitzer expression.
pub fn pitzer_b1(tr: f64) -> f64 {
    0.139 - 0.172 / tr.powf(4.2)
}

/// Second virial coefficient `B(T)` for a pure component, in **cm³/mol**.
///
/// `B·Pc/(R·Tc) = B⁰ + ω·B¹` →  `B = (R·Tc/Pc)·(B⁰ + ω·B¹)`.
///
/// The 1e3 factor converts the canonical-unit product
/// `R[kJ/(kmol·K)] · Tc[K] / Pc[kPa] = m³/kmol` to **cm³/mol**.
///
/// # Arguments
/// * `comp` — Component (uses `tc`, `pc`, `omega`).
/// * `t` — Temperature in **K**.
///
/// # Returns
/// `B(T)` in **cm³/mol**.
pub fn pitzer_b(comp: &Component, t: f64) -> f64 {
    let tr = t / comp.tc;
    let b_reduced = pitzer_b0(tr) + comp.omega * pitzer_b1(tr);
    // R·Tc/Pc in m³/kmol → ·1000 → cm³/mol  (1 m³/kmol = 1000 cm³/mol).
    1000.0 * R_GAS * comp.tc / comp.pc * b_reduced
}

/// dB/dT for a pure component, in **cm³/(mol·K)** — needed for H^R and S^R.
///
/// Differentiate the Pitzer correlation with respect to T (Tr = T/Tc):
///   dB⁰/dTr = 0.422 · 1.6 · Tr^(−2.6)
///   dB¹/dTr = 0.172 · 4.2 · Tr^(−5.2)
///   dB/dT  = (R·Tc/Pc)·(dB⁰/dTr + ω·dB¹/dTr) · (1/Tc)
///          = (R/Pc)·(dB⁰/dTr + ω·dB¹/dTr)
/// with the same 1e3 m³/kmol → cm³/mol unit factor.
pub fn pitzer_d_b_d_t(comp: &Component, t: f64) -> f64 {
    let tr = t / comp.tc;
    let db0_dtr = 0.422 * 1.6 * tr.powf(-2.6);
    let db1_dtr = 0.172 * 4.2 * tr.powf(-5.2);
    1000.0 * R_GAS / comp.pc * (db0_dtr + comp.omega * db1_dtr)
}

/// Compressibility factor `Z` from the truncated virial equation, pure.
///
/// `Z = 1 + B·P/(R·T)`. Returns a single value (the virial truncated at
/// the second coefficient has no liquid root — it's a vapor-only model).
///
/// # Arguments
/// * `t` — Temperature in **K**.
/// * `p` — Pressure in **kPa absolute**.
///
/// # Returns
/// Z, dimensionless.
pub fn z_factor_virial(comp: &Component, t: f64, p: f64) -> f64 {
    let b = pitzer_b(comp, t); // cm³/mol
    // P (kPa) · B (cm³/mol) / (R (kJ/(kmol·K)) · T (K))
    // Units: kPa·cm³/(mol) / (kJ/(kmol·K)·K) = kPa·cm³/mol / (kJ/kmol)
    //        = kPa·cm³/mol / (1000 J/mol·1) = kPa·cm³/(1000·J/mol).
    // Since 1 kPa·cm³ = 1 J/mol·1e-3 ... let's just multiply out:
    // R_GAS = 8.31451 kJ/(kmol·K) = 8.31451 J/(mol·K) (since kJ/kmol = J/mol).
    // So R·T in J/mol. B·P with B in cm³/mol = (1e-6 m³/mol), P in kPa = 1000 Pa:
    //   B·P → 1e-6 m³/mol · 1000 Pa = 1e-3 Pa·m³/mol = 1e-3 J/mol.
    // We want B·P/(R·T) dimensionless, R·T in J/mol → multiply B·P (cm³/mol·kPa)
    // by 1e-3 to convert to J/mol. Equivalently divide by 1000.
    1.0 + b * p / (1000.0 * R_GAS * t)
}

/// Fugacity coefficient `ln(φ)` from the truncated virial equation, pure.
///
/// Derivation (textbook): `ln(φ) = B·P/(R·T)`.
/// (Same units gymnastics as [`z_factor_virial`] — note the 1/1000.)
pub fn ln_phi_pure_virial(comp: &Component, t: f64, p: f64) -> f64 {
    let b = pitzer_b(comp, t);
    b * p / (1000.0 * R_GAS * t)
}

/// Departure enthalpy `H^R/(R·T)`, dimensionless.
///
/// `H^R/RT = (B − T·dB/dT) · P / (R·T)`.
pub fn h_departure_rt_virial(comp: &Component, t: f64, p: f64) -> f64 {
    let b = pitzer_b(comp, t);
    let db_dt = pitzer_d_b_d_t(comp, t);
    (b - t * db_dt) * p / (1000.0 * R_GAS * t)
}

/// Departure entropy `S^R/R`, dimensionless.
///
/// `S^R/R = −P · dB/dT / R`. (Multiplied by 1e−3 for the cm³ → m³ unit shift.)
pub fn s_departure_r_virial(comp: &Component, t: f64, p: f64) -> f64 {
    let db_dt = pitzer_d_b_d_t(comp, t);
    -p * db_dt / (1000.0 * R_GAS)
}

// ===========================================================================
// Multicomponent virial — second virial mixing.
// ===========================================================================

/// Build the binary B_{ij} matrix for a mixture using Pitzer combining
/// rules (Tcij = √(Tci·Tcj)·(1 − k_ij), Pcij = Zcij·R·Tcij/Vcij, …).
///
/// For M7.1 simplicity we use the simplest rules: Tcij = √(Tci·Tcj),
/// Pcij = (Pci + Pcj)/2, ωij = (ωi + ωj)/2. The Tsonopoulos/Hayden
/// O'Connell refinements (which the VB6 code optionally uses) are
/// deferred to a later milestone.
pub fn b_mix_matrix(components: &[Component], t: f64) -> Vec<Vec<f64>> {
    let n = components.len();
    let flat = b_mix_matrix_flat(components, t);
    (0..n).map(|i| flat[i * n..(i + 1) * n].to_vec()).collect()
}

/// One cross second-virial coefficient `Bᵢⱼ` in **cm³/mol**, at `t` in **K**.
///
/// Factored out so the nested-`Vec` and flat matrix builders share one copy of
/// the combining rules.
fn b_cross(components: &[Component], i: usize, j: usize, t: f64) -> f64 {
    if i == j {
        return pitzer_b(&components[i], t);
    }
    // Mixing rules for cross-term cubic properties.
    let tc_ij = (components[i].tc * components[j].tc).sqrt();
    let pc_ij = 0.5 * (components[i].pc + components[j].pc);
    let omega_ij = 0.5 * (components[i].omega + components[j].omega);
    let tr = t / tc_ij;
    let b_reduced = pitzer_b0(tr) + omega_ij * pitzer_b1(tr);
    1000.0 * R_GAS * tc_ij / pc_ij * b_reduced
}

/// [`b_mix_matrix`] in **flat row-major** storage (`mat[i*n + j]`) — audit
/// Part 2 §9.
///
/// One allocation and contiguous rows, versus `Vec<Vec<f64>>`'s one allocation
/// *per row* plus a pointer chase on every element access. Bᵢⱼ depends only on
/// temperature, so a flash at fixed T builds this once and feeds it to
/// [`ln_phi_mix_virial_flat_into`] on every iteration.
///
/// `t` in **K**; entries in **cm³/mol**.
pub fn b_mix_matrix_flat(components: &[Component], t: f64) -> Vec<f64> {
    let n = components.len();
    let mut mat = vec![0.0_f64; n * n];
    for i in 0..n {
        for j in i..n {
            let bij = b_cross(components, i, j, t);
            mat[i * n + j] = bij;
            mat[j * n + i] = bij;
        }
    }
    mat
}

/// ln φ̂ᵢ from a flat Bᵢⱼ matrix, written into a caller-owned slice — audit
/// Part 2 §9.
///
/// `ln φ̂ᵢ = P/(R·T)·(2·Σⱼ xⱼBᵢⱼ − B_mix)`. The row dot products `Σⱼ xⱼBᵢⱼ` are
/// computed **once** and then reused for `B_mix = Σᵢ xᵢ·(Bx)ᵢ`, so the matrix
/// is traversed once instead of twice (the previous route summed `B_mix` over
/// the full n² matrix and then walked it again per component).
///
/// `mat` must be `n×n` row-major from [`b_mix_matrix_flat`] at the same `t`;
/// `t` in **K**, `p` in **kPa absolute**. `row_dot` and `out` are length-n
/// scratch/output slices.
pub fn ln_phi_mix_virial_flat_into(
    mat: &[f64],
    x: &[f64],
    t: f64,
    p: f64,
    row_dot: &mut [f64],
    out: &mut [f64],
) {
    let n = x.len();
    let mut bmix = 0.0_f64;
    for i in 0..n {
        let row = &mat[i * n..(i + 1) * n];
        let dot: f64 = row.iter().zip(x).map(|(&b, &xj)| b * xj).sum();
        row_dot[i] = dot;
        bmix += x[i] * dot;
    }
    let factor = p / (1000.0 * R_GAS * t);
    for i in 0..n {
        out[i] = factor * (2.0 * row_dot[i] - bmix);
    }
}

/// Mixture second virial coefficient `B_mix(T, x)` in **cm³/mol**.
///
/// `B_mix = ΣᵢΣⱼ xᵢ xⱼ Bᵢⱼ` (Lewis-Randall quadratic mixing).
pub fn b_mix(components: &[Component], mole_fractions: &[f64], t: f64) -> Result<f64, VirialError> {
    if components.len() != mole_fractions.len() {
        return Err(VirialError::InvalidInput(format!(
            "components.len()={} but mole_fractions.len()={}",
            components.len(),
            mole_fractions.len()
        )));
    }
    let mat = b_mix_matrix(components, t);
    let n = components.len();
    let mut acc = 0.0_f64;
    for i in 0..n {
        for j in 0..n {
            acc += mole_fractions[i] * mole_fractions[j] * mat[i][j];
        }
    }
    Ok(acc)
}

/// Partial fugacity coefficient ln(φᵢ) for component i in a mixture.
///
/// `ln(φᵢ) = P/(R·T) · (2·Σⱼ xⱼ·Bᵢⱼ − B_mix)`.
///
/// Returns one ln(φᵢ) per component in the same index order as `components`.
pub fn ln_phi_mix_virial(
    components: &[Component],
    mole_fractions: &[f64],
    t: f64,
    p: f64,
) -> Result<Vec<f64>, VirialError> {
    let n = components.len();
    if n != mole_fractions.len() {
        return Err(VirialError::InvalidInput(format!(
            "components.len()={} but mole_fractions.len()={}",
            n,
            mole_fractions.len()
        )));
    }
    let mat = b_mix_matrix_flat(components, t);
    let mut row_dot = vec![0.0; n];
    let mut out = vec![0.0; n];
    ln_phi_mix_virial_flat_into(&mat, mole_fractions, t, p, &mut row_dot, &mut out);
    Ok(out)
}

/// [`ln_phi_mix_virial`] with a caller-provided Bᵢⱼ matrix.
///
/// The Bᵢⱼ matrix depends only on temperature — NOT on composition or
/// pressure — so a flash iteration updating (x, P) at fixed T should
/// build the matrix once with [`b_mix_matrix`] and call this variant
/// (M8.2 cache rule, PERFORMANCE_PROPOSAL §C2). This also fixes the
/// double matrix construction the one-shot wrapper used to do (it built
/// the matrix itself and then `b_mix` built it again).
///
/// `mat` must be the N×N matrix from [`b_mix_matrix`] at the same `t`;
/// `t` in **K**, `p` in **kPa absolute**. Returns ln(φᵢ) per component.
pub fn ln_phi_mix_virial_with_matrix(
    mat: &[Vec<f64>],
    mole_fractions: &[f64],
    t: f64,
    p: f64,
) -> Vec<f64> {
    // B_mix = ΣᵢΣⱼ xᵢxⱼBᵢⱼ, computed from the same matrix (no rebuild).
    let mut bmix = 0.0_f64;
    for (row, x_i) in mat.iter().zip(mole_fractions.iter()) {
        for (b_ij, x_j) in row.iter().zip(mole_fractions.iter()) {
            bmix += x_i * x_j * b_ij;
        }
    }
    let factor = p / (1000.0 * R_GAS * t);
    mat.iter()
        .map(|row| {
            let sum_j: f64 = row
                .iter()
                .zip(mole_fractions.iter())
                .map(|(b_ij, x_j)| b_ij * x_j)
                .sum();
            factor * (2.0 * sum_j - bmix)
        })
        .collect()
}

#[cfg(test)]
mod tests {
    use super::*;

    fn dummy_methane() -> Component {
        Component {
            name: "methane".into(),
            tc: 190.564,
            pc: 4599.0, // kPa
            omega: 0.0115,
            ..Component::default()
        }
    }

    /// The flat row-major Bᵢⱼ must be the same matrix as the nested-`Vec`
    /// build, and the single-pass fugacity kernel must agree with the
    /// two-pass one it replaced.
    #[test]
    fn flat_virial_matches_nested() {
        let ethane = Component {
            name: "ethane".into(),
            tc: 305.32,
            pc: 4872.0,
            omega: 0.099,
            ..Component::default()
        };
        let propane = Component {
            name: "propane".into(),
            tc: 369.83,
            pc: 4248.0,
            omega: 0.152,
            ..Component::default()
        };
        let comps = [dummy_methane(), ethane, propane];
        let n = comps.len();
        let (t, p) = (300.0, 200.0);
        let nested = b_mix_matrix(&comps, t);
        let flat = b_mix_matrix_flat(&comps, t);
        for i in 0..n {
            for j in 0..n {
                assert_eq!(flat[i * n + j], nested[i][j], "B[{i}][{j}]");
            }
        }
        let x = [0.5, 0.3, 0.2];
        let want = ln_phi_mix_virial_with_matrix(&nested, &x, t, p);
        let mut row_dot = vec![0.0; n];
        let mut got = vec![0.0; n];
        ln_phi_mix_virial_flat_into(&flat, &x, t, p, &mut row_dot, &mut got);
        for i in 0..n {
            assert!(
                (got[i] - want[i]).abs() <= 1e-13 * want[i].abs().max(1.0),
                "comp {i}: flat={} nested={}",
                got[i],
                want[i]
            );
        }
    }

    #[test]
    fn pitzer_b0_b1_at_critical() {
        // At Tr = 1, B⁰ = 0.083 - 0.422 = -0.339; B¹ = 0.139 - 0.172 = -0.033.
        assert!((pitzer_b0(1.0) + 0.339).abs() < 1e-6);
        assert!((pitzer_b1(1.0) + 0.033).abs() < 1e-6);
    }

    #[test]
    fn virial_z_ideal_limit() {
        // At very low pressure, Z → 1 regardless of T.
        let c = dummy_methane();
        let z = z_factor_virial(&c, 300.0, 1.0); // 1 kPa
        assert!((z - 1.0).abs() < 1e-3);
    }

    #[test]
    fn ln_phi_mix_with_matrix_matches_one_shot() {
        // M8.2: the cached-matrix variant must be bit-identical to the
        // one-shot wrapper (same matrix, same math, one construction).
        let ethane = Component {
            name: "ethane".into(),
            tc: 305.3,
            pc: 4872.0,
            omega: 0.0995,
            ..Component::default()
        };
        let comps = [dummy_methane(), ethane];
        let x = [0.6, 0.4];
        let (t, p) = (280.0, 800.0);
        let one_shot = ln_phi_mix_virial(&comps, &x, t, p).unwrap();
        let mat = b_mix_matrix(&comps, t);
        let cached = ln_phi_mix_virial_with_matrix(&mat, &x, t, p);
        for (a, b) in one_shot.iter().zip(cached.iter()) {
            assert!((a - b).abs() < 1e-15, "one-shot {a} vs cached {b}");
        }
    }
}