R7 composite: cut flag threshold, P2 direct comparison against the uniformly fine solution on the patch cells

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
2026-09-25 17:55:23 -05:00
parent e19cad03b2
commit d49f2787d9
2 changed files with 95 additions and 2 deletions
@@ -222,7 +222,7 @@ fn run_uniform(
                let v = frac[idx] * h * h * h;
                let inside = x.iter().all(|&c| c > region.0 && c < region.1);
                err[usize::from(!inside)].add(e, v);
-                if frac[idx] < 1.0 {
+                if frac[idx] < 1.0 - 1e-9 {
                    err[2].add(e, v);
                }
                err[3].add(e, v);
@@ -638,3 +638,96 @@ fn composite_quadratic_consistency() {
        }
    }
 }
+
+/// P2, direct: the composite solution on the patch's fine cells against
+/// the uniformly fine solution on the SAME cells with the same cut (the
+/// fine cut geometry of the patch vs the uniform grid's, compared too).
+#[test]
+#[ignore = "P2 composite vs uniform fine on the patch cells (minutes)"]
+fn composite_sphere_vs_fine() {
+    let mut csv = out_file("p2_vs_fine.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "scheme,n,max_dfrac,l2_diff_patch,linf_diff_patch,l2_diff_cut,linf_diff_cut,l2_err_fine_patch,linf_err_fine_patch"
+        )
+        .ok();
+    }
+    let body = sphere();
+    for &n in &[16usize, 32, 64] {
+        let m = 2 * n;
+        let gf = Grid::cubic(m, m, m, 1.0 / m as f64);
+        let (prob, frac) = uniform_problem(gf, Some(&body), &sph_f, &sph_u);
+        let mut pf = vec![0.0; gf.cells()];
+        let l1: f64 = prob.rhs.iter().map(|v| v.abs()).sum();
+        let params = MultigridParameters {
+            max_iterations: 2000,
+            ..MultigridParameters::default()
+        };
+        assert!(solve_pcg(&prob, &mut pf, &params, 1e-12 * l1, None).converged);
+        for iface in [Interface::Octree, Interface::Quadratic] {
+            let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+            let spec = CompositeSpec {
+                coarse: g,
+                lo: [n / 4; 3],
+                hi: [3 * n / 4; 3],
+                interface: iface,
+                body: Some(body.clone()),
+                source: &sph_f,
+                dirichlet: &sph_u,
+            };
+            let c = Composite::build(&spec);
+            let mut x = vec![0.0; c.unknowns()];
+            assert!(solve_bicgstab(&c, &mut x, 1e-12, 400, 2).converged);
+            let (mut d_all, mut d_cut, mut e_fine) =
+                (Err::default(), Err::default(), Err::default());
+            let mut max_dfrac: f64 = 0.0;
+            let off = n / 2; // the patch's first fine index on the uniform fine grid
+            for k in 0..c.fine.nz {
+                for j in 0..c.fine.ny {
+                    for i in 0..c.fine.nx {
+                        let u = c.fine_id[c.fine.cell(k, j, i)];
+                        let idx = gf.cell(k + off, j + off, i + off);
+                        if u == usize::MAX {
+                            assert!(!prob.active[idx] || frac[idx] < 1e-12, "activity differs");
+                            continue;
+                        }
+                        let hf = c.fine.dx;
+                        let fr = c.volume[u] / (hf * hf * hf);
+                        max_dfrac = max_dfrac.max((fr - frac[idx]).abs());
+                        let d = x[u] - pf[idx];
+                        d_all.add(d, c.volume[u]);
+                        if c.cut[u] {
+                            d_cut.add(d, c.volume[u]);
+                        }
+                        let p = c.centre[u];
+                        e_fine.add(pf[idx] - sph_u(p[0], p[1], p[2]), c.volume[u]);
+                    }
+                }
+            }
+            eprintln!(
+                "{iface:?} n {n} vs uniform {m}: max |Δfrac| {max_dfrac:.1e}; composite − fine on the patch L2 {:.3e} Linf {:.3e}; on cut cells ({}) L2 {:.3e} Linf {:.3e}; the fine grid's own error there L2 {:.3e} Linf {:.3e}",
+                d_all.l2(),
+                d_all.max,
+                d_cut.count,
+                d_cut.l2(),
+                d_cut.max,
+                e_fine.l2(),
+                e_fine.max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{iface:?},{n},{max_dfrac:.3e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e}",
+                    d_all.l2(),
+                    d_all.max,
+                    d_cut.l2(),
+                    d_cut.max,
+                    e_fine.l2(),
+                    e_fine.max
+                )
+                .ok();
+            }
+        }
+    }
+}