rustytorch/crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/balance.rs

//! P4 option B, patch side (`docs/overset_metal_campaign.md` §5.11): the
//! patch's own momentum balance on its SOLVED cells (everything but the
//! acceptor row), in the scheme's own fluxes.
//!
//! The marched equation on a solved cell is the face-flux form of
//! `predictor.rs` (outward `ρ F_f u_f` with the scheme's face value,
//! outward `μ L_f(u)` in Laplacian form) plus the projection's pressure
//! force, which is the least-squares CELL gradient `−A_c ∇p_c` — not a
//! face pressure, so the pressure term does not telescope: summed over
//! the solved cells it need not equal the boundary integral of any face
//! pressure, and that difference `δP` is the patch's momentum
//! non-conservation. Everything here is a plain sum of the scheme's
//! terms; the steady balance `−conv_acc + visc_acc + visc_wall − p_ls = 0`
//! is the gate (the unsteady term is not stored by `PatchField`; at the
//! settled state it is small and reported as omitted).

use super::{CurvilinearPisoSolver, PatchConvection, PatchField, SideBc};
use crate::mesh::{PatchMesh, PatchSide};

/// The pieces of the patch's momentum balance, x / y, N/m.
#[derive(Debug, Clone, Copy, Default)]
pub struct PatchBalance {
    /// Outward `ρ F_f u_f` through the solved/acceptor interface faces.
    pub conv_acc: [f64; 2],
    /// Outward `μ ∇u · S` (Laplacian form) through the interface faces.
    pub visc_acc: [f64; 2],
    /// `Σ p_f S_out` on the interface faces, `p_f` linear in the two cells.
    pub p_face_acc: [f64; 2],
    /// Outward `μ ∇u · S` through the wall faces (Dirichlet wall value).
    pub visc_wall: [f64; 2],
    /// `Σ p_f S_out` on the wall faces, `p_f` extrapolated as
    /// `surface_force` does (cell value + least-squares gradient).
    pub p_face_wall: [f64; 2],
    /// `Σ_solved A_c ∇p_c` — the pressure force the scheme applied.
    pub p_ls: [f64; 2],
    /// Solved cells.
    pub cells: usize,
    /// Interface faces.
    pub acc_faces: usize,
    /// Wall faces.
    pub wall_faces: usize,
}

impl PatchBalance {
    /// Steady balance residual of the marched equation on the solved
    /// region: `−conv_acc + visc_acc + visc_wall − p_ls` (the gate).
    pub fn balance(&self) -> [f64; 2] {
        [0, 1].map(|k| -self.conv_acc[k] + self.visc_acc[k] + self.visc_wall[k] - self.p_ls[k])
    }
    /// The body force read through the interface in flux form:
    /// `∮ (σ·n − ρ u u·n)` with `σ` in the scheme's Laplacian form.
    pub fn flux_force(&self) -> [f64; 2] {
        [0, 1].map(|k| -self.p_face_acc[k] + self.visc_acc[k] - self.conv_acc[k])
    }
    /// The wall force in the scheme's own wall fluxes (Laplacian form).
    pub fn wall_force(&self) -> [f64; 2] {
        [0, 1].map(|k| self.p_face_wall[k] - self.visc_wall[k])
    }
    /// The pressure non-conservation `p_ls − p_face_acc − p_face_wall`
    /// (= `flux_force − wall_force` when the balance holds).
    pub fn pressure_defect(&self) -> [f64; 2] {
        [0, 1].map(|k| self.p_ls[k] - self.p_face_acc[k] - self.p_face_wall[k])
    }
}

impl CurvilinearPisoSolver {
    /// The momentum balance of the solved cells at time `t` (see the
    /// module doc). Stationary mesh only.
    pub fn momentum_balance(&self, field: &PatchField, t: f64) -> PatchBalance {
        let mesh: &PatchMesh = &self.mesh;
        let rho = self.config.density;
        let mu = self.config.viscosity;
        let ops = &self.ops;
        let bvel = |side: PatchSide, xy: [f64; 2]| -> Option<(f64, f64)> {
            match self.params.boundaries.get(side) {
                SideBc::Velocity => Some(self.boundary_velocity(side, xy[0], xy[1], t)),
                SideBc::Outlet => None,
            }
        };
        let un = ops.node_values(mesh, &field.u, &|s, xy| bvel(s, xy).map(|v| v.0));
        let vn = ops.node_values(mesh, &field.v, &|s, xy| bvel(s, xy).map(|v| v.1));
        let mut b = PatchBalance::default();
        for c in 0..mesh.cell_count() {
            if self.is_acceptor(c) {
                continue;
            }
            b.cells += 1;
            let g = self.pressure_gradient(&field.p, c);
            let a = mesh.area(c);
            b.p_ls[0] += a * g[0];
            b.p_ls[1] += a * g[1];
            for (f, sign) in mesh.cell_faces(c) {
                let face = &mesh.faces()[f];
                let s_out = [sign * face.s[0], sign * face.s[1]];
                match (face.owner, face.neigh) {
                    (Some(p), Some(q)) => {
                        let other = if p == c { q } else { p };
                        if !self.is_acceptor(other) {
                            continue;
                        }
                        b.acc_faces += 1;
                        let out = sign * field.flux[f];
                        let (uf, vf) = match self.params.convection {
                            PatchConvection::Upwind => {
                                let up = if out >= 0.0 { c } else { other };
                                (field.u[up], field.v[up])
                            }
                            PatchConvection::TvdVanAlbada => {
                                let (up, dn) = if out >= 0.0 { (c, other) } else { (other, c) };
                                let (du, dv) = self.tvd_correction(field, f, up, dn);
                                (field.u[up] + du, field.v[up] + dv)
                            }
                            PatchConvection::None => (0.0, 0.0),
                        };
                        b.conv_acc[0] += rho * out * uf;
                        b.conv_acc[1] += rho * out * vf;
                        b.visc_acc[0] +=
                            mu * sign * ops.face_gradient_flux(mesh, f, &field.u, &un, None);
                        b.visc_acc[1] +=
                            mu * sign * ops.face_gradient_flux(mesh, f, &field.v, &vn, None);
                        let pf = face.w * field.p[p] + (1.0 - face.w) * field.p[q];
                        b.p_face_acc[0] += pf * s_out[0];
                        b.p_face_acc[1] += pf * s_out[1];
                    }
                    _ => {
                        let side = mesh.side(f).expect("boundary face has a side");
                        if side != PatchSide::Inner {
                            continue;
                        }
                        b.wall_faces += 1;
                        let bv = bvel(side, face.centre);
                        b.visc_wall[0] += mu
                            * sign
                            * ops.face_gradient_flux(mesh, f, &field.u, &un, bv.map(|v| v.0));
                        b.visc_wall[1] += mu
                            * sign
                            * ops.face_gradient_flux(mesh, f, &field.v, &vn, bv.map(|v| v.1));
                        let xc = mesh.centre(c);
                        let dxf = [face.centre[0] - xc[0], face.centre[1] - xc[1]];
                        let pf = field.p[c] + g[0] * dxf[0] + g[1] * dxf[1];
                        b.p_face_wall[0] += pf * s_out[0];
                        b.p_face_wall[1] += pf * s_out[1];
                    }
                }
            }
        }
        b
    }
}