diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs index f6aee3d..c4ee0ec 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/device_cg.rs +++ b/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::::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::>() + }; + let up_u = |v: &[u32]| -> CudaSlice { + rt.stream + .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v }) + .expect("upload") + }; + let up_f = |v: &[f64]| -> CudaSlice { 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::(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) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs index 14a813e..4cb3cf9 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs +++ b/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, + /// 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, } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs index 941d442..57d28d9 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs +++ b/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, ¶ms); + } rt.stream .memcpy_dtod(&self.u, &mut self.u_star) .expect("u*"); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs b/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs index a0807cb..755eaf7 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_device_moving.rs @@ -99,7 +99,9 @@ fn march(tight: bool, steps: usize, bound: Option) { 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") }; +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs index 87864cd..c777f89 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_geometry.rs +++ b/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() + ); + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index 3c041f4..2e6b43d 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/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,