From 338d66f85eac4e7696d7ccaad5554f61cd399cbd Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 17:24:38 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20item=2012:=20the=20DFG=203D-2Z=20de?= =?UTF-8?q?vice=20driver=20(loads=20by=20both=20routes,=20=CE=94p,=20settl?= =?UTF-8?q?e=20criterion,=20CSV)=20and=20the=20VTK=20instant=20export?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/export_vtk.rs | 55 +++++ .../solvers/incompressible/embedded3/mod.rs | 2 + .../rtx-cfd/tests/embedded3_dfg_2z.rs | 189 ++++++++++++++++++ 3 files changed, 246 insertions(+) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/export_vtk.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/export_vtk.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/export_vtk.rs new file mode 100644 index 0000000..2510d5d --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/export_vtk.rs @@ -0,0 +1,55 @@ +//! VTK legacy `STRUCTURED_POINTS` export of an instant: cell-centred +//! velocity (face values averaged), pressure, and the fluid fraction / +//! mask (1 fluid, 0 solid) — for ParaView and the viewer's 3D wake. + +use super::field::Field; +use super::wall::Mask; +use std::io::{BufWriter, Write}; + +/// Write `field` (and the mask's fluid fraction when given) at `path`. +pub fn write_vtk( + path: &std::path::Path, + field: &Field, + mask: Option<&Mask>, +) -> std::io::Result<()> { + let g = field.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let mut out = BufWriter::new(std::fs::File::create(path)?); + writeln!(out, "# vtk DataFile Version 3.0")?; + writeln!(out, "embedded3 instant")?; + writeln!(out, "ASCII")?; + writeln!(out, "DATASET STRUCTURED_POINTS")?; + writeln!(out, "DIMENSIONS {nx} {ny} {nz}")?; + writeln!(out, "ORIGIN {} {} {}", 0.5 * g.dx, 0.5 * g.dy, 0.5 * g.dz)?; + writeln!(out, "SPACING {} {} {}", g.dx, g.dy, g.dz)?; + writeln!(out, "POINT_DATA {}", g.cells())?; + writeln!(out, "VECTORS velocity double")?; + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let uc = 0.5 * (field.u[g.uface(k, j, i)] + field.u[g.uface(k, j, i + 1)]); + let vc = 0.5 * (field.v[g.vface(k, j, i)] + field.v[g.vface(k, j + 1, i)]); + let wc = 0.5 * (field.w[g.wface(k, j, i)] + field.w[g.wface(k + 1, j, i)]); + writeln!(out, "{uc:.6e} {vc:.6e} {wc:.6e}")?; + } + } + } + writeln!(out, "SCALARS pressure double 1")?; + writeln!(out, "LOOKUP_TABLE default")?; + for &p in &field.p { + writeln!(out, "{p:.6e}")?; + } + if let Some(m) = mask { + writeln!(out, "SCALARS fluid double 1")?; + writeln!(out, "LOOKUP_TABLE default")?; + for idx in 0..g.cells() { + let f = if m.is_fluid_cell(idx) { + m.vol(idx) + } else { + 0.0 + }; + writeln!(out, "{f:.4}")?; + } + } + out.flush() +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs index c83cdb9..c46facf 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -9,6 +9,7 @@ pub mod body; pub mod cut; pub mod cutwall; +pub mod export_vtk; pub mod field; pub mod grid; pub mod impose; @@ -19,6 +20,7 @@ pub mod wall; pub use body::{Body, SurfaceSample}; pub use cut::CutGeometry; +pub use export_vtk::write_vtk; pub use field::Field; pub use grid::Grid; pub use loads::SurfaceForce; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs new file mode 100644 index 0000000..622b894 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs @@ -0,0 +1,189 @@ +//! embedded3 item 12: Schäfer–Turek DFG 3D-2Z (the laminar cylinder in the +//! square channel, Re 20, steady) on the device cut-cell wall. Channel +//! 2.5 × 0.41 × 0.41, cylinder D 0.1 at (0.5, 0.2) across the width, +//! inflow U(y, z) = 16 U_m y z (H − y)(H − z)/H⁴ with U_m 0.45 (Ū = 0.2), +//! ρ 1, ν 1e-3. Coefficients `c = 2F/(ρ Ū² D H)`, Δp between the front +//! and back stagnation points. Reference (the DFG bar): c_D 6.05–6.25, +//! c_L 0.008–0.010, Δp 0.165–0.175. +//! +//! `RTX_E3_DFG_NY=62 RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_dfg_2z -- --ignored --nocapture` +//! `RTX_E3_DFG_T` sets the flow time marched (default 10 s); `RTX_E3_DFG_VTK=` +//! writes the final instant; `RTX_E3_DFG_CSV=` the load history. +#![cfg(feature = "cuda")] + +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, write_vtk, +}; +use std::io::Write as _; + +const H: f64 = 0.41; +const L: f64 = 2.5; +const D: f64 = 0.1; +const CX: f64 = 0.5; +const CY: f64 = 0.2; +const U_M: f64 = 0.45; +const U_BAR: f64 = 4.0 / 9.0 * U_M; +const RHO: f64 = 1.0; +const NU: f64 = 1e-3; + +fn env_f(name: &str, default: f64) -> f64 { + std::env::var(name) + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(default) +} + +fn inflow(y: f64, z: f64) -> f64 { + 16.0 * U_M * y * z * (H - y) * (H - z) / (H * H * H * H) +} + +#[test] +#[ignore = "item 12: the DFG rung on the device (minutes at ny 62, hours at ny 123)"] +fn dfg_3d_2z_on_the_device() { + let ny = env_f("RTX_E3_DFG_NY", 62.0) as usize; + let h = H / ny as f64; + let nx = (L / h).round() as usize; + let nz = ny; + let t_end = env_f("RTX_E3_DFG_T", 10.0); + // Explicit stability: CFL 0.3 on U_m and half the viscous limit. + let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU)); + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: RHO * NU, + reference_velocity: U_BAR, + reference_length: D, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + convection_scheme: ConvectionScheme::TvdVanAlbada, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|x, y, z, _t| { + if x <= 0.0 { + (inflow(y, z), 0.0, 0.0) + } else { + (0.0, 0.0, 0.0) + } + }); + solver.set_body(Body::cylinder_z(CX, CY, 0.5 * D)); + let g = Grid::cubic(nx, ny, nz, h); + let mut field = Field::new(g); + // Start from the inflow profile everywhere (a faster approach to steady). + 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); + for i in 0..=nx { + field.u[g.uface(k, j, i)] = u0; + } + } + } + solver.initialize(&mut field); + let mask_cells = solver.mask().map_or(0, |m| m.fluid_cells()); + println!( + " DFG 3D-2Z ny {ny}: {nx}×{ny}×{nz} = {} cells ({mask_cells} fluid), h {h:.4e}, dt {dt:.3e}, t_end {t_end}", + g.cells() + ); + let mut device = DeviceStep::new(solver, g); + device.upload(&field); + let steps = (t_end / dt).ceil() as usize; + let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H); + let csv = std::env::var("RTX_E3_DFG_CSV").ok().map(|p| { + let mut f = std::fs::File::create(p).expect("csv"); + writeln!(f, "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg").unwrap(); + f + }); + let mut csv = csv; + let sample_every = (steps / 100).max(1); + let start = std::time::Instant::now(); + let mut last: Option<(f64, f64, f64, f64, f64)> = None; + let mut settled = false; + for step in 0..steps { + let r = device.advance(dt); + if (step + 1) % sample_every == 0 || step + 1 == steps { + device.download(&mut field); + let solver = &device.solver; + let mask = solver.mask().expect("mask"); + let body = solver.body().expect("body"); + let t = solver.time(); + let fw = mask + .cut_wall_force(body, &field, RHO * NU, t) + .expect("wall"); + let margin = 3.0 * D; + let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2); + let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2); + let bx = ( + ci(CX - margin), + ci(CX + margin), + cj(CY - 0.15), + cj(CY + 0.15), + 0, + nz, + ); + let fcv = mask.control_volume_force(&field, dt, RHO, RHO * NU, None, bx); + let zc = 0.5 * H; + let p_front = mask + .pressure_at(&field.p, CX - 0.5 * D, CY, zc) + .unwrap_or(f64::NAN); + let p_back = mask + .pressure_at(&field.p, CX + 0.5 * D, CY, zc) + .unwrap_or(f64::NAN); + let dp = p_front - p_back; + let (cd, cl, cd_cv, cl_cv) = (coef * fw[0], coef * fw[1], coef * fcv[0], coef * fcv[1]); + println!( + " t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}) c_L {cl:.5} (CV {cl_cv:.5}) Δp {dp:.4} residual {:.1e} CG {} [{:.0} s]", + r.final_residual, + r.poisson_iterations, + start.elapsed().as_secs_f64() + ); + if let Some(f) = csv.as_mut() { + writeln!( + f, + "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{}", + r.final_residual, r.poisson_iterations + ) + .unwrap(); + } + if let Some((pcd, pcl, _, _, pdp)) = last { + let rel = ((cd - pcd) / cd) + .abs() + .max(((dp - pdp) / dp).abs()) + .max((cl - pcl).abs() / 0.01); + if rel < 1e-4 && t > 2.0 { + settled = true; + } + } + last = Some((cd, cl, cd_cv, cl_cv, dp)); + if settled { + println!(" settled (relative change < 1e-4 between samples)"); + break; + } + } + } + device.download(&mut field); + let solver = &device.solver; + let mask = solver.mask().expect("mask"); + if let Ok(dir) = std::env::var("RTX_E3_DFG_VTK") { + let path = std::path::Path::new(&dir).join(format!("dfg_2z_ny{ny}.vtk")); + write_vtk(&path, &field, Some(mask)).expect("vtk"); + println!(" instant written to {}", path.display()); + } + let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples"); + println!( + " FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart) c_L {cl:.5} (CV {cl_cv:.5}) Δp {dp:.4} — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s", + ((cd - cd_cv) / cd).abs(), + start.elapsed().as_secs_f64() + ); + if let Some(t) = device.timers() { + println!(" timers: {t:?}"); + } +}