diff --git a/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs b/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs index dfa2a2b..708abf4 100644 --- a/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs +++ b/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs @@ -6,9 +6,9 @@ //! and ≈ 1 with upwind; every step is divergence-free to the solver's //! tolerance. -use rtx_cfd::mesh::patch_gen::{annulus_skewed, cartesian}; use rtx_cfd::mesh::PatchMesh; use rtx_cfd::mesh::PatchSide; +use rtx_cfd::mesh::patch_gen::{annulus_skewed, cartesian}; use rtx_cfd::solvers::incompressible::{ CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, RobinWall, @@ -433,3 +433,95 @@ async fn skewed_annulus_stokes_with_a_robin_inner_wall_keeps_second_order() -> C ); Ok(()) } + +/// The Robin wall's SIGN — the steady MMS pin above is blind to it (both +/// parts vanish at its fixed point). A still CLOSED annulus under a +/// uniform pressure `p₀` above a zero datum: the explicit part must move +/// the wall INTO the body (offset · n < 0 with n = S/|S| into the fluid, +/// the magnitude `p₀/(α (1 + g))`, `g = μ/(α d)` the viscous damping); +/// the fluid cannot follow (the outer wall is fixed), so the implicit +/// compliant term `|S| p'/α` must answer with a pressure DROP that stops +/// the recession — `p' ≈ −p₀/(1 + g)`, the net wall flux absorbed. With +/// the wrong sign the pressure would rise by the same amount. +#[tokio::test] +async fn robin_wall_recedes_under_pressure_on_both_parts() -> CfdResult<()> { + let mesh = annulus_skewed([0.0, 0.0], 0.5, 1.5, 48, 12, 0.3, 3.0)?; + let config = CfdConfig::new() + .with_density(RHO) + .with_viscosity(MU) + .with_reference_velocity(1.0) + .with_reference_length(1.0); + let params = CurvilinearParameters { + tolerance: 1e-12, + convection: PatchConvection::None, + normal_diffusion: NormalDiffusion::Explicit, + ..CurvilinearParameters::default() + }; + let mut solver = CurvilinearPisoSolver::new(config, params, mesh)?; + solver.set_boundary_velocity(|_, _, _| (0.0, 0.0)); + let n_inner = (0..solver.mesh().faces().len()) + .filter(|&f| solver.mesh().side(f) == Some(PatchSide::Inner)) + .count(); + let alpha = 4.0; + solver.set_robin_wall(Some(RobinWall { + alpha, + datum: vec![[0.0, 0.0]; n_inner], + })); + let mut field = PatchField::new(solver.mesh()); + solver.initialize(&mut field, |_, _| (0.0, 0.0)); + let p0 = 1.0; + field.p.iter_mut().for_each(|p| *p = p0); + let h = min_spacing(solver.mesh()); + let dt = 0.1 * (h * h / (4.0 * MU / RHO)).min(h); + solver.advance(&mut field, dt).await?; + let mesh = solver.mesh(); + let inner: Vec = (0..mesh.faces().len()) + .filter(|&f| mesh.side(f) == Some(PatchSide::Inner)) + .collect(); + let offsets = solver.robin_offsets(); + assert_eq!(offsets.len(), inner.len()); + let (mut worst_rel, mut net_flux, mut recession, mut p_expect) = (0.0_f64, 0.0, 0.0, 0.0); + for (j, &f) in inner.iter().enumerate() { + let face = &mesh.faces()[f]; + let s = face.s; // owner None, neigh Some: S points into the fluid + let len = (s[0] * s[0] + s[1] * s[1]).sqrt(); + let n = [s[0] / len, s[1] / len]; + let c = mesh.boundary_cell(f); + let xc = mesh.centre(c); + let d = ((face.centre[0] - xc[0]).powi(2) + (face.centre[1] - xc[1]).powi(2)).sqrt(); + let g = MU / (alpha * d); + let expect = -p0 / (alpha * (1.0 + g)); + let got = offsets[j][0] * n[0] + offsets[j][1] * n[1]; + assert!( + got < 0.0, + "face {f}: explicit offset · n = {got:.3e} (must recede)" + ); + worst_rel = worst_rel.max(((got - expect) / expect).abs()); + net_flux += field.flux[f]; + recession += got.abs() * len; + p_expect += (p0 - p0 / (1.0 + g)) / inner.len() as f64; + } + let p_mean = field.p.iter().sum::() / field.p.len() as f64; + let p_max = field.p.iter().cloned().fold(f64::MIN, f64::max); + println!( + " Robin sign pin: {} Inner faces, explicit offset −p₀/(α(1+g)) to {worst_rel:.2e}; net wall flux {net_flux:.2e} vs the recession's {recession:.2e}; pressure after the step mean {p_mean:.4} max {p_max:.4} (p₀ {p0}, compliant answer ≈ {p_expect:.4})", + inner.len() + ); + assert!( + worst_rel < 1e-9, + "explicit offset magnitude off by {worst_rel:.2e}" + ); + assert!( + net_flux.abs() < 1e-6 * recession, + "the compliant wall did not absorb the recession: net {net_flux:.3e} of {recession:.3e}" + ); + assert!( + p_max < p0, + "the pressure ROSE (max {p_max:.4} ≥ p₀): the implicit term advances the wall" + ); + assert!( + (p_mean - p_expect).abs() < 0.3 * p0, + "pressure drop {p_mean:.4} far from the compliant answer {p_expect:.4}" + ); + Ok(()) +}