From c3dbd5040adc215ba93a86696aa6f94202b69289 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Tue, 15 Sep 2026 10:46:35 -0500 Subject: [PATCH] =?UTF-8?q?P6-b=20design=202:=20the=20Robin=20wall=20on=20?= =?UTF-8?q?the=20patch's=20Inner=20side=20=E2=80=94=20RobinWall=20{=20alph?= =?UTF-8?q?a,=20datum=20}=20with=20the=20wall=20velocity=20u=5Fs=20+=20(t?= =?UTF-8?q?=5Ff=20=E2=88=92=20datum)/=CE=B1=20(explicit,=20damped=20by=20t?= =?UTF-8?q?he=20wall's=20own=20viscous=20gain=20=CE=BC/(=CE=B1=20d))=20and?= =?UTF-8?q?=20the=20compliant-wall=20pressure=20term=20|S|=20p'/=CE=B1=20i?= =?UTF-8?q?mplicit=20in=20the=20projection=20(the=20added-mass=20operator?= =?UTF-8?q?=20in=20the=20fluid's=20own=20response);=20Dirichlet=20bit=20fo?= =?UTF-8?q?r=20bit=20when=20off;=20MMS=20pin=20on=20the=20skewed=20annulus?= =?UTF-8?q?=20(orders=202.46/2.13=20at=20=CE=B1=20=3D=2010=20and=20100=20?= =?UTF-8?q?=CE=BC/h,=20the=20Dirichlet=20values;=20a=201e12=20wall=20repro?= =?UTF-8?q?duces=20Dirichlet=20to=202e-6);=20OversetPisoSolver::patch=5Fmu?= =?UTF-8?q?t;=20harness=20knob=20RTX=5FFSI2O=5FROBIN=5FALPHA=20with=20the?= =?UTF-8?q?=20previous=20pass's=20tractions=20as=20the=20datum,=20printed?= =?UTF-8?q?=20in=20the=20header?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01YJPeT6WA2e7YvAnS875AHL --- .../solvers/incompressible/curvilinear/mod.rs | 111 ++++++++++++++++- .../incompressible/curvilinear/projection.rs | 25 +++- .../solvers/incompressible/embedded/mod.rs | 2 +- .../incompressible/embedded/projection.rs | 28 ++++- .../solvers/incompressible/embedded_body.rs | 6 +- .../rtx-cfd/src/solvers/incompressible/mod.rs | 6 +- .../src/solvers/incompressible/overset/mod.rs | 4 + .../src/solvers/incompressible/piso.rs | 2 +- .../src/solvers/incompressible/polygon_sdf.rs | 6 +- .../rtx-cfd/tests/curvilinear_mms.rs | 115 +++++++++++++++++- .../rtx-fsi/tests/fsi2_harness/overset.rs | 20 ++- .../tests/fsi2_harness/overset_march.rs | 28 ++++- 12 files changed, 333 insertions(+), 20 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/mod.rs index d98507e..2e1501d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/mod.rs @@ -232,6 +232,33 @@ pub struct CurvilinearPisoSolver { acceptors: Option, time: f64, matrix: Option, + /// The Robin wall on the Inner side, if any. + robin: Option, + /// Global face index of every Inner face (Inner-face order). + robin_faces: Vec, + /// `(t_f − datum) / alpha` per Inner face for the current step. + robin_offset: Vec<[f64; 2]>, +} + +/// A Robin wall on the `Inner` side (P6-b, the coupler with the added +/// mass built in — `docs/overset_metal_campaign.md` §5.19 in omni-cortex): +/// the wall velocity is the prescribed one plus `(t_f − datum) / alpha`, +/// with `t_f` the fluid's traction on the body (the [`Self::wall_tractions`] +/// convention) and `datum` the traction the structure was loaded with — +/// a wall of impedance `alpha` (Pa·s/m) that recedes when the fluid +/// pushes harder than the structure expects. Explicit in the predictor +/// (the start-of-step traction); IMPLICIT in the pressure: the wall flux +/// answers the pressure correction with `|S| p' / alpha` (a compliant +/// wall), which is the term that carries the added-mass operator into +/// the fluid's own response per subiterate. At the coupled fixed point +/// `t_f = datum` and the wall is the Dirichlet one. +#[derive(Debug, Clone)] +pub struct RobinWall { + /// Impedance, Pa·s/m (`ρ_s h_s / Δt` for a plate of thickness `h_s`). + pub alpha: f64, + /// The structure's traction per Inner face, in Inner-face order (the + /// order [`CurvilinearPisoSolver::wall_tractions`] returns). + pub datum: Vec<[f64; 2]>, } /// The assembled pressure-correction system for one `dt` and geometry. @@ -277,9 +304,80 @@ impl CurvilinearPisoSolver { acceptors: None, time: 0.0, matrix: None, + robin: None, + robin_faces: Vec::new(), + robin_offset: Vec::new(), }) } + /// Put a [`RobinWall`] on the Inner side (`datum` per Inner face, in + /// [`Self::wall_tractions`] order) or replace its datum; the pressure + /// matrix is rebuilt. `None` restores the Dirichlet wall bit for bit. + pub fn set_robin_wall(&mut self, wall: Option) { + let faces: Vec = (0..self.mesh.faces().len()) + .filter(|&f| self.mesh.side(f) == Some(PatchSide::Inner)) + .collect(); + if let Some(w) = &wall { + assert_eq!( + w.datum.len(), + faces.len(), + "Robin datum must have one entry per Inner face" + ); + } + if self.robin_offset.len() != faces.len() { + self.robin_offset = vec![[0.0, 0.0]; faces.len()]; + } + self.robin_faces = faces; + self.robin = wall; + self.matrix = None; + } + /// The Robin wall's current velocity offsets per Inner face. + pub fn robin_offsets(&self) -> &[[f64; 2]] { + &self.robin_offset + } + /// `(t_f − datum) / alpha` on every Inner face from the field's current + /// tractions (the explicit part of the Robin wall), called at the start + /// of a step; with no Robin wall the offsets stay zero. + fn refresh_robin_offsets(&mut self, field: &PatchField, t: f64) { + let Some(w) = &self.robin else { return }; + let alpha = w.alpha; + let datum = w.datum.clone(); + let mu = self.config.viscosity; + let tractions = self.wall_tractions(field, PatchSide::Inner, t); + let faces: Vec = self.robin_faces.clone(); + for (j, (_, _, _, tf)) in tractions.iter().enumerate() { + // The wall's own viscous stress answers the wall velocity as + // μ/d; taken implicitly in the update (a plain explicit + // (t_f − datum)/α has gain (μ/d)/α and blew up at α = μ/h), the + // fixed point unchanged: offset = (t_f − datum)/α. + let f = faces[j]; + let c = self.mesh.boundary_cell(f); + let xc = self.mesh.centre(c); + let xf = self.mesh.faces()[f].centre; + let d = ((xf[0] - xc[0]).powi(2) + (xf[1] - xc[1]).powi(2)) + .sqrt() + .max(1e-300); + let g = mu / (alpha * d); + let old = self.robin_offset[j]; + self.robin_offset[j] = [ + ((tf[0] - datum[j][0]) / alpha + g * old[0]) / (1.0 + g), + ((tf[1] - datum[j][1]) / alpha + g * old[1]) / (1.0 + g), + ]; + } + } + /// The Robin offset at a point of the Inner side (the nearest face). + fn robin_offset_at(&self, x: f64, y: f64) -> (f64, f64) { + let mut best = (f64::INFINITY, [0.0, 0.0]); + for (j, &f) in self.robin_faces.iter().enumerate() { + let c = self.mesh.faces()[f].centre; + let d = (c[0] - x).powi(2) + (c[1] - y).powi(2); + if d < best.0 { + best = (d, self.robin_offset[j]); + } + } + (best.1[0], best.1[1]) + } + /// Velocity on every `Velocity` side, `(x, y, t) -> (u, v)`. pub fn set_boundary_velocity(&mut self, f: F) where @@ -475,10 +573,16 @@ impl CurvilinearPisoSolver { } pub(crate) fn boundary_velocity(&self, side: PatchSide, x: f64, y: f64, t: f64) -> (f64, f64) { - self.side_velocity[side_index(side)] + let (u, v) = self.side_velocity[side_index(side)] .as_ref() .or(self.boundary_velocity.as_ref()) - .map_or((0.0, 0.0), |f| f(x, y, t)) + .map_or((0.0, 0.0), |f| f(x, y, t)); + if side == PatchSide::Inner && self.robin.is_some() { + let (ou, ov) = self.robin_offset_at(x, y); + (u + ou, v + ov) + } else { + (u, v) + } } /// Dirichlet `p'` of acceptor cell `c`, if it is one. pub(crate) fn acceptor_correction(&self, c: usize) -> Option { @@ -568,6 +672,9 @@ impl CurvilinearPisoSolver { Some(o) => StepGeometry::new(o, &self.mesh, self.params.swept_face_rule), None => StepGeometry::stationary(&self.mesh), }; + if self.robin.is_some() { + self.refresh_robin_offsets(field, t_old); + } let old_mesh = old.as_ref().unwrap_or(&self.mesh); let mesh = &self.mesh; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/projection.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/projection.rs index 1d69c8e..156e38c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/projection.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/projection.rs @@ -6,7 +6,7 @@ use super::{CurvilinearPisoSolver, PatchField, PressureSystem, SideBc, StepGeometry}; use crate::mesh::PatchSide; use crate::solvers::incompressible::sparse_bicgstab::{ - BicgstabResult, CsrMatrix, bicgstab_jacobi, project_mean, + bicgstab_jacobi, project_mean, BicgstabResult, CsrMatrix, }; impl CurvilinearPisoSolver { @@ -116,7 +116,9 @@ impl CurvilinearPisoSolver { ] .iter() .any(|&s| self.params.boundaries.get(s) == SideBc::Outlet); - if has_outlet || self.acceptors.is_some() { + // A Robin wall absorbs the net flux through its compliance (the + // pressure system is then not pure Neumann). + if has_outlet || self.acceptors.is_some() || self.robin.is_some() { return 0.0; } let (mut net, mut total_len) = (0.0, 0.0); @@ -158,6 +160,15 @@ impl CurvilinearPisoSolver { continue; } for (f, sign) in mesh.cell_faces(c) { + if let Some(w) = &self.robin { + if mesh.side(f) == Some(PatchSide::Inner) { + // The compliant wall: outward flux `+|S| p'_c / alpha`. + let sv = mesh.faces()[f].s; + let len = (sv[0] * sv[0] + sv[1] * sv[1]).sqrt(); + tri.push((c, c, len / w.alpha)); + any_dirichlet = true; + } + } self.ops .face_gradient_coeffs(mesh, &self.params.boundaries, f, &mut coefs); if mesh.side(f).is_some() && !coefs.is_empty() { @@ -257,6 +268,16 @@ impl CurvilinearPisoSolver { for f in 0..mesh.faces().len() { field.flux[f] -= dt / rho * lp[f]; } + if let Some(w) = &self.robin { + // The compliant wall's flux answer: `δu_b = −p' S / (alpha |S|)`, + // `δF = δu_b · S = −p' |S| / alpha` in the face's own orientation. + for &f in &self.robin_faces { + let c = mesh.boundary_cell(f); + let sv = mesh.faces()[f].s; + let len = (sv[0] * sv[0] + sv[1] * sv[1]).sqrt(); + field.flux[f] -= pc[c] * len / w.alpha; + } + } for c in 0..mesh.cell_count() { if self.is_acceptor(c) { continue; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs index c54ecb4..1d9735b 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs @@ -41,7 +41,7 @@ mod projection; use super::ale::{AleBoundaries, SideBoundary}; use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind}; use super::poisson::{ - MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg, + solve_multigrid_pcg, MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, }; use super::simple::ConvectionScheme; use super::{FlowField, SolverResult}; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs index fc54806..89a8e70 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/projection.rs @@ -4,12 +4,12 @@ //! meshes before applying the correction once. use super::EmbeddedPisoSolver; -use crate::CfdResult; use crate::solvers::incompressible::ale::SideBoundary; use crate::solvers::incompressible::poisson::{ - MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg, + solve_multigrid_pcg, MultigridParameters, PoissonProblem, PoissonSolverKind, }; use crate::solvers::incompressible::{EmbeddedMask, FlowField}; +use crate::CfdResult; impl EmbeddedPisoSolver { /// The pressure-correction system of one projection as a @@ -337,28 +337,44 @@ impl EmbeddedPisoSolver { // prescribed: a domain side with velocity data, or a // non-fluid interior face. let ae = if i + 1 == nx { - if b.right == outlet { ae_outlet } else { 0.0 } + if b.right == outlet { + ae_outlet + } else { + 0.0 + } } else if self.u_is_fluid(j, i + 1) { ae_interior } else { 0.0 }; let aw = if i == 0 { - if b.left == outlet { ae_outlet } else { 0.0 } + if b.left == outlet { + ae_outlet + } else { + 0.0 + } } else if self.u_is_fluid(j, i) { ae_interior } else { 0.0 }; let an = if j + 1 == ny { - if b.top == outlet { an_outlet } else { 0.0 } + if b.top == outlet { + an_outlet + } else { + 0.0 + } } else if self.v_is_fluid(j + 1, i) { an_interior } else { 0.0 }; let as_ = if j == 0 { - if b.bottom == outlet { an_outlet } else { 0.0 } + if b.bottom == outlet { + an_outlet + } else { + 0.0 + } } else if self.v_is_fluid(j, i) { an_interior } else { diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded_body.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded_body.rs index 7e01a98..20ef8db 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded_body.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded_body.rs @@ -322,7 +322,11 @@ pub fn polygon_signed_distance(vertices: &[(f64, f64)], x: f64, y: f64) -> f64 { } } let dist = dist2.sqrt(); - if inside { -dist } else { dist } + if inside { + -dist + } else { + dist + } } /// Velocity of the point on a closed polygon nearest to `(x, y)`, where diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs index 046f3b3..c5a9fee 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs @@ -49,12 +49,12 @@ pub use boundary_conditions::{ pub use curvilinear::{ CurvilinearParameters, CurvilinearPisoSolver, CurvilinearResult, CurvilinearSolverState, NormalDiffusion, Operators, PatchBalance, PatchBoundaries, PatchConvection, PatchField, - PatchLoad, SideBc, StepGeometry, + PatchLoad, RobinWall, SideBc, StepGeometry, }; pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState}; pub use embedded_body::{ - EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_interface_velocity, - polygon_signed_distance, + polygon_interface_velocity, polygon_signed_distance, EmbeddedBody, EmbeddedMask, FaceKind, + SurfaceForce, SurfaceSample, }; pub use flow_field::FlowField; pub use overset::{ diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/overset/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/overset/mod.rs index 0d57f0d..e59280b 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/overset/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/overset/mod.rs @@ -305,6 +305,10 @@ impl OversetPisoSolver { pub fn patch(&self) -> &CurvilinearPisoSolver { &self.patch } + /// The patch solver, mutably (the Robin wall's datum per pass). + pub fn patch_mut(&mut self) -> &mut CurvilinearPisoSolver { + &mut self.patch + } /// The current overlap map. pub fn overlap(&self) -> &OverlapMap { &self.overlap diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/piso.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/piso.rs index 28f53d2..5ed4a87 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/piso.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/piso.rs @@ -38,7 +38,7 @@ //! instead of using the prescribed boundary faces that exist there. use super::poisson::{ - MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg, + solve_multigrid_pcg, MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, }; use super::{BoundaryConditions, FlowField, IncompressibleSolver, SolverResult}; use crate::{CfdConfig, CfdResult}; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/polygon_sdf.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/polygon_sdf.rs index b8679e9..4f2b907 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/polygon_sdf.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/polygon_sdf.rs @@ -206,7 +206,11 @@ impl PolygonSdf { } let dist = dist2.sqrt(); - if inside { -dist } else { dist } + if inside { + -dist + } else { + dist + } } } diff --git a/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs b/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs index 592e0e3..dfa2a2b 100644 --- a/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs +++ b/crates/specialized/rtx-cfd/tests/curvilinear_mms.rs @@ -6,10 +6,12 @@ //! and ≈ 1 with upwind; every step is divergence-free to the solver's //! tolerance. -use rtx_cfd::mesh::PatchMesh; use rtx_cfd::mesh::patch_gen::{annulus_skewed, cartesian}; +use rtx_cfd::mesh::PatchMesh; +use rtx_cfd::mesh::PatchSide; use rtx_cfd::solvers::incompressible::{ CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, + RobinWall, }; use rtx_cfd::{CfdConfig, CfdResult}; use std::f64::consts::PI; @@ -54,11 +56,47 @@ fn min_spacing(mesh: &PatchMesh) -> f64 { h } +/// The manufactured traction on the body per Inner face (the +/// `wall_tractions` convention: `S` into the fluid, `(−p S + μ (∇u + ∇uᵀ) S)/|S|`). +fn robin_datum(mesh: &PatchMesh) -> Vec<[f64; 2]> { + let mut out = Vec::new(); + for (f, face) in mesh.faces().iter().enumerate() { + if mesh.side(f) != Some(PatchSide::Inner) { + continue; + } + let [x, y] = face.centre; + let sign = if face.neigh.is_some() { 1.0 } else { -1.0 }; + let s = [sign * face.s[0], sign * face.s[1]]; + let len = (s[0] * s[0] + s[1] * s[1]).sqrt(); + let p = (PI * x).sin() * (PI * y).sin(); + let ux = PI * (PI * x).cos() * (PI * y).cos(); + let uy = -PI * (PI * x).sin() * (PI * y).sin(); + let vx = PI * (PI * x).sin() * (PI * y).sin(); + let vy = -PI * (PI * x).cos() * (PI * y).cos(); + let tx = MU * (2.0 * ux * s[0] + (uy + vx) * s[1]); + let ty = MU * ((uy + vx) * s[0] + 2.0 * vy * s[1]); + out.push([(-p * s[0] + tx) / len, (-p * s[1] + ty) / len]); + } + out +} + async fn march( mesh: PatchMesh, convection: PatchConvection, diffusion: NormalDiffusion, steady_tol: f64, +) -> CfdResult { + march_with(mesh, convection, diffusion, steady_tol, None).await +} + +/// `robin_alpha`: put a Robin wall of that impedance on the Inner side with +/// the manufactured traction as its datum (P6-b's MMS pin). +async fn march_with( + mesh: PatchMesh, + convection: PatchConvection, + diffusion: NormalDiffusion, + steady_tol: f64, + robin_alpha: Option, ) -> CfdResult { let nu = MU / RHO; let h = min_spacing(&mesh); @@ -85,6 +123,10 @@ async fn march( let mut solver = CurvilinearPisoSolver::new(config, params, mesh)?; solver.set_boundary_velocity(|x, y, _t| (u_exact(x, y), v_exact(x, y))); solver.set_momentum_source(move |x, y, _t| source(x, y, convecting)); + if let Some(alpha) = robin_alpha { + let datum = robin_datum(solver.mesh()); + solver.set_robin_wall(Some(RobinWall { alpha, datum })); + } let mut field = PatchField::new(solver.mesh()); solver.initialize(&mut field, |_, _| (0.0, 0.0)); @@ -320,3 +362,74 @@ async fn snapshot_restore_rerun_is_bit_identical() -> CfdResult<()> { assert!(max_diff == 0.0, "re-run differs by {max_diff:.3e}"); Ok(()) } + +/// P6-b's MMS pin: the annulus in the Stokes limit with a Robin wall of +/// impedance `alpha` on the Inner side (datum = the manufactured traction) +/// keeps the Dirichlet wall's order (≥ 1.8) at two impedances of the +/// viscous scale, and an effectively rigid wall (`alpha` = 1e12) +/// reproduces the Dirichlet march to rounding. +#[tokio::test] +async fn skewed_annulus_stokes_with_a_robin_inner_wall_keeps_second_order() -> CfdResult<()> { + let steady_tol = 1e-4; + let mut dirichlet = Vec::new(); + for ns in [24, 48, 96] { + let m = march( + annulus_skewed([0.0, 0.0], 0.5, 1.5, ns, ns / 4, 0.3, 3.0)?, + PatchConvection::None, + NormalDiffusion::LineImplicit, + steady_tol, + ) + .await?; + dirichlet.push(m.l2_velocity); + } + println!("dirichlet L2 {dirichlet:?} orders {:?}", orders(&dirichlet)); + for &scale in &[10.0, 100.0] { + let mut errs = Vec::new(); + for ns in [24, 48, 96] { + let mesh = annulus_skewed([0.0, 0.0], 0.5, 1.5, ns, ns / 4, 0.3, 3.0)?; + let h = min_spacing(&mesh); + let alpha = scale * MU / h; + let m = march_with( + mesh, + PatchConvection::None, + NormalDiffusion::LineImplicit, + steady_tol, + Some(alpha), + ) + .await?; + println!( + "robin alpha = {scale} μ/h ns={ns}: L2 {:.6e}, max div {:.2e}, {} steps", + m.l2_velocity, m.max_div_rel, m.steps + ); + errs.push(m.l2_velocity); + } + let o = orders(&errs); + println!("robin alpha = {scale} μ/h orders {o:?}"); + assert!( + o.iter().all(|&x| x > 1.8), + "Robin ({scale} μ/h) orders {o:?} (gate >= 1.8)" + ); + } + let mesh = annulus_skewed([0.0, 0.0], 0.5, 1.5, 48, 12, 0.3, 3.0)?; + let rigid = march_with( + mesh, + PatchConvection::None, + NormalDiffusion::LineImplicit, + steady_tol, + Some(1e12), + ) + .await?; + let rel = (rigid.l2_velocity - dirichlet[1]).abs() / dirichlet[1]; + println!( + "rigid Robin (1e12) vs Dirichlet at ns 48: L2 {:.6e} vs {:.6e} (rel {rel:.2e})", + rigid.l2_velocity, dirichlet[1] + ); + // 1e-5: the two marches stop at different steps under the 1e-4 steady + // tolerance and the Robin path skips the closed-patch flux adjustment + // (measured 1.9e-6 at ns 48). + assert!( + rel < 1e-5, + "an effectively rigid Robin wall differs from Dirichlet by {rel:.2e}" + ); + Ok(()) +} diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs index 35b8d80..56afef3 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset.rs @@ -18,7 +18,7 @@ use rtx_cfd::solvers::incompressible::{ AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField, OversetParameters, OversetPisoSolver, OversetResult, OversetSolverState, PatchConvection, - PatchField, PoissonSolverKind, SideBoundary, + PatchField, PoissonSolverKind, RobinWall, SideBoundary, }; use rtx_cfd::{CfdConfig, CfdResult}; use rtx_fea::mesh::{Mesh, NodeId}; @@ -967,6 +967,24 @@ impl OversetFluid { Ok(()) } + /// The fluid's traction on every Inner (wall) face, in the patch's + /// Inner-face order — the Robin wall's datum for the next pass. + pub fn inner_tractions(&self) -> Vec<[f64; 2]> { + self.solver + .patch() + .wall_tractions(&self.field.patch, PatchSide::Inner, self.solver.time()) + .into_iter() + .map(|(_, _, _, t)| t) + .collect() + } + /// Put the Robin wall (impedance `alpha`, datum per Inner face) on the + /// patch, or remove it. + pub fn set_robin(&mut self, alpha: f64, datum: Option>) { + self.solver + .patch_mut() + .set_robin_wall(datum.map(|d| RobinWall { alpha, datum: d })); + } + /// Drag and lift on cylinder + flag from the patch's wall stress. pub fn measure_force(&self) -> (f64, f64) { let f = self diff --git a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs index 111a8ec..7c19f38 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs @@ -100,7 +100,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov // Every setting the acceptance rule reads, printed once: P5-3 lost a // day to a floor of 2e-4 against the overnight marches' 1e-6. println!( - " coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m, fict mass α {}", + " coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows, patch convection {:?}, bg convection {:?}, patch stretch {}, fillet {} m, tip corner {} m, fict mass α {}, robin α {}", cfg.coupler, cfg.reuse, cfg.initial_relaxation, @@ -119,6 +119,7 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov super::overset::fillet(), super::overset::tip_corner(), std::env::var("RTX_FSI2O_FICT_MASS").unwrap_or_else(|_| "0".into()), + std::env::var("RTX_FSI2O_ROBIN_ALPHA").unwrap_or_else(|_| "0".into()), ); // Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the @@ -295,6 +296,22 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov } a }; + // P6-b design 2 (`docs/overset_metal_campaign.md` §5.19): the Robin wall + // on the patch — `RTX_FSI2O_ROBIN_ALPHA` (Pa·s/m; 0 = the Dirichlet wall + // bit for bit; the plate's impedance is ρ_s h_s / Δt). Each pass gives + // the fluid the tractions of the previous pass as the datum, so at the + // coupled fixed point the wall is the Dirichlet one. + let robin_alpha: f64 = std::env::var("RTX_FSI2O_ROBIN_ALPHA") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(0.0); + if robin_alpha > 0.0 { + println!( + " Robin wall: α = {robin_alpha:.3e} Pa·s/m (ρ_s h_s / Δt = {:.3e})", + case.rho_s * 0.02 / dt + ); + } + let robin_datum: RefCell> = RefCell::new(Vec::new()); // The load with the compensating term for a given previous-subiterate acceleration. let with_fict = |nodal: &[(NodeId, Vector3)], accel: &[f64]| -> Vec<(NodeId, Vector3)> { @@ -383,6 +400,9 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov let (predicted, _) = flag.borrow_mut().step(&flag_state).unwrap(); extract(&predicted) }; + if robin_alpha > 0.0 { + robin_datum.replace(fluid.borrow().inner_tractions()); + } let saved = fluid.borrow().snapshot(); type PassResult = (DynamicState, Vec<(NodeId, Vector3)>, f64, usize); let latest: RefCell> = RefCell::new(None); @@ -390,9 +410,15 @@ pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> Ov let fs = std::time::Instant::now(); let mut fl = fluid.borrow_mut(); fl.restore(&saved); + if robin_alpha > 0.0 { + fl.set_robin(robin_alpha, Some(robin_datum.borrow().clone())); + } fl.advance_subcycled(&d_n, d_candidate, cfg.subcycle, v_n.as_deref()) .expect("fluid pass"); let (nodal, conservation, faces) = fl.sample_load(d_candidate); + if robin_alpha > 0.0 { + robin_datum.replace(fl.inner_tractions()); + } t_fluid.set(t_fluid.get() + fs.elapsed().as_secs_f64()); let ss = std::time::Instant::now(); let mut flag_ref = flag.borrow_mut();