embedded3 item 12: the DFG 3D-2Z device driver (loads by both routes, Δp, settle criterion, CSV) and the VTK instant export

Co-Authored-By: Claude Fable 5.1 <[email protected]>

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()
+}

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;

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=<dir>`
+//! writes the final instant; `RTX_E3_DFG_CSV=<path>` 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:?}");
+    }
+}