embedded3 S2-2b-ii: DeviceCg::refresh (fine operator per step, hierarchy every K steps, z zeroed before the scatter) + gate test; the flag driver's correctors 3 / inner 1e-3 (residual remedy measured 7.6e-9 on the moving circle); residual study test

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

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs

@@ -177,6 +177,73 @@ impl DeviceCg {
         self.key.matches(problem, params)
     }
 
+    /// Refresh the FINE operator (cells, coefficients, links, components)
+    /// for `problem`, keeping the multigrid hierarchy as it was: the CG's
+    /// operator is exact, the preconditioner stale (a moving body's
+    /// operator changes little per step; the hierarchy's rebuild is the
+    /// cost). `z` is zeroed before every V-cycle scatter, so cells absent
+    /// from the stale hierarchy get no correction rather than a stale one.
+    pub fn refresh(&mut self, problem: &Problem, params: &MultigridParameters) {
+        let rt = runtime();
+        let fine = Level::<f64>::new(problem.clone());
+        let components = Components::find(problem, &fine.cells);
+        let singular_count = components.singular.iter().filter(|&&s| s).count();
+        assert!(
+            singular_count == 0 || components.members.len() == 1,
+            "DeviceCg::refresh: {} components with {} singular",
+            components.members.len(),
+            singular_count
+        );
+        let to_u32 = |v: &[usize]| {
+            v.iter()
+                .map(|&i| if i == usize::MAX { u32::MAX } else { i as u32 })
+                .collect::<Vec<u32>>()
+        };
+        let up_u = |v: &[u32]| -> CudaSlice<u32> {
+            rt.stream
+                .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v })
+                .expect("upload")
+        };
+        let up_f = |v: &[f64]| -> CudaSlice<f64> { rt.stream.memcpy_stod(v).expect("upload") };
+        let lists = problem.link_lists();
+        let mut link_ptr = Vec::with_capacity(self.n + 1);
+        let mut link_idx = Vec::new();
+        let mut link_coef = Vec::new();
+        link_ptr.push(0u32);
+        for list in &lists {
+            for &(other, c) in list {
+                link_idx.push(other as u32);
+                link_coef.push(c);
+            }
+            link_ptr.push(link_idx.len() as u32);
+        }
+        self.n_cells = fine.cells.len();
+        self.n_blocks = self.n_cells.div_ceil(256).max(1);
+        if self.partial.len() < self.n_blocks {
+            self.partial = rt.stream.alloc_zeros::<f64>(self.n_blocks).expect("alloc");
+        }
+        self.cells = up_u(&to_u32(&fine.cells));
+        self.top = up_u(&to_u32(&fine.top));
+        self.bot = up_u(&to_u32(&fine.bot));
+        self.ae = up_f(&fine.ae);
+        self.aw = up_f(&fine.aw);
+        self.an = up_f(&fine.an);
+        self.as_ = up_f(&fine.as_);
+        self.at = up_f(&fine.at);
+        self.ab = up_f(&fine.ab);
+        self.ap = up_f(&fine.ap);
+        self.link_ptr = up_u(&link_ptr);
+        self.link_idx = up_u(&link_idx);
+        self.link_coef = up_f(if link_coef.is_empty() {
+            &[0.0]
+        } else {
+            &link_coef
+        });
+        self.singular = singular_count > 0;
+        self.active_host = fine.active.clone();
+        self.key = OperatorKey::of(problem, params);
+    }
+
     pub fn n_cells(&self) -> usize {
         self.n_cells
     }
@@ -367,6 +434,7 @@ impl DeviceCg {
         let rt = runtime();
         let k = kernels();
         let n_i = self.n_cells as i32;
+        rt.stream.memset_zeros(&mut self.z).expect("z = 0");
         unsafe {
             rt.stream
                 .launch_builder(&k.gather_f32)

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs

@@ -155,6 +155,8 @@ pub struct DeviceStep {
     initialized: bool,
     /// A static cut-cell mask's tables (item 9b), when the solver has one.
     cut: Option<cut::DeviceCut>,
+    /// Steps since the multigrid hierarchy was last rebuilt (moving bodies).
+    steps_since_hierarchy: usize,
 }
 
 impl DeviceStep {
@@ -210,6 +212,7 @@ impl DeviceStep {
             timers,
             initialized: false,
             cut,
+            steps_since_hierarchy: 0,
         }
     }
 

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs

@@ -329,7 +329,26 @@ impl DeviceStep {
             fresh_cells = self.solver.rebuild_moving_mask(&mut field, dt, t_new);
             self.upload(&field);
             self.cut = DeviceCut::build(&self.solver, g, Phase::Projection, t_new);
-            self.cg = None;
+            // The operator: refreshed every step, the hierarchy every
+            // `RTX_E3_PRECOND_REFRESH` steps (default 10; 1 = rebuild always).
+            let every: usize = std::env::var("RTX_E3_PRECOND_REFRESH")
+                .ok()
+                .and_then(|v| v.parse().ok())
+                .unwrap_or(10)
+                .max(1);
+            self.steps_since_hierarchy += 1;
+            if self.steps_since_hierarchy >= every || self.cg_dt != dt {
+                self.cg = None;
+                self.steps_since_hierarchy = 0;
+            } else if let Some(cg) = self.cg.as_mut() {
+                let problem = self.solver.poisson_operator(g, dt);
+                let params = MultigridParameters {
+                    precision: self.solver.params.poisson_precision,
+                    smoother: self.solver.params.poisson_smoother,
+                    ..MultigridParameters::default()
+                };
+                cg.refresh(&problem, &params);
+            }
             rt.stream
                 .memcpy_dtod(&self.u, &mut self.u_star)
                 .expect("u*");

crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs

@@ -99,7 +99,9 @@ fn march(tight: bool, steps: usize, bound: Option<f64>) {
     let sp = scale(&fh.p).max(1000.0);
     let t = device.timers().expect("timers");
     println!(
-        "  moving circle 152²×4 CutCell (tight {tight}): {steps} steps; host vs device max |Δu| {du:.3e} on {su:.3e}, max |Δp| {dp:.3e} on {sp:.3e}; fresh cells host {fresh_h} device {fresh_d}; corrector counts differ on {differ} steps; {:.1} ms per step (host + device), device step {:.1} ms of which rebuild {:.1} ms",
+        "  moving circle 152²×4 CutCell (tight {tight}, hierarchy every {} steps): {steps} steps; host vs device max |Δu| {du:.3e} on {su:.3e}, max |Δp| {dp:.3e} on {sp:.3e}; fresh cells host {fresh_h} device {fresh_d}; corrector counts differ on {differ} steps; device CG {:.1}/step; {:.1} ms per step (host + device), device step {:.1} ms of which rebuild {:.1} ms",
+        std::env::var("RTX_E3_PRECOND_REFRESH").unwrap_or_else(|_| "10".into()),
+        t.cg_iterations as f64 / steps as f64,
         1e3 * seconds / steps as f64,
         1e-6 * (t.predictor_ns + t.poisson_ns + t.apply_ns + t.transfer_ns) as f64 / steps as f64,
         1e-6 * t.transfer_ns as f64 / steps as f64
@@ -116,3 +118,16 @@ fn moving_circle_host_equals_device() {
     march(false, 100, None);
     march(true, 100, Some(1e-9));
 }
+
+/// S2-2b-ii: the hierarchy kept for ten steps (the fine operator exact
+/// every step) — the tight identity to the host still holds (the CG
+/// converges to the tolerance under any preconditioner), the CG count
+/// and the step time recorded against the every-step rebuild.
+#[test]
+fn stale_hierarchy_keeps_the_identity() {
+    unsafe { std::env::set_var("RTX_E3_PRECOND_REFRESH", "1") };
+    march(true, 60, Some(1e-9));
+    unsafe { std::env::set_var("RTX_E3_PRECOND_REFRESH", "10") };
+    march(true, 60, Some(1e-9));
+    unsafe { std::env::remove_var("RTX_E3_PRECOND_REFRESH") };
+}

crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs

@@ -249,3 +249,68 @@ fn operator_export_profile_at_ny_62() {
         problem.links.len()
     );
 }
+
+/// The residual clause's remedy, measured on the moving circle (host):
+/// the worst mass residual over 100 steps with two correctors at the
+/// 1e-2 inner stop against three correctors at 1e-3.
+#[test]
+#[ignore = "residual study on the moving circle (host, a minute)"]
+fn moving_circle_residual_remedy() {
+    use rtx_cfd::solvers::incompressible::ConvectionScheme;
+    use rtx_cfd::solvers::incompressible::embedded3::{
+        Field, Fluid, Parameters, Solver, WallScheme,
+    };
+    let n = 152;
+    let h = 1.0 / n as f64;
+    let dt = 3.24e-4;
+    for (correctors, inner) in [(2usize, 1e-2), (3, 1e-3), (4, 1e-4)] {
+        let mut solver = Solver::new(
+            Fluid {
+                density: 1000.0,
+                viscosity: 1.0,
+                reference_velocity: 1.0,
+                reference_length: 0.1,
+            },
+            Parameters {
+                corrector_steps: correctors,
+                tolerance: 1e-8,
+                inner_stop_factor: inner,
+                convection_scheme: ConvectionScheme::Upwind,
+                wall_scheme: WallScheme::CutCell,
+                boundaries: Boundaries {
+                    z0: Side::Periodic,
+                    z1: Side::Periodic,
+                    ..Boundaries::default()
+                },
+                max_surface_speed: Some(1.0),
+                ..Parameters::default()
+            },
+        );
+        solver.set_boundary_velocity(|_, _, _, _| (0.0, 0.0, 0.0));
+        let yc = |t: f64| 0.5 + 0.08 * (t / 0.08).sin();
+        let vc = |t: f64| (t / 0.08).cos();
+        solver.set_moving_body(
+            Body::from_sdf(move |x, y, _z, t| {
+                ((x - 0.5_f64).powi(2) + (y - yc(t)).powi(2)).sqrt() - 0.05
+            })
+            .with_surface_velocity(move |_, _, _, t| (0.0, vc(t), 0.0)),
+        );
+        let g = Grid::cubic(n, n, 4, h);
+        let mut f = Field::new(g);
+        solver.initialize(&mut f);
+        let (mut worst, mut sum_cg, mut sum_corr) = (0.0_f64, 0usize, 0usize);
+        let start = std::time::Instant::now();
+        for _ in 0..100 {
+            let r = solver.advance(&mut f, dt);
+            worst = worst.max(r.final_residual);
+            sum_cg += r.poisson_iterations;
+            sum_corr += r.corrector_steps_performed;
+        }
+        println!(
+            "  correctors {correctors} inner {inner:.0e}: worst residual {worst:.2e}, CG {:.1}/step, correctors {:.2}/step, {:.2} s",
+            sum_cg as f64 / 100.0,
+            sum_corr as f64 / 100.0,
+            start.elapsed().as_secs_f64()
+        );
+    }
+}

crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs

@@ -141,7 +141,12 @@ fn flag_wake_on_the_device() {
             reference_length: 2.0 * R_CYL,
         },
         Parameters {
-            corrector_steps: 2,
+            // Three correctors at a 1e-3 inner stop hold the moving cut
+            // wall's mass residual under 1e-8 (the moving circle: 7.6e-9
+            // against 1.5e-6 with two at 1e-2); `RTX_E3_FLAG_CORRECTORS`
+            // overrides for the comparison runs.
+            corrector_steps: env_f("RTX_E3_FLAG_CORRECTORS", 3.0) as usize,
+            inner_stop_factor: env_f("RTX_E3_FLAG_INNER", 1e-3),
             tolerance: 1e-8,
             convection_scheme: ConvectionScheme::TvdVanAlbada,
             wall_scheme: WallScheme::CutCell,