rustytorch/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs

//! embedded3 gates 9a and 10: the manufactured solution with an embedded
//! sphere (centre (0.6, 0.45, 0.5), r 0.2, off-centre so the exact force is
//! not zero by symmetry) carrying the exact field as its surface velocity,
//! on the binary ghost wall (item 9) and the apertured cut-cell wall (item
//! 10). The velocity error falls at the scheme's order, every fluid cell
//! is divergence-free (apertured, with the porous surface's flux, on the
//! cut wall), the compatibility correction shrinks, and both load routes
//! converge to the exact surface integral of the manufactured stress (the
//! control-volume route measures F − M with M the momentum flux through
//! the porous manufactured surface). Item 10's gate: the cut wall's errors
//! are at most the binary wall's at every n, its loads within 10 % at the
//! finest rung.

use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, FaceKind, Field, Fluid, Grid, Parameters, Solver, WallScheme,
};
use std::f64::consts::PI;

const RHO: f64 = 1.0;
const MU: f64 = 0.05;
const C: (f64, f64, f64) = (0.6, 0.45, 0.5);
const R: f64 = 0.2;

fn u3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos() * (PI * z).cos()
}
fn v3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).cos() * (PI * y).sin() * (PI * z).cos()
}
fn w3(x: f64, y: f64, z: f64) -> f64 {
    -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin()
}
fn p3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin() * (PI * z).sin()
}
/// The velocity gradient ∂u_i/∂x_j and the pressure gradient.
fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) {
    let (sx, cx) = (PI * x).sin_cos();
    let (sy, cy) = (PI * y).sin_cos();
    let (sz, cz) = (PI * z).sin_cos();
    (
        [
            [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz],
            [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz],
            [
                2.0 * PI * sx * cy * sz,
                2.0 * PI * cx * sy * sz,
                -2.0 * PI * cx * cy * cz,
            ],
        ],
        [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz],
    )
}
fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let (g, gp) = grads(x, y, z);
    let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
    let lap = -3.0 * PI * PI;
    let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2];
    (
        RHO * conv(0) + gp[0] - MU * lap * u[0],
        RHO * conv(1) + gp[1] - MU * lap * u[1],
        RHO * conv(2) + gp[2] - MU * lap * u[2],
    )
}
fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let u = if x <= 0.0 || x >= 1.0 {
        0.0
    } else {
        u3(x, y, z)
    };
    let v = if y <= 0.0 || y >= 1.0 {
        0.0
    } else {
        v3(x, y, z)
    };
    let w = if z <= 0.0 || z >= 1.0 {
        0.0
    } else {
        w3(x, y, z)
    };
    (u, v, w)
}

/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA`
/// over the sphere by a fine Fibonacci quadrature.
fn exact_force_and_flux() -> ([f64; 3], [f64; 3]) {
    let n = 200_000;
    let golden = PI * (3.0 - 5.0_f64.sqrt());
    let (mut f, mut m) = ([0.0; 3], [0.0; 3]);
    let da = 4.0 * PI * R * R / n as f64;
    for k in 0..n {
        let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64;
        let rr = (1.0 - zz * zz).sqrt();
        let th = golden * k as f64;
        let nrm = [rr * th.cos(), rr * th.sin(), zz];
        let (x, y, z) = (C.0 + R * nrm[0], C.1 + R * nrm[1], C.2 + R * nrm[2]);
        let (g, _) = grads(x, y, z);
        let p = p3(x, y, z);
        let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
        let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2];
        for i in 0..3 {
            let mut t = -p * nrm[i];
            for j in 0..3 {
                t += MU * (g[i][j] + g[j][i]) * nrm[j];
            }
            f[i] += t * da;
            m[i] += RHO * u[i] * un * da;
        }
    }
    (f, m)
}

struct Measurement {
    l2_velocity: f64,
    max_div: f64,
    ghost_correction: f64,
    force_surface: [f64; 3],
    skipped: usize,
    force_cv: [f64; 3],
}

fn measure(n: usize, scheme: WallScheme) -> Measurement {
    let h = 1.0 / n as f64;
    let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h);
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            convection_scheme: ConvectionScheme::Upwind,
            wall_scheme: scheme,
            ..Parameters::default()
        },
    );
    solver.set_momentum_source(|x, y, z, _t| source3(x, y, z));
    solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z));
    solver.set_body(
        Body::sphere(|_t| C, R)
            .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))),
    );
    let g = Grid::cubic(n, n, n, h);
    let mut f = Field::new(g);
    solver.initialize(&mut f);
    let mut last = solver.advance(&mut f, dt);
    for _ in 0..200_000 {
        let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone());
        last = solver.advance(&mut f, dt);
        let mut change = 0.0_f64;
        for (a, b) in
            f.u.iter()
                .zip(&bu)
                .chain(f.v.iter().zip(&bv))
                .chain(f.w.iter().zip(&bw))
        {
            change = change.max((a - b).abs());
        }
        if change / dt < 1e-6 {
            break;
        }
    }
    let mask = solver.mask().expect("mask");
    let (mut sq, mut vol) = (0.0, 0.0);
    let dv = h * h * h;
    for k in 0..n {
        for j in 0..n {
            for i in 1..n {
                if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid {
                    let e = f.u[g.uface(k, j, i)]
                        - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
        for j in 1..n {
            for i in 0..n {
                if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid {
                    let e = f.v[g.vface(k, j, i)]
                        - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    for k in 1..n {
        for j in 0..n {
            for i in 0..n {
                if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid {
                    let e = f.w[g.wface(k, j, i)]
                        - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    let body = solver.body().expect("body");
    let t = solver.time();
    // The apertured divergence per unit volume, the porous surface's flux
    // through the wall included (the plain divergence on the binary wall).
    let mut max_div = 0.0_f64;
    let mut at_vol = 1.0;
    let mut sum_flux = 0.0;
    let (wall_fluxes, _) = mask.wall_flux_table(body, t);
    for k in 0..n {
        for j in 0..n {
            for i in 0..n {
                let idx = g.cell(k, j, i);
                if mask.is_fluid_cell(idx) {
                    let flux = (mask.a_u(g.uface(k, j, i + 1)) * f.u[g.uface(k, j, i + 1)]
                        - mask.a_u(g.uface(k, j, i)) * f.u[g.uface(k, j, i)])
                        * h
                        * h
                        + (mask.a_v(g.vface(k, j + 1, i)) * f.v[g.vface(k, j + 1, i)]
                            - mask.a_v(g.vface(k, j, i)) * f.v[g.vface(k, j, i)])
                            * h
                            * h
                        + (mask.a_w(g.wface(k + 1, j, i)) * f.w[g.wface(k + 1, j, i)]
                            - mask.a_w(g.wface(k, j, i)) * f.w[g.wface(k, j, i)])
                            * h
                            * h
                        + wall_fluxes[idx];
                    sum_flux += flux.abs();
                    if (flux / (h * h * h)).abs() > max_div {
                        max_div = (flux / (h * h * h)).abs();
                        at_vol = mask.vol(idx);
                    }
                }
            }
        }
    }
    println!(
        "    [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}",
        last.final_residual
    );
    let surface = match scheme {
        WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h),
        WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce {
            f: mask.cut_wall_force(body, &f, MU, t).expect("cut wall"),
            samples: 0,
            skipped: 0,
        },
    };
    let (i0, i1) = (n / 8, n - n / 8);
    let src = |x: f64, y: f64, z: f64| source3(x, y, z);
    let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1));
    Measurement {
        l2_velocity: (sq / vol).sqrt(),
        max_div,
        ghost_correction: solver.ghost_correction().abs(),
        force_surface: surface.f,
        skipped: surface.skipped,
        force_cv,
    }
}

fn norm(a: [f64; 3]) -> f64 {
    (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
}

/// The velocity errors and the two routes' relative force errors per rung.
struct Ladder {
    errors: Vec<f64>,
    surface: Vec<f64>,
    cv: Vec<f64>,
}

fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder {
    let (fe, m) = exact_force_and_flux();
    let f_scale = norm(fe);
    let fcv = [fe[0] - m[0], fe[1] - m[1], fe[2] - m[2]];
    println!(
        "  {scheme:?}: exact force {fe:.5?}; momentum flux {m:.5?}; the control-volume route measures {fcv:.5?}"
    );
    let ms: Vec<Measurement> = resolutions.iter().map(|&n| measure(n, scheme)).collect();
    let errors: Vec<f64> = ms.iter().map(|x| x.l2_velocity).collect();
    let mut se = Vec::new();
    let mut ce = Vec::new();
    for (k, (mm, &n)) in ms.iter().zip(resolutions).enumerate() {
        let rate = if k == 0 {
            "    -".to_string()
        } else {
            format!("{:5.2}", (errors[k - 1] / errors[k]).log2())
        };
        let s = norm([
            mm.force_surface[0] - fe[0],
            mm.force_surface[1] - fe[1],
            mm.force_surface[2] - fe[2],
        ]) / f_scale;
        let c = norm([
            mm.force_cv[0] - fcv[0],
            mm.force_cv[1] - fcv[1],
            mm.force_cv[2] - fcv[2],
        ]) / f_scale;
        println!(
            "  n = {n:3}  L2 u {:.4e} (order {rate})  max div {:.2e}  ghost corr {:.2e}  F_surface {:.4?} rel {s:.3e} (skipped {})  F_cv {:.4?} rel {c:.3e}",
            mm.l2_velocity,
            mm.max_div,
            mm.ghost_correction,
            mm.force_surface,
            mm.skipped,
            mm.force_cv
        );
        se.push(s);
        ce.push(c);
    }
    assert!(
        errors.windows(2).all(|w| w[1] < w[0]),
        "errors not monotone {errors:?}"
    );
    for w in errors.windows(2) {
        let rate = (w[0] / w[1]).log2();
        assert!(
            rate > 0.75 && rate < 2.3,
            "order {rate:.3} outside [0.75, 2.3]"
        );
    }
    for mm in &ms {
        assert!(mm.max_div < 1e-5, "max div {:.3e}", mm.max_div);
    }
    assert!(
        se.windows(2).all(|w| w[1] < w[0]),
        "surface-route error not falling {se:?}"
    );
    assert!(
        ce.windows(2).all(|w| w[1] < w[0]),
        "control-volume-route error not falling {ce:?}"
    );
    Ladder {
        errors,
        surface: se,
        cv: ce,
    }
}

/// Item 10's comparison: the cut wall's velocity error at most the binary
/// wall's at every rung; both routes within `load_bound` at the finest.
fn compare(resolutions: &[usize], load_bound: f64) {
    let ghost = ladder(resolutions, WallScheme::GhostBinary);
    let cut = ladder(resolutions, WallScheme::CutCell);
    for (k, &n) in resolutions.iter().enumerate() {
        println!(
            "  n = {n:3}  L2 u ghost {:.4e} cut {:.4e}  (ratio {:.3})",
            ghost.errors[k],
            cut.errors[k],
            cut.errors[k] / ghost.errors[k]
        );
        assert!(
            cut.errors[k] <= ghost.errors[k],
            "cut-cell error above the binary wall's at n = {n}"
        );
    }
    let last = resolutions.len() - 1;
    assert!(
        cut.surface[last] < load_bound && cut.cv[last] < load_bound,
        "cut-cell loads at the finest rung: surface {:.3e}, control volume {:.3e} (bound {load_bound})",
        cut.surface[last],
        cut.cv[last]
    );
}

#[test]
fn embedded_sphere_recovers_the_manufactured_solution() {
    ladder(&[12, 24], WallScheme::GhostBinary);
}

#[test]
fn cut_cell_wall_recovers_the_manufactured_solution() {
    compare(&[12, 24], 0.2);
}

#[test]
#[ignore = "the three-rung ladder to n = 48 (minutes on the host)"]
fn embedded_sphere_three_rungs() {
    ladder(&[12, 24, 48], WallScheme::GhostBinary);
}

#[test]
#[ignore = "item 10's finest rung: the cut wall's loads within 10 % at n = 48"]
fn cut_cell_three_rungs() {
    compare(&[12, 24, 48], 0.1);
}