feat(cfd-python): expose blake threshold and giersiepen models via PyO3

.gitignore

@@ -39,3 +39,8 @@ validation/references/
 # Example/optimiser output directories (all crates)
 outputs/
 report/
+.venv/
+__pycache__/
+venv/
+.venv/
+__pycache__/

crates/cfd-core/src/physics/cavitation/regimes/classifier.rs

@@ -17,13 +17,13 @@ use super::types::CavitationRegime;
 #[derive(Debug, Clone, Serialize, Deserialize)]
 pub struct CavitationRegimeClassifier<T: RealField + Copy> {
     /// Rayleigh-Plesset bubble model
-    pub(super) bubble_model: RayleighPlesset<T>,
+    pub bubble_model: RayleighPlesset<T>,
     /// Ambient pressure (Pa)
-    pub(super) ambient_pressure: T,
+    pub ambient_pressure: T,
     /// Acoustic pressure amplitude (Pa), if applicable
-    pub(super) acoustic_pressure: Option<T>,
+    pub acoustic_pressure: Option<T>,
     /// Acoustic frequency (Hz), if applicable
-    pub(super) acoustic_frequency: Option<T>,
+    pub acoustic_frequency: Option<T>,
 }
 
 impl<T: RealField + Copy + FromPrimitive> CavitationRegimeClassifier<T> {

crates/cfd-python/src/cavitation.rs

@@ -0,0 +1,66 @@
+//! Cavitation regime models wrapper for `PyO3`
+
+use cfd_core::physics::cavitation::RayleighPlesset;
+use cfd_core::physics::cavitation::CavitationRegimeClassifier;
+use pyo3::prelude::*;
+
+/// Cavitation Regime Classifier for predicting cavitation behavior
+#[pyclass(name = "CavitationRegimeClassifier")]
+pub struct PyCavitationRegimeClassifier {
+    inner: CavitationRegimeClassifier<f64>,
+}
+
+#[pymethods]
+impl PyCavitationRegimeClassifier {
+    /// Create a new Cavitation Regime Classifier
+    #[new]
+    #[pyo3(signature = (initial_radius, liquid_density, liquid_viscosity, surface_tension, vapor_pressure, polytropic_index, ambient_pressure))]
+    fn new(
+        initial_radius: f64,
+        liquid_density: f64,
+        liquid_viscosity: f64,
+        surface_tension: f64,
+        vapor_pressure: f64,
+        polytropic_index: f64,
+        ambient_pressure: f64,
+    ) -> Self {
+        let bubble_model = RayleighPlesset {
+            initial_radius,
+            liquid_density,
+            liquid_viscosity,
+            surface_tension,
+            vapor_pressure,
+            polytropic_index,
+        };
+
+        Self {
+            inner: CavitationRegimeClassifier::new(
+                bubble_model,
+                ambient_pressure,
+                None, // no acoustic pressure by default
+                None, // no acoustic frequency by default
+            ),
+        }
+    }
+
+    /// Calculate Blake threshold pressure [Pa]
+    fn blake_threshold(&self) -> f64 {
+        self.inner.blake_threshold()
+    }
+
+    /// Calculate critical Blake radius for unstable growth [m]
+    fn blake_critical_radius(&self) -> f64 {
+        self.inner.bubble_model.blake_critical_radius(self.inner.ambient_pressure)
+    }
+
+    fn __str__(&self) -> String {
+        format!(
+            "CavitationRegimeClassifier(P_ambient={:.1} Pa)",
+            self.inner.ambient_pressure
+        )
+    }
+
+    fn __repr__(&self) -> String {
+        self.__str__()
+    }
+}

crates/cfd-python/src/hemolysis.rs

@@ -0,0 +1,43 @@
+//! Hemolysis model wrappers for `PyO3`
+
+use cfd_core::physics::api::HemolysisModel;
+use pyo3::prelude::*;
+
+/// Giersiepen hemolysis power law model
+#[pyclass(name = "GiersiepenModel")]
+pub struct PyGiersiepenModel {
+    inner: HemolysisModel,
+}
+
+#[pymethods]
+impl PyGiersiepenModel {
+    /// Create standard Giersiepen model
+    #[new]
+    fn new() -> Self {
+        PyGiersiepenModel {
+            inner: HemolysisModel::giersiepen_standard(),
+        }
+    }
+
+    /// Calculate hemolysis damage index from shear stress and exposure time
+    ///
+    /// # Arguments
+    /// - `shear_stress`: Shear stress [Pa]
+    /// - `exposure_time`: Exposure time [s]
+    ///
+    /// # Returns
+    /// - Blood damage index [-]
+    fn calculate_damage(&self, shear_stress: f64, exposure_time: f64) -> PyResult<f64> {
+        self.inner
+            .damage_index(shear_stress, exposure_time)
+            .map_err(|e: cfd_core::error::Error| pyo3::exceptions::PyValueError::new_err(e.to_string()))
+    }
+
+    fn __str__(&self) -> String {
+        "GiersiepenModel()".to_string()
+    }
+
+    fn __repr__(&self) -> String {
+        self.__str__()
+    }
+}

crates/cfd-python/src/lib.rs

@@ -72,6 +72,8 @@ use pyo3::prelude::*;
 
 mod bifurcation;
 mod blood;
+mod cavitation;
+mod hemolysis;
 mod poiseuille_2d;
 mod result_types;
 mod solver_2d;
@@ -80,6 +82,8 @@ mod womersley;
 
 pub use bifurcation::{PyBifurcationSolver, PyTrifurcationResult, PyTrifurcationSolver};
 pub use blood::*;
+pub use cavitation::PyCavitationRegimeClassifier;
+pub use hemolysis::PyGiersiepenModel;
 pub use poiseuille_2d::{PyPoiseuilleConfig, PyPoiseuilleResult, PyPoiseuilleSolver};
 pub use result_types::PyBifurcationResult;
 pub use solver_2d::*;
@@ -108,6 +112,10 @@ fn cfd_python(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<PyCrossBlood>()?;
     m.add_class::<PyFahraeuasLindqvist>()?;
 
+    // Cavitation and Hemolysis
+    m.add_class::<PyCavitationRegimeClassifier>()?;
+    m.add_class::<PyGiersiepenModel>()?;
+
     // Womersley pulsatile flow
     m.add_class::<PyWomersleyNumber>()?;
     m.add_class::<PyWomersleyProfile>()?;

validation/cross_validate_rust_python.py

@@ -54,12 +54,38 @@
 print(f"  P_Blake = {P_Blake_python:.2f} Pa = {P_Blake_python/1000:.2f} kPa")
 
 if has_cfd_python:
-    # TODO: Check if cfd_python exposes Blake threshold calculation
-    # For now, document that Rust implementation is in regimes.rs
     print(f"\nRust implementation:")
     print(f"  Located in: crates/cfd-core/src/physics/cavitation/regimes.rs")
-    print(f"  Method: blake_threshold() and blake_critical_radius()")
-    print(f"  Formula matches Python implementation ✓")
+
+    if hasattr(cfd_python, 'CavitySolver2D'):
+        # Check if the required API is exposed
+        if hasattr(cfd_python, 'CavitationRegimeClassifier'):
+            classifier = cfd_python.CavitationRegimeClassifier(
+                initial_radius=R_0,
+                liquid_density=WATER_DENSITY,
+                liquid_viscosity=WATER_VISCOSITY,
+                surface_tension=WATER_SURFACE_TENSION,
+                vapor_pressure=WATER_VAPOR_PRESSURE,
+                polytropic_index=1.4,
+                ambient_pressure=P_inf
+            )
+
+            R_c_rust = classifier.blake_critical_radius()
+            P_Blake_rust = classifier.blake_threshold()
+
+            print(f"  R_c (Rust) = {R_c_rust*1e6:.4f} μm")
+            print(f"  P_Blake (Rust) = {P_Blake_rust:.2f} Pa")
+
+            # Assert equivalence
+            rc_error = abs(R_c_rust - R_c_python) / R_c_python * 100
+            pb_error = abs(P_Blake_rust - P_Blake_python) / P_Blake_python * 100
+
+            if rc_error < 0.01 and pb_error < 0.01:
+                print(f"  Cross-validation: ✓ MATCH (<0.01% diff)")
+            else:
+                print(f"  Cross-validation: ✗ MISMATCH (Rc err: {rc_error:.4f}%, Pb err: {pb_error:.4f}%)")
+        else:
+            print(f"  Skipping test: CavitationRegimeClassifier not found in cfd_python")
 else:
     print(f"\nRust verification skipped (cfd_python not available)")
 
@@ -101,12 +127,25 @@ def carreau_yasuda_python(shear_rate):
 if has_cfd_python:
     print(f"\nRust implementation:")
     print(f"  Located in: crates/cfd-core/src/physics/fluid/blood.rs")
-    print(f"  Type: CarreauYasudaBlood")
-    print(f"  Method: apparent_viscosity(shear_rate)")
 
-    # Try to test if we can create a blood model
-    # Note: This depends on cfd_python API structure
-    print(f"\n  TODO: Add cfd_python API test if blood model is exposed")
+    if hasattr(cfd_python, 'CarreauYasudaBlood'):
+        blood = cfd_python.CarreauYasudaBlood()
+        print(f"  Cross-validation:")
+
+        all_match = True
+        for gamma_dot in test_shear_rates:
+            mu_python = carreau_yasuda_python(gamma_dot)
+            mu_rust = blood.apparent_viscosity(gamma_dot)
+            error_pct = abs(mu_rust - mu_python) / mu_python * 100
+
+            if error_pct > 0.01:
+                all_match = False
+                print(f"    Mismatch at γ̇ = {gamma_dot}: Python={mu_python*1000:.4f}, Rust={mu_rust*1000:.4f} (err: {error_pct:.2f}%)")
+
+        if all_match:
+            print(f"  ✓ MATCH for all shear rates (<0.01% diff)")
+    else:
+        print(f"  Skipping test: CarreauYasudaBlood not found in cfd_python")
 else:
     print(f"\nRust verification skipped (cfd_python not available)")
 
@@ -145,9 +184,29 @@ def giersiepen_python(shear_stress, exposure_time):
 
 if has_cfd_python:
     print(f"\nRust implementation:")
-    print(f"  Located in: crates/cfd-core/src/physics/hemolysis/giersiepen.rs")
-    print(f"  Method: calculate_damage(shear_stress, exposure_time)")
-    print(f"\n  TODO: Add cfd_python API test if hemolysis model is exposed")
+    print(f"  Located in: crates/cfd-core/src/physics/hemolysis/models.rs")
+
+    if hasattr(cfd_python, 'GiersiepenModel'):
+        model = cfd_python.GiersiepenModel()
+        print(f"  Cross-validation:")
+
+        all_match = True
+        for tau, t in test_cases:
+            damage_python = giersiepen_python(tau, t)
+            damage_rust = model.calculate_damage(tau, t)
+
+            # Note: Rust model might return slightly different results based on float precision or slight variations in constants
+            # so we use a reasonable tolerance
+            error_pct = abs(damage_rust - damage_python) / damage_python * 100
+
+            if error_pct > 0.01:
+                all_match = False
+                print(f"    Mismatch at τ={tau}, t={t}: Python={damage_python:.6f}, Rust={damage_rust:.6f} (err: {error_pct:.2f}%)")
+
+        if all_match:
+            print(f"  ✓ MATCH for all test cases (<0.01% diff)")
+    else:
+        print(f"  Skipping test: GiersiepenModel not found in cfd_python")
 else:
     print(f"\nRust verification skipped (cfd_python not available)")