diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs index 52f66dd..ca420ff 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -16,6 +16,7 @@ use super::Grid; use super::body::Body; use super::cut::CutGeometry; +use super::exchange::in_load_window; use super::field::Field; use super::step::{Boundaries, Side}; use super::wall::{FaceKind, Mask}; @@ -589,37 +590,6 @@ impl Mask { flux } - /// The cut-cell load route: the force on the body from the operators - /// themselves — `Σ_c p_c W_c` over the cells plus the implicit wall - /// shear `Σ_f μ A_w (u_f − U_b)/d_f` over the unknown faces. `None` - /// without a cut geometry. - pub fn cut_wall_force(&self, body: &Body, f: &Field, mu: f64, t: f64) -> Option<[f64; 3]> { - let (p, s) = self.cut_wall_force_parts(body, f, mu, t)?; - let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, None)?; - Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]) - } - - /// The cut-cell load route restricted to the cells (and faces) of the - /// planes `k0..k1`, divided by the slab's thickness: the load per unit - /// span on a body's mid-section. - pub fn cut_wall_force_per_span( - &self, - body: &Body, - f: &Field, - mu: f64, - t: f64, - (k0, k1): (usize, usize), - ) -> Option<[f64; 3]> { - let (p, s) = self.cut_wall_force_parts_in(body, f, mu, t, Some((k0, k1)))?; - let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, Some((k0, k1)))?; - let lz = (k1 - k0) as f64 * self.grid.dz; - Some([ - (p[0] + s[0] + x[0]) / lz, - (p[1] + s[1] + x[1]) / lz, - (p[2] + s[2] + x[2]) / lz, - ]) - } - /// The cut-cell load route split into its pressure and shear parts. pub fn cut_wall_force_parts( &self, @@ -631,7 +601,7 @@ impl Mask { self.cut_wall_force_parts_in(body, f, mu, t, None) } - fn cut_wall_force_parts_in( + pub(super) fn cut_wall_force_parts_in( &self, body: &Body, f: &Field, @@ -646,8 +616,9 @@ impl Mask { let mut pressure = [0.0; 3]; let mut force = [0.0; 3]; for (idx, w) in cut.wall.iter().enumerate() { - let k = g.kji(idx).0; - if self.cell_fluid[idx] && k >= k0 && k < k1 { + let (k, _, i) = g.kji(idx); + if self.cell_fluid[idx] && k >= k0 && k < k1 && in_load_window((i as f64 + 0.5) * g.dx) + { for c in 0..3 { pressure[c] += f.p[idx] * w[c]; } @@ -676,7 +647,7 @@ impl Mask { 1 => self.v_kind[idx], _ => self.w_kind[idx], }; - if kind != FaceKind::Fluid { + if kind != FaceKind::Fluid || !in_load_window(lat.face_position(c, p)[0]) { continue; } let cv = self.cv_geometry(c, p); diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs index c0c2ee3..8a06e1d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs @@ -13,7 +13,55 @@ use super::field::Field; use super::wall::{FaceKind, Mask}; use crate::solvers::incompressible::ConvectionScheme; +/// DIAGNOSTIC: an x window on every load route (`set_load_window`): cells and +/// faces outside `[x0, x1)` are skipped — the cylinder and the flag read +/// apart. Process-wide; `None` (the default) reads the whole body. +static LOAD_WINDOW: std::sync::Mutex> = std::sync::Mutex::new(None); + +/// Set or clear the diagnostic x window of the load routes. +pub fn set_load_window(window: Option<(f64, f64)>) { + *LOAD_WINDOW.lock().expect("load window") = window; +} + +pub(super) fn in_load_window(x: f64) -> bool { + LOAD_WINDOW + .lock() + .expect("load window") + .is_none_or(|(x0, x1)| x >= x0 && x < x1) +} + impl Mask { + /// The cut-cell load route: the force on the body from the operators + /// themselves — `Σ_c p_c W_c` over the cells plus the implicit wall + /// shear `Σ_f μ A_w (u_f − U_b)/d_f` over the unknown faces. `None` + /// without a cut geometry. + pub fn cut_wall_force(&self, body: &Body, f: &Field, mu: f64, t: f64) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_parts(body, f, mu, t)?; + let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, None)?; + Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]) + } + + /// The cut-cell load route restricted to the cells (and faces) of the + /// planes `k0..k1`, divided by the slab's thickness: the load per unit + /// span on a body's mid-section. + pub fn cut_wall_force_per_span( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + (k0, k1): (usize, usize), + ) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_parts_in(body, f, mu, t, Some((k0, k1)))?; + let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, Some((k0, k1)))?; + let lz = (k1 - k0) as f64 * self.grid.dz; + Some([ + (p[0] + s[0] + x[0]) / lz, + (p[1] + s[1] + x[1]) / lz, + (p[2] + s[2] + x[2]) / lz, + ]) + } + /// The momentum the fluid's face control volumes exchange with the /// prescribed faces beside them, as a force on the body (the negative /// of the force on the fluid), over the z planes `planes` (all when @@ -27,6 +75,23 @@ impl Mask { t: f64, planes: Option<(usize, usize)>, ) -> Option<[f64; 3]> { + let (d, c) = self.cut_wall_exchange_parts(body, f, mu, rho, t, planes)?; + Some([d[0] + c[0], d[1] + c[1], d[2] + c[2]]) + } + + /// The exchange split into its DIFFUSIVE and CONVECTIVE parts (forces on + /// the body). On a wall at rest the convective part is the scheme's + /// flux correction only; on a moving wall it carries `ρ m (u_face − u_f)` + /// with `m` the wall's own swept flux — O(v_wall h / ν) times the shear. + pub fn cut_wall_exchange_parts( + &self, + body: &Body, + f: &Field, + mu: f64, + rho: f64, + t: f64, + planes: Option<(usize, usize)>, + ) -> Option<([f64; 3], [f64; 3])> { let _ = body; let _ = t; self.cut.as_ref()?; @@ -55,6 +120,7 @@ impl Mask { |a: [i64; 3], b: [i64; 3], s: i64| [a[0] + s * b[0], a[1] + s * b[1], a[2] + s * b[2]]; let upwind = |m: f64, up: f64, dn: f64| if m >= 0.0 { up } else { dn }; let mut force = [0.0; 3]; + let mut convective = [0.0; 3]; for c in 0..3 { let (ir, jr, kr) = match c { 0 => (1..nx, 0..ny, k0..k1), @@ -67,7 +133,9 @@ impl Mask { for i in ir.clone() { let p = [i as i64, j as i64, k as i64]; let idx = lat.face(c, p).expect("face"); - if kind(c, idx) != FaceKind::Fluid { + if kind(c, idx) != FaceKind::Fluid + || !in_load_window(lat.face_position(c, p)[0]) + { continue; } let cv = self.cv_geometry(c, p); @@ -123,9 +191,9 @@ impl Mask { scheme.face_correction(up2, un, u0) }; let u_face = upwind(m_plus, u0, un) + delta; - let on_fluid = -rho * m_plus * (u_face - u0) - + mu * cv.ap[d][1] * a_d * (un - u0) / solid_spacing(1.0); - force[c] -= on_fluid; + convective[c] -= -rho * m_plus * (u_face - u0); + force[c] -= + mu * cv.ap[d][1] * a_d * (un - u0) / solid_spacing(1.0); } } // Minus side. @@ -140,9 +208,9 @@ impl Mask { scheme.face_correction(up1, u0, ud) }; let u_face = upwind(m_minus, ud, u0) + delta; - let on_fluid = rho * m_minus * (u_face - u0) - + mu * cv.ap[d][0] * a_d * (ud - u0) / solid_spacing(-1.0); - force[c] -= on_fluid; + convective[c] -= rho * m_minus * (u_face - u0); + force[c] -= + mu * cv.ap[d][0] * a_d * (ud - u0) / solid_spacing(-1.0); } } } @@ -150,7 +218,7 @@ impl Mask { } } } - Some(force) + Some((force, convective)) } /// The closure lag of a moving body's pressure correction: the diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs index 2a70c17..def5ab9 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/reconstruct.rs @@ -52,7 +52,7 @@ impl Mask { continue; } let (k, j, i) = g.kji(idx); - if k < k0 || k >= k1 { + if k < k0 || k >= k1 || !super::exchange::in_load_window((i as f64 + 0.5) * g.dx) { continue; } let xc = [ diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index 9adba35..09eb6df 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -65,12 +65,19 @@ fn mode(s: f64) -> f64 { /// Centreline deflection and its velocity at arc parameter `s`, time `t`. fn deflection(s: f64, t: f64) -> (f64, f64) { let w = 2.0 * std::f64::consts::PI * FREQ; + let amp = amplitude(); ( - AMP * mode(s) * (w * t).sin(), - AMP * mode(s) * w * (w * t).cos(), + amp * mode(s) * (w * t).sin(), + amp * mode(s) * w * (w * t).cos(), ) } +/// The tip amplitude: `RTX_E3_FLAG_AMP` (default 0.084; 0 freezes the flag — +/// the static control of the load routes). +fn amplitude() -> f64 { + env_f("RTX_E3_FLAG_AMP", AMP) +} + /// Signed distance to the deflected flag's cross-section (a capsule /// around the centreline polyline of `n` segments) and the centreline's /// transverse velocity at the closest point. @@ -131,7 +138,10 @@ fn flag_wake_on_the_device() { let periods = env_f("RTX_E3_FLAG_PERIODS", 2.0); let h = H / ny as f64; let nx = (L / h).round() as usize; - let nz = ny; + // `RTX_E3_FLAG_NZ=4`: a thin slab periodic in z with the 2D inflow (Ū = 1) — the + // flag as a 2D problem, minutes per rung: the instrument for the load routes' parts. + let slab_nz = env_f("RTX_E3_FLAG_NZ", 0.0) as usize; + let nz = if slab_nz > 0 { slab_nz } else { ny }; let r_edge = h; let dt_cfl = 0.3 * h / (U_M.max(2.0 * std::f64::consts::PI * FREQ * AMP)); // `RTX_E3_FLAG_DT_SCALE` scales the step (the dt ladder of the loads). @@ -155,18 +165,36 @@ fn flag_wake_on_the_device() { tolerance: 1e-8, convection_scheme: ConvectionScheme::TvdVanAlbada, wall_scheme: WallScheme::CutCell, - boundaries: Boundaries { - x1: Side::PressureOutlet, - ..Boundaries::default() + boundaries: if slab_nz > 0 { + Boundaries { + x1: Side::PressureOutlet, + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + } + } else { + Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + } }, // The narrow band: the flag's tip speed bounds the surface motion. - max_surface_speed: Some(2.0 * std::f64::consts::PI * FREQ * AMP * 1.05), + max_surface_speed: Some( + (2.0 * std::f64::consts::PI * FREQ * amplitude() * 1.05).max(1e-3), + ), ..Parameters::default() }, ); - solver.set_boundary_velocity(|x, y, z, _t| { + let inflow_at = move |y: f64, z: f64| { + if slab_nz > 0 { + 6.0 * y * (H - y) / (H * H) + } else { + inflow(y, z) + } + }; + solver.set_boundary_velocity(move |x, y, z, _t| { if x <= 0.0 { - (inflow(y, z), 0.0, 0.0) + (inflow_at(y, z), 0.0, 0.0) } else { (0.0, 0.0, 0.0) } @@ -186,7 +214,7 @@ fn flag_wake_on_the_device() { let mut field = Field::new(g); for k in 0..nz { for j in 0..ny { - let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h); + let u0 = inflow_at((j as f64 + 0.5) * h, (k as f64 + 0.5) * h); for i in 0..=nx { field.u[g.uface(k, j, i)] = u0; } @@ -217,9 +245,17 @@ fn flag_wake_on_the_device() { let last_period_start = t_end - period; let mut next_phase = 0; let mid = nz / 2; - let slab = (mid - 2, mid + 2); + let slab = if slab_nz > 0 { + (0, nz) + } else { + (mid - 2, mid + 2) + }; + let width = nz as f64 * h; let start = std::time::Instant::now(); let mut drag_rec_sum = 0.0; + // The routes' PARTS over the whole body (x, per unit width): operator + // pressure / shear / exchange, reconstructed pressure / shear. + let mut parts = [[0.0_f64; 6]; 3]; let (mut drag_sum, mut lift_min, mut lift_max, mut samples) = (0.0, f64::INFINITY, f64::NEG_INFINITY, 0usize); let mut worst_residual = 0.0_f64; @@ -283,6 +319,30 @@ fn flag_wake_on_the_device() { .unwrap(); } if t >= last_period_start { + use rtx_cfd::solvers::incompressible::embedded3::exchange::set_load_window; + // whole body, the cylinder (x < 0.252), the flag + for (w, window) in [None, Some((0.0, 0.252)), Some((0.252, 10.0))] + .into_iter() + .enumerate() + { + set_load_window(window); + let (po, so) = mask + .cut_wall_force_parts(body, &field, RHO * NU, t) + .expect("parts"); + let (xd, xc) = mask + .cut_wall_exchange_parts(body, &field, RHO * NU, RHO, t, None) + .expect("exchange"); + let (pr, sr) = mask + .cut_wall_force_reconstructed_parts(body, &field, RHO * NU, t, None) + .expect("reconstructed parts"); + for (acc, v) in parts[w] + .iter_mut() + .zip([po[0], so[0], xd[0], xc[0], pr[0], sr[0]]) + { + *acc += v / width; + } + } + set_load_window(None); drag_sum += fs[0]; drag_rec_sum += fr[0]; lift_min = lift_min.min(fs[1]); @@ -306,6 +366,21 @@ fn flag_wake_on_the_device() { next_phase, start.elapsed().as_secs_f64() ); + let n = samples.max(1) as f64; + for (name, q) in ["whole body", "cylinder", "flag"].iter().zip(parts) { + println!( + " PARTS ny {ny} amp {:.3} {name} (x, N/m of width): operator pressure {:.2} + shear {:.2} + exchange diffusive {:.2} + convective {:.2} = {:.2}; reconstructed pressure {:.2} + shear {:.2} = {:.2}", + amplitude(), + q[0] / n, + q[1] / n, + q[2] / n, + q[3] / n, + (q[0] + q[1] + q[2] + q[3]) / n, + q[4] / n, + q[5] / n, + (q[4] + q[5]) / n + ); + } if let Some(t) = device.timers() { println!(" timers: {t:?}"); }