R7 phase 1: host prototype of the composite Poisson with one nested ratio-2 patch (embedded3::composite; three coarse-fine fluxes, FAC-preconditioned BiCGStab; P1/P1b/P2/P3 gates as ignored tests). Additive module, nothing calls it.
Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
co-authored by
Claude Opus 5.5
Claude Opus 5.5
parent
d63806c0e6
commit
e19cad03b2
+519
@@ -0,0 +1,519 @@
|
||||
//! The composite operator: coarse–coarse, fine–fine and coarse–fine faces,
|
||||
//! the Dirichlet outer boundary, the cut apertures of an optional body on
|
||||
//! the fine level, and the uniform coarse [`Problem`] the preconditioner's
|
||||
//! coarse correction solves.
|
||||
|
||||
use super::Csr;
|
||||
use crate::solvers::incompressible::embedded3::poisson::Problem;
|
||||
use crate::solvers::incompressible::embedded3::{Body, CutGeometry, Grid};
|
||||
|
||||
/// The coarse–fine face flux (see the module docs).
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub enum Interface {
|
||||
Direct,
|
||||
Octree,
|
||||
Quadratic,
|
||||
}
|
||||
|
||||
/// What to build. The coarse grid must be cubic; the patch `lo ≤ (i, j, k)
|
||||
/// < hi` (coarse indices) keeps two coarse cells from the domain boundary.
|
||||
pub struct CompositeSpec<'a> {
|
||||
pub coarse: Grid,
|
||||
pub lo: [usize; 3],
|
||||
pub hi: [usize; 3],
|
||||
pub interface: Interface,
|
||||
/// A body cutting the fine level (φ > 0 fluid), Neumann on its wall.
|
||||
pub body: Option<Sdf>,
|
||||
/// `f = −Δp` (integrated by the midpoint rule over the fluid volume).
|
||||
pub source: &'a (dyn Fn(f64, f64, f64) -> f64 + Sync),
|
||||
/// The Dirichlet value on the outer boundary.
|
||||
pub dirichlet: &'a (dyn Fn(f64, f64, f64) -> f64 + Sync),
|
||||
}
|
||||
|
||||
/// The assembled composite problem.
|
||||
pub struct Composite {
|
||||
pub coarse: Grid,
|
||||
pub fine: Grid,
|
||||
pub lo: [usize; 3],
|
||||
pub hi: [usize; 3],
|
||||
pub interface: Interface,
|
||||
/// Per coarse cell: its unknown, or `usize::MAX` (covered / solid).
|
||||
pub coarse_id: Vec<usize>,
|
||||
/// Per fine cell: its unknown, or `usize::MAX` (solid).
|
||||
pub fine_id: Vec<usize>,
|
||||
/// Unknowns `0..n_coarse` are coarse, the rest fine.
|
||||
pub n_coarse: usize,
|
||||
/// Per unknown: the coarse cell holding it (itself, or the fine
|
||||
/// cell's parent) — the restriction / prolongation map.
|
||||
pub parent: Vec<usize>,
|
||||
/// Per unknown: the cell centre and the fluid volume.
|
||||
pub centre: Vec<[f64; 3]>,
|
||||
pub volume: Vec<f64>,
|
||||
/// Per unknown: the cell owns a coarse–fine face.
|
||||
pub at_interface: Vec<bool>,
|
||||
/// Per unknown: the cell is cut by the body (0 < fluid fraction < 1).
|
||||
pub cut: Vec<bool>,
|
||||
pub a: Csr,
|
||||
pub rhs: Vec<f64>,
|
||||
/// The uniform coarse operator (coarse apertures under the patch too).
|
||||
pub coarse_problem: Problem,
|
||||
}
|
||||
|
||||
const NONE: usize = usize::MAX;
|
||||
/// (axis, sign) of the six faces: −x, +x, −y, +y, −z, +z.
|
||||
const DIRS: [(usize, isize); 6] = [(0, -1), (0, 1), (1, -1), (1, 1), (2, -1), (2, 1)];
|
||||
|
||||
/// A signed distance (φ > 0 fluid) shared by the coarse and fine builds.
|
||||
pub type Sdf = std::sync::Arc<dyn Fn(f64, f64, f64) -> f64 + Send + Sync>;
|
||||
|
||||
/// The body of `phi` on a grid starting at `origin`.
|
||||
fn body_at(phi: &Sdf, origin: [f64; 3]) -> Body {
|
||||
let phi = phi.clone();
|
||||
Body::from_sdf(move |x, y, z, _t| phi(x + origin[0], y + origin[1], z + origin[2]))
|
||||
}
|
||||
|
||||
/// The fluid aperture of face `dir` of cell `c = (i, j, k)` of `cut`.
|
||||
fn aperture(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3], dir: usize) -> f64 {
|
||||
let Some(cut) = cut else { return 1.0 };
|
||||
let [i, j, k] = c;
|
||||
match dir {
|
||||
0 => cut.a_u[g.uface(k, j, i)],
|
||||
1 => cut.a_u[g.uface(k, j, i + 1)],
|
||||
2 => cut.a_v[g.vface(k, j, i)],
|
||||
3 => cut.a_v[g.vface(k, j + 1, i)],
|
||||
4 => cut.a_w[g.wface(k, j, i)],
|
||||
_ => cut.a_w[g.wface(k + 1, j, i)],
|
||||
}
|
||||
}
|
||||
|
||||
fn fraction(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3]) -> f64 {
|
||||
cut.map_or(1.0, |cut| cut.vol[g.cell(c[2], c[1], c[0])])
|
||||
}
|
||||
|
||||
fn active(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3]) -> bool {
|
||||
fraction(cut, g, c) > 0.0 && (0..6).any(|d| aperture(cut, g, c, d) > 0.0)
|
||||
}
|
||||
|
||||
impl Composite {
|
||||
#[must_use]
|
||||
pub fn build(spec: &CompositeSpec<'_>) -> Self {
|
||||
let cg = spec.coarse;
|
||||
let (lo, hi) = (spec.lo, spec.hi);
|
||||
let hc = cg.dx;
|
||||
assert!(
|
||||
(cg.dy - hc).abs() < 1e-14 * hc && (cg.dz - hc).abs() < 1e-14 * hc,
|
||||
"cubic cells"
|
||||
);
|
||||
let dims = [cg.nx, cg.ny, cg.nz];
|
||||
for d in 0..3 {
|
||||
assert!(
|
||||
lo[d] >= 2 && hi[d] + 2 <= dims[d] && hi[d] > lo[d],
|
||||
"patch margin"
|
||||
);
|
||||
}
|
||||
let hf = 0.5 * hc;
|
||||
let fg = Grid::cubic(
|
||||
2 * (hi[0] - lo[0]),
|
||||
2 * (hi[1] - lo[1]),
|
||||
2 * (hi[2] - lo[2]),
|
||||
hf,
|
||||
);
|
||||
let fdims = [fg.nx, fg.ny, fg.nz];
|
||||
let origin = [lo[0] as f64 * hc, lo[1] as f64 * hc, lo[2] as f64 * hc];
|
||||
let fine_cut = spec
|
||||
.body
|
||||
.as_ref()
|
||||
.map(|b| CutGeometry::build(&body_at(b, origin), fg, 0.0));
|
||||
let coarse_cut = spec
|
||||
.body
|
||||
.as_ref()
|
||||
.map(|b| CutGeometry::build(&body_at(b, [0.0; 3]), cg, 0.0));
|
||||
let (fcut, ccut) = (fine_cut.as_ref(), coarse_cut.as_ref());
|
||||
let covered = |c: [usize; 3]| (0..3).all(|d| c[d] >= lo[d] && c[d] < hi[d]);
|
||||
|
||||
// Unknowns.
|
||||
let mut coarse_id = vec![NONE; cg.cells()];
|
||||
let mut fine_id = vec![NONE; fg.cells()];
|
||||
let (mut parent, mut centre, mut volume, mut cut_flag) =
|
||||
(Vec::new(), Vec::new(), Vec::new(), Vec::new());
|
||||
for k in 0..cg.nz {
|
||||
for j in 0..cg.ny {
|
||||
for i in 0..cg.nx {
|
||||
let c = [i, j, k];
|
||||
if covered(c) {
|
||||
continue;
|
||||
}
|
||||
assert!(
|
||||
fraction(ccut, cg, c) > 1.0 - 1e-12,
|
||||
"the body must stay inside the patch"
|
||||
);
|
||||
coarse_id[cg.cell(k, j, i)] = centre.len();
|
||||
parent.push(cg.cell(k, j, i));
|
||||
centre.push([
|
||||
(i as f64 + 0.5) * hc,
|
||||
(j as f64 + 0.5) * hc,
|
||||
(k as f64 + 0.5) * hc,
|
||||
]);
|
||||
volume.push(hc * hc * hc);
|
||||
cut_flag.push(false);
|
||||
}
|
||||
}
|
||||
}
|
||||
let n_coarse = centre.len();
|
||||
for k in 0..fg.nz {
|
||||
for j in 0..fg.ny {
|
||||
for i in 0..fg.nx {
|
||||
let c = [i, j, k];
|
||||
if !active(fcut, fg, c) {
|
||||
continue;
|
||||
}
|
||||
let frac = fraction(fcut, fg, c);
|
||||
fine_id[fg.cell(k, j, i)] = centre.len();
|
||||
parent.push(cg.cell(lo[2] + k / 2, lo[1] + j / 2, lo[0] + i / 2));
|
||||
centre.push([
|
||||
origin[0] + (i as f64 + 0.5) * hf,
|
||||
origin[1] + (j as f64 + 0.5) * hf,
|
||||
origin[2] + (k as f64 + 0.5) * hf,
|
||||
]);
|
||||
volume.push(frac * hf * hf * hf);
|
||||
cut_flag.push(frac < 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
let n = centre.len();
|
||||
let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
|
||||
let mut rhs = vec![0.0; n];
|
||||
let mut at_interface = vec![false; n];
|
||||
|
||||
let step = |c: [usize; 3], d: usize, s: isize| -> [isize; 3] {
|
||||
let mut o = [c[0] as isize, c[1] as isize, c[2] as isize];
|
||||
o[d] += s;
|
||||
o
|
||||
};
|
||||
let cid = |c: [isize; 3]| coarse_id[cg.cell(c[2] as usize, c[1] as usize, c[0] as usize)];
|
||||
let fid = |c: [usize; 3]| fine_id[fg.cell(c[2], c[1], c[0])];
|
||||
// The fine children of covered coarse cell `k` (coarse coords).
|
||||
let children = |kc: [usize; 3]| -> Vec<[usize; 3]> {
|
||||
let b = [
|
||||
2 * (kc[0] - lo[0]),
|
||||
2 * (kc[1] - lo[1]),
|
||||
2 * (kc[2] - lo[2]),
|
||||
];
|
||||
let mut v = Vec::with_capacity(8);
|
||||
for dk in 0..2 {
|
||||
for dj in 0..2 {
|
||||
for di in 0..2 {
|
||||
v.push([b[0] + di, b[1] + dj, b[2] + dk]);
|
||||
}
|
||||
}
|
||||
}
|
||||
v
|
||||
};
|
||||
// The quadratic ghost of fine cell `f` across its face (d, s): the
|
||||
// weights on unknowns.
|
||||
let ghost = |f: [usize; 3], d: usize, s: isize| -> Vec<(usize, f64)> {
|
||||
let mut w = Vec::with_capacity(12);
|
||||
let mut fin = f;
|
||||
fin[d] = (f[d] as isize - s) as usize;
|
||||
w.push((fid(fin), -0.2));
|
||||
w.push((fid(f), 2.0 / 3.0));
|
||||
// The coarse cell holding the ghost point.
|
||||
let gf = [
|
||||
(2 * lo[0] + f[0]) as isize,
|
||||
(2 * lo[1] + f[1]) as isize,
|
||||
(2 * lo[2] + f[2]) as isize,
|
||||
];
|
||||
let mut gc = gf;
|
||||
gc[d] += s;
|
||||
let cc = [
|
||||
gc[0].div_euclid(2),
|
||||
gc[1].div_euclid(2),
|
||||
gc[2].div_euclid(2),
|
||||
];
|
||||
let (t1, t2) = match d {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
let off = |t: usize| if gf[t] % 2 == 0 { -0.25 } else { 0.25 };
|
||||
let (a, b) = (off(t1), off(t2));
|
||||
let wn = 8.0 / 15.0;
|
||||
let at = |dt1: isize, dt2: isize| {
|
||||
let mut c = cc;
|
||||
c[t1] += dt1;
|
||||
c[t2] += dt2;
|
||||
let id = cid(c);
|
||||
assert!(
|
||||
id != NONE,
|
||||
"interpolation stencil on a covered/solid coarse cell"
|
||||
);
|
||||
id
|
||||
};
|
||||
w.push((at(0, 0), wn * (1.0 - a * a - b * b)));
|
||||
w.push((at(-1, 0), wn * (0.5 * a * a - 0.5 * a)));
|
||||
w.push((at(1, 0), wn * (0.5 * a * a + 0.5 * a)));
|
||||
w.push((at(0, -1), wn * (0.5 * b * b - 0.5 * b)));
|
||||
w.push((at(0, 1), wn * (0.5 * b * b + 0.5 * b)));
|
||||
let x = wn * 0.25 * a * b;
|
||||
w.push((at(1, 1), x));
|
||||
w.push((at(-1, -1), x));
|
||||
w.push((at(1, -1), -x));
|
||||
w.push((at(-1, 1), -x));
|
||||
w
|
||||
};
|
||||
|
||||
// Coarse rows.
|
||||
for k in 0..cg.nz {
|
||||
for j in 0..cg.ny {
|
||||
for i in 0..cg.nx {
|
||||
let row = coarse_id[cg.cell(k, j, i)];
|
||||
if row == NONE {
|
||||
continue;
|
||||
}
|
||||
let c = [i, j, k];
|
||||
let x = centre[row];
|
||||
rhs[row] += volume[row] * (spec.source)(x[0], x[1], x[2]);
|
||||
for &(d, s) in &DIRS {
|
||||
let nb = step(c, d, s);
|
||||
if nb[d] < 0 || nb[d] >= dims[d] as isize {
|
||||
let coef = 2.0 * hc;
|
||||
let mut fx = x;
|
||||
fx[d] += s as f64 * 0.5 * hc;
|
||||
rows[row].push((row, coef));
|
||||
rhs[row] += coef * (spec.dirichlet)(fx[0], fx[1], fx[2]);
|
||||
continue;
|
||||
}
|
||||
let nbu = [nb[0] as usize, nb[1] as usize, nb[2] as usize];
|
||||
if !covered(nbu) {
|
||||
let other = cid(nb);
|
||||
if other != NONE {
|
||||
rows[row].push((row, hc));
|
||||
rows[row].push((other, -hc));
|
||||
}
|
||||
continue;
|
||||
}
|
||||
at_interface[row] = true;
|
||||
let kids = children(nbu);
|
||||
// The four children on the face (normal index nearest C).
|
||||
let near: Vec<[usize; 3]> = kids
|
||||
.iter()
|
||||
.copied()
|
||||
.filter(|f| {
|
||||
let want = 2 * (nbu[d] - lo[d]) + usize::from(s < 0);
|
||||
f[d] == want
|
||||
})
|
||||
.collect();
|
||||
match spec.interface {
|
||||
Interface::Direct => {
|
||||
let coef = hf / 1.5;
|
||||
for f in near {
|
||||
rows[row].push((row, coef));
|
||||
rows[row].push((fid(f), -coef));
|
||||
}
|
||||
}
|
||||
Interface::Octree => {
|
||||
rows[row].push((row, hc));
|
||||
for f in kids {
|
||||
rows[row].push((fid(f), -hc / 8.0));
|
||||
}
|
||||
}
|
||||
Interface::Quadratic => {
|
||||
for f in near {
|
||||
// Out of the fine cell through its face −s.
|
||||
rows[row].push((fid(f), -hf));
|
||||
for (u, w) in ghost(f, d, -s) {
|
||||
rows[row].push((u, hf * w));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Fine rows.
|
||||
for k in 0..fg.nz {
|
||||
for j in 0..fg.ny {
|
||||
for i in 0..fg.nx {
|
||||
let row = fine_id[fg.cell(k, j, i)];
|
||||
if row == NONE {
|
||||
continue;
|
||||
}
|
||||
let f = [i, j, k];
|
||||
let x = centre[row];
|
||||
rhs[row] += volume[row] * (spec.source)(x[0], x[1], x[2]);
|
||||
for (dir, &(d, s)) in DIRS.iter().enumerate() {
|
||||
let nb = step(f, d, s);
|
||||
if nb[d] >= 0 && nb[d] < fdims[d] as isize {
|
||||
let nbu = [nb[0] as usize, nb[1] as usize, nb[2] as usize];
|
||||
let a = aperture(fcut, fg, f, dir);
|
||||
let other = fid(nbu);
|
||||
if a > 0.0 && other != NONE {
|
||||
rows[row].push((row, a * hf));
|
||||
rows[row].push((other, -a * hf));
|
||||
}
|
||||
continue;
|
||||
}
|
||||
at_interface[row] = true;
|
||||
assert!(
|
||||
aperture(fcut, fg, f, dir) > 1.0 - 1e-12,
|
||||
"the body must not reach the interface"
|
||||
);
|
||||
let mut gc = [
|
||||
(2 * lo[0] + f[0]) as isize,
|
||||
(2 * lo[1] + f[1]) as isize,
|
||||
(2 * lo[2] + f[2]) as isize,
|
||||
];
|
||||
gc[d] += s;
|
||||
let cc = [
|
||||
gc[0].div_euclid(2),
|
||||
gc[1].div_euclid(2),
|
||||
gc[2].div_euclid(2),
|
||||
];
|
||||
let c_row = cid(cc);
|
||||
match spec.interface {
|
||||
Interface::Direct => {
|
||||
let coef = hf / 1.5;
|
||||
rows[row].push((row, coef));
|
||||
rows[row].push((c_row, -coef));
|
||||
}
|
||||
Interface::Octree => {
|
||||
let coef = hf * hf / hc;
|
||||
let kc = [lo[0] + i / 2, lo[1] + j / 2, lo[2] + k / 2];
|
||||
for sib in children(kc) {
|
||||
rows[row].push((fid(sib), coef / 8.0));
|
||||
}
|
||||
rows[row].push((c_row, -coef));
|
||||
}
|
||||
Interface::Quadratic => {
|
||||
rows[row].push((row, hf));
|
||||
for (u, w) in ghost(f, d, s) {
|
||||
rows[row].push((u, -hf * w));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let a = Csr::from_rows(rows);
|
||||
let coarse_problem = Self::coarse_operator(cg, ccut);
|
||||
Self {
|
||||
coarse: cg,
|
||||
fine: fg,
|
||||
lo,
|
||||
hi,
|
||||
interface: spec.interface,
|
||||
coarse_id,
|
||||
fine_id,
|
||||
n_coarse,
|
||||
parent,
|
||||
centre,
|
||||
volume,
|
||||
at_interface,
|
||||
cut: cut_flag,
|
||||
a,
|
||||
rhs,
|
||||
coarse_problem,
|
||||
}
|
||||
}
|
||||
|
||||
/// The uniform coarse operator of the same problem (the coarse cut
|
||||
/// under the patch), Dirichlet outside: the FAC coarse correction.
|
||||
fn coarse_operator(cg: Grid, ccut: Option<&CutGeometry>) -> Problem {
|
||||
let h = cg.dx;
|
||||
let mut p = Problem::new(cg.nx, cg.ny, cg.nz);
|
||||
let dims = [cg.nx, cg.ny, cg.nz];
|
||||
for k in 0..cg.nz {
|
||||
for j in 0..cg.ny {
|
||||
for i in 0..cg.nx {
|
||||
let idx = cg.cell(k, j, i);
|
||||
p.active[idx] = active(ccut, cg, [i, j, k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
for k in 0..cg.nz {
|
||||
for j in 0..cg.ny {
|
||||
for i in 0..cg.nx {
|
||||
let idx = cg.cell(k, j, i);
|
||||
if !p.active[idx] {
|
||||
continue;
|
||||
}
|
||||
let c = [i, j, k];
|
||||
for (dir, &(d, s)) in DIRS.iter().enumerate() {
|
||||
let mut nb = [i as isize, j as isize, k as isize];
|
||||
nb[d] += s;
|
||||
if nb[d] < 0 || nb[d] >= dims[d] as isize {
|
||||
p.extra_diag[idx] += 2.0 * h;
|
||||
continue;
|
||||
}
|
||||
let other = cg.cell(nb[2] as usize, nb[1] as usize, nb[0] as usize);
|
||||
let coef = if p.active[other] {
|
||||
aperture(ccut, cg, c, dir) * h
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
match dir {
|
||||
0 => p.aw[idx] = coef,
|
||||
1 => p.ae[idx] = coef,
|
||||
2 => p.as_[idx] = coef,
|
||||
3 => p.an[idx] = coef,
|
||||
4 => p.ab[idx] = coef,
|
||||
_ => p.at[idx] = coef,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
p
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn unknowns(&self) -> usize {
|
||||
self.centre.len()
|
||||
}
|
||||
}
|
||||
|
||||
/// A uniform grid's operator and right-hand side in the same form (the
|
||||
/// reference / cost comparator), optionally cut by `body`; Dirichlet
|
||||
/// outside. Returns the problem and the per-cell fluid fraction.
|
||||
#[must_use]
|
||||
pub fn uniform_problem(
|
||||
g: Grid,
|
||||
body: Option<&Sdf>,
|
||||
source: &(dyn Fn(f64, f64, f64) -> f64 + Sync),
|
||||
dirichlet: &(dyn Fn(f64, f64, f64) -> f64 + Sync),
|
||||
) -> (Problem, Vec<f64>) {
|
||||
let cut = body.map(|b| CutGeometry::build(&body_at(b, [0.0; 3]), g, 0.0));
|
||||
let mut p = Composite::coarse_operator(g, cut.as_ref());
|
||||
let h = g.dx;
|
||||
let dims = [g.nx, g.ny, g.nz];
|
||||
let mut frac = vec![0.0; g.cells()];
|
||||
for k in 0..g.nz {
|
||||
for j in 0..g.ny {
|
||||
for i in 0..g.nx {
|
||||
let idx = g.cell(k, j, i);
|
||||
if !p.active[idx] {
|
||||
continue;
|
||||
}
|
||||
let fr = fraction(cut.as_ref(), g, [i, j, k]);
|
||||
frac[idx] = fr;
|
||||
let x = [
|
||||
(i as f64 + 0.5) * h,
|
||||
(j as f64 + 0.5) * h,
|
||||
(k as f64 + 0.5) * h,
|
||||
];
|
||||
p.rhs[idx] = fr * h * h * h * source(x[0], x[1], x[2]);
|
||||
for &(d, s) in &DIRS {
|
||||
let c = [i, j, k][d] as isize + s;
|
||||
if c < 0 || c >= dims[d] as isize {
|
||||
let mut fx = x;
|
||||
fx[d] += s as f64 * 0.5 * h;
|
||||
p.rhs[idx] += 2.0 * h * dirichlet(fx[0], fx[1], fx[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
(p, frac)
|
||||
}
|
||||
@@ -0,0 +1,103 @@
|
||||
//! R7 phase 1 (omni-cortex roadmap R7, A3-ii): the composite pressure
|
||||
//! Poisson problem on a uniform coarse grid with ONE nested ratio-2 box
|
||||
//! patch, cell-centred, in the integrated (finite-volume) form of
|
||||
//! [`super::poisson::Problem`]: `Σ_faces c_f (p_i − p_nb) = rhs_i`.
|
||||
//!
|
||||
//! A host prototype, not wired into the solver (nothing in the step calls
|
||||
//! it): it exists to measure what the coarse–fine interface costs in order
|
||||
//! and in work before the patch is carried into the predictor and onto the
|
||||
//! device. Three interface fluxes are built so the choice is a measurement:
|
||||
//!
|
||||
//! - [`Interface::Direct`]: each fine sub-face couples its fine cell to the
|
||||
//! coarse cell by a two-point flux over the centre distance `1.5 h_f`
|
||||
//! (symmetric, conservative, not consistent tangentially);
|
||||
//! - [`Interface::Octree`]: 2604.18886 eq. 13 — the coarse face's flux
|
||||
//! `(p_C − mean of the covered coarse cell's 8 children) / h_c`, shared
|
||||
//! equally by the 4 fine sub-faces (conservative, non-symmetric);
|
||||
//! - [`Interface::Quadratic`]: the fine ghost by quadratic interpolation
|
||||
//! (tangential quadratic with the cross term on the coarse layer, then
|
||||
//! quadratic in the normal direction through two fine cells), the coarse
|
||||
//! face's flux = the sum of its four fine fluxes (refluxing; conservative,
|
||||
//! non-symmetric) — the AMReX / Martin–Cartwright construction.
|
||||
//!
|
||||
//! Layout: coarse cells `(k, j, i)` of [`Grid`]; the patch covers the
|
||||
//! coarse box `lo ≤ (i, j, k) < hi` and is refined by 2 in every direction.
|
||||
//! Unknowns: the uncovered active coarse cells first, then the active fine
|
||||
//! cells. A body (φ > 0 fluid) may cut the fine level; it must not reach
|
||||
//! the interface or the coarse uncovered cells.
|
||||
|
||||
mod assemble;
|
||||
mod solve;
|
||||
|
||||
pub use assemble::{Composite, CompositeSpec, Interface, Sdf, uniform_problem};
|
||||
pub use solve::{CompositeSolve, SolveStats, solve_bicgstab};
|
||||
|
||||
/// Compressed sparse rows.
|
||||
#[derive(Debug, Clone, Default)]
|
||||
pub struct Csr {
|
||||
pub start: Vec<usize>,
|
||||
pub col: Vec<usize>,
|
||||
pub val: Vec<f64>,
|
||||
}
|
||||
|
||||
impl Csr {
|
||||
/// Rows from per-row `(column, value)` lists; duplicates are summed in
|
||||
/// the order they appear, zero entries dropped.
|
||||
#[must_use]
|
||||
pub fn from_rows(rows: Vec<Vec<(usize, f64)>>) -> Self {
|
||||
let mut start = Vec::with_capacity(rows.len() + 1);
|
||||
let mut col = Vec::new();
|
||||
let mut val = Vec::new();
|
||||
start.push(0);
|
||||
for mut row in rows {
|
||||
row.sort_by_key(|e| e.0);
|
||||
let mut last: Option<usize> = None;
|
||||
for (c, v) in row {
|
||||
if last == Some(c) {
|
||||
*val.last_mut().expect("entry") += v;
|
||||
} else {
|
||||
col.push(c);
|
||||
val.push(v);
|
||||
last = Some(c);
|
||||
}
|
||||
}
|
||||
start.push(col.len());
|
||||
}
|
||||
Self { start, col, val }
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn rows(&self) -> usize {
|
||||
self.start.len() - 1
|
||||
}
|
||||
|
||||
pub fn apply(&self, x: &[f64], y: &mut [f64]) {
|
||||
use rayon::prelude::*;
|
||||
y.par_iter_mut().enumerate().for_each(|(r, yr)| {
|
||||
let mut s = 0.0;
|
||||
for e in self.start[r]..self.start[r + 1] {
|
||||
s += self.val[e] * x[self.col[e]];
|
||||
}
|
||||
*yr = s;
|
||||
});
|
||||
}
|
||||
|
||||
/// The largest `|a_ij − a_ji|` relative to the largest |a_ij|.
|
||||
#[must_use]
|
||||
pub fn asymmetry(&self) -> f64 {
|
||||
let mut map = std::collections::HashMap::with_capacity(self.val.len());
|
||||
let mut amax: f64 = 0.0;
|
||||
for r in 0..self.rows() {
|
||||
for e in self.start[r]..self.start[r + 1] {
|
||||
map.insert((r, self.col[e]), self.val[e]);
|
||||
amax = amax.max(self.val[e].abs());
|
||||
}
|
||||
}
|
||||
let mut d: f64 = 0.0;
|
||||
for (&(r, c), &v) in &map {
|
||||
let t = map.get(&(c, r)).copied().unwrap_or(0.0);
|
||||
d = d.max((v - t).abs());
|
||||
}
|
||||
d / amax.max(f64::MIN_POSITIVE)
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,195 @@
|
||||
//! The composite solve: BiCGStab (the Quadratic and Octree operators are
|
||||
//! not symmetric) right-preconditioned by one FAC cycle (McCormick's fast
|
||||
//! adaptive composite grid method in correction form):
|
||||
//!
|
||||
//! 1. `ν` forward Gauss–Seidel sweeps of the composite operator from zero;
|
||||
//! 2. the composite residual, restricted to the whole coarse grid — an
|
||||
//! uncovered coarse cell keeps its own, a covered coarse cell receives
|
||||
//! the SUM of its children's (integrated form);
|
||||
//! 3. one V-cycle of the uniform coarse operator (embedded3's aggregation
|
||||
//! [`Hierarchy`], unchanged) on that residual;
|
||||
//! 4. the correction added to the coarse unknowns and injected (piecewise
|
||||
//! constant) into the fine ones;
|
||||
//! 5. `ν` backward Gauss–Seidel sweeps.
|
||||
//!
|
||||
//! Correction form over the COMPOSITE residual is the linear counterpart of
|
||||
//! 2604.18886's FAS right-hand side `b = β R r + A_c u*` (Algorithm 4): the
|
||||
//! interface flux enters the coarse problem exactly once, through the
|
||||
//! composite residual, so nothing is double-counted at the T-junctions.
|
||||
|
||||
use super::Composite;
|
||||
use crate::solvers::incompressible::MultigridParameters;
|
||||
use crate::solvers::incompressible::embedded3::poisson::Hierarchy;
|
||||
|
||||
/// The prepared preconditioner.
|
||||
pub struct CompositeSolve {
|
||||
hier: Hierarchy<f64>,
|
||||
diag: Vec<f64>,
|
||||
sweeps: usize,
|
||||
rc: Vec<f64>,
|
||||
ec: Vec<f64>,
|
||||
res: Vec<f64>,
|
||||
}
|
||||
|
||||
/// Outcome of [`solve_bicgstab`].
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub struct SolveStats {
|
||||
pub iterations: usize,
|
||||
/// `‖b − A x‖₂ / ‖b‖₂` at exit (recomputed).
|
||||
pub rel_residual: f64,
|
||||
pub converged: bool,
|
||||
pub setup_s: f64,
|
||||
pub solve_s: f64,
|
||||
}
|
||||
|
||||
impl CompositeSolve {
|
||||
#[must_use]
|
||||
pub fn new(c: &Composite, sweeps: usize) -> Self {
|
||||
let params = MultigridParameters::default();
|
||||
let hier = Hierarchy::<f64>::build(&c.coarse_problem, ¶ms);
|
||||
let n = c.unknowns();
|
||||
let mut diag = vec![0.0; n];
|
||||
for (r, d) in diag.iter_mut().enumerate() {
|
||||
for e in c.a.start[r]..c.a.start[r + 1] {
|
||||
if c.a.col[e] == r {
|
||||
*d = c.a.val[e];
|
||||
}
|
||||
}
|
||||
assert!(*d > 0.0, "row {r} has no positive diagonal");
|
||||
}
|
||||
let nc = c.coarse.cells();
|
||||
Self {
|
||||
hier,
|
||||
diag,
|
||||
sweeps: sweeps.max(1),
|
||||
rc: vec![0.0; nc],
|
||||
ec: vec![0.0; nc],
|
||||
res: vec![0.0; n],
|
||||
}
|
||||
}
|
||||
|
||||
fn gs_row(&self, c: &Composite, b: &[f64], x: &mut [f64], r: usize) {
|
||||
let mut s = b[r];
|
||||
for e in c.a.start[r]..c.a.start[r + 1] {
|
||||
let col = c.a.col[e];
|
||||
if col != r {
|
||||
s -= c.a.val[e] * x[col];
|
||||
}
|
||||
}
|
||||
x[r] = s / self.diag[r];
|
||||
}
|
||||
|
||||
/// `z = M⁻¹ r` (one FAC cycle from zero).
|
||||
pub fn precondition(&mut self, c: &Composite, r: &[f64], z: &mut [f64]) {
|
||||
let n = c.unknowns();
|
||||
z.iter_mut().for_each(|v| *v = 0.0);
|
||||
for _ in 0..self.sweeps {
|
||||
for row in 0..n {
|
||||
self.gs_row(c, r, z, row);
|
||||
}
|
||||
}
|
||||
c.a.apply(z, &mut self.res);
|
||||
self.rc.iter_mut().for_each(|v| *v = 0.0);
|
||||
for u in 0..n {
|
||||
self.rc[c.parent[u]] += r[u] - self.res[u];
|
||||
}
|
||||
for (idx, v) in self.rc.iter_mut().enumerate() {
|
||||
if !c.coarse_problem.active[idx] {
|
||||
*v = 0.0;
|
||||
}
|
||||
}
|
||||
self.hier.apply_preconditioner(&self.rc, &mut self.ec);
|
||||
for u in 0..n {
|
||||
let p = c.parent[u];
|
||||
if c.coarse_problem.active[p] {
|
||||
z[u] += self.ec[p];
|
||||
}
|
||||
}
|
||||
for _ in 0..self.sweeps {
|
||||
for row in (0..n).rev() {
|
||||
self.gs_row(c, r, z, row);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn dot(a: &[f64], b: &[f64]) -> f64 {
|
||||
a.iter().zip(b).map(|(x, y)| x * y).sum()
|
||||
}
|
||||
|
||||
/// Solve `A x = rhs` to `‖r‖₂ ≤ tol ‖rhs‖₂` from the given `x`.
|
||||
pub fn solve_bicgstab(
|
||||
c: &Composite,
|
||||
x: &mut [f64],
|
||||
tol: f64,
|
||||
max_it: usize,
|
||||
sweeps: usize,
|
||||
) -> SolveStats {
|
||||
let t0 = std::time::Instant::now();
|
||||
let mut pre = CompositeSolve::new(c, sweeps);
|
||||
let setup_s = t0.elapsed().as_secs_f64();
|
||||
let t1 = std::time::Instant::now();
|
||||
let n = c.unknowns();
|
||||
let b = &c.rhs;
|
||||
let bn = dot(b, b).sqrt().max(f64::MIN_POSITIVE);
|
||||
let mut r = vec![0.0; n];
|
||||
c.a.apply(x, &mut r);
|
||||
for i in 0..n {
|
||||
r[i] = b[i] - r[i];
|
||||
}
|
||||
let rhat = r.clone();
|
||||
let (mut rho, mut alpha, mut omega) = (1.0, 1.0, 1.0);
|
||||
let mut v = vec![0.0; n];
|
||||
let mut p = vec![0.0; n];
|
||||
let mut ph = vec![0.0; n];
|
||||
let mut sh = vec![0.0; n];
|
||||
let mut s = vec![0.0; n];
|
||||
let mut t = vec![0.0; n];
|
||||
let mut it = 0;
|
||||
let mut rel = dot(&r, &r).sqrt() / bn;
|
||||
while rel > tol && it < max_it {
|
||||
it += 1;
|
||||
let rho_n = dot(&rhat, &r);
|
||||
let beta = (rho_n / rho) * (alpha / omega);
|
||||
for i in 0..n {
|
||||
p[i] = r[i] + beta * (p[i] - omega * v[i]);
|
||||
}
|
||||
pre.precondition(c, &p, &mut ph);
|
||||
c.a.apply(&ph, &mut v);
|
||||
alpha = rho_n / dot(&rhat, &v);
|
||||
for i in 0..n {
|
||||
s[i] = r[i] - alpha * v[i];
|
||||
}
|
||||
if dot(&s, &s).sqrt() / bn <= tol {
|
||||
for i in 0..n {
|
||||
x[i] += alpha * ph[i];
|
||||
}
|
||||
r.copy_from_slice(&s);
|
||||
break;
|
||||
}
|
||||
pre.precondition(c, &s, &mut sh);
|
||||
c.a.apply(&sh, &mut t);
|
||||
omega = dot(&t, &s) / dot(&t, &t);
|
||||
for i in 0..n {
|
||||
x[i] += alpha * ph[i] + omega * sh[i];
|
||||
r[i] = s[i] - omega * t[i];
|
||||
}
|
||||
rho = rho_n;
|
||||
rel = dot(&r, &r).sqrt() / bn;
|
||||
}
|
||||
// The true residual.
|
||||
c.a.apply(x, &mut r);
|
||||
let mut rr = 0.0;
|
||||
for i in 0..n {
|
||||
let d = b[i] - r[i];
|
||||
rr += d * d;
|
||||
}
|
||||
let rel = rr.sqrt() / bn;
|
||||
SolveStats {
|
||||
iterations: it,
|
||||
rel_residual: rel,
|
||||
converged: rel <= 10.0 * tol,
|
||||
setup_s,
|
||||
solve_s: t1.elapsed().as_secs_f64(),
|
||||
}
|
||||
}
|
||||
@@ -8,6 +8,7 @@
|
||||
|
||||
pub mod body;
|
||||
pub mod closure;
|
||||
pub mod composite;
|
||||
pub mod cut;
|
||||
pub mod cutwall;
|
||||
pub mod exchange;
|
||||
|
||||
Reference in New Issue
Block a user