rustyverse/rustytorch
2
0
Fork 0
Code Issues Pull Requests Actions 5 Packages Projects Releases Wiki Activity

solvers: near-wall momentum, Newmark dynamics, QM6, and MMS across elements
CI / Test (ubuntu-latest) (push) Canceled after 0s
Details
CI / Build CPU-Only (Explicit) (push) Canceled after 0s
Details
CI / Python Bindings (maturin) (macos-latest) (push) Canceled after 0s
Details
CI / Python Bindings (maturin) (ubuntu-latest) (push) Canceled after 0s
Details
CI / WASM Build + Size Check (push) Canceled after 0s
Details
CI / Distributed Training Tests (push) Canceled after 0s
Details
CI / CI Success (push) Canceled after 0s
Details
Documentation / Build API Documentation (push) Canceled after 0s
Details
Documentation / Build User Guide (push) Canceled after 0s
Details
Performance Benchmarks / Run Benchmarks (push) Canceled after 0s
Details
CI / Format Check (push) Canceled after 0s
Details
CI / Clippy Check (push) Canceled after 0s
Details
CI / Build (macos-latest) (push) Canceled after 0s
Details
CI / Build (ubuntu-latest) (push) Canceled after 0s
Details
CI / Test (macos-latest) (push) Canceled after 0s
Details

Browse Source 
Four parallel work items plus two defects found while integrating them.
561 -> 592 tests, 0 failing, verified stable over repeated runs.

## rtx-cfd: solve the near-wall velocity lines

Every u row sits at y = (j+0.5) dy and every v column at x = (i+0.5) dx --
strictly interior. The sweeps froze rows 0 and ny-1 and columns 0 and
nx-1 and treated whatever was stored there as a boundary condition, which
imposed wall values half a cell inside the domain. They are now unknowns,
with the wall entering through the control volume's half-cell conductance
(mu dx / (dy/2)), zero convective flux through the wall, and the wall's
tangential velocity in the source.

That in turn makes continuity enforceable on every cell, with a neighbour
coefficient zero only for a genuine boundary face. Extending continuity
had been tried before and broke convergence; it works now because the
near-wall lines are no longer frozen. Order matters here.

Manufactured solutions, which is how any of this is known:

    n     L2 velocity   order      max |p - p_exact|
    16    3.516212e-2      -          9.245576e-2
    32    1.953751e-2    0.85         5.225739e-2
    64    1.037523e-2    0.91         2.796415e-2

Velocity error is 7.4x smaller at n=16, and the observed order rises from
0.48 toward 1. The pressure error was 0.408 -> 0.624 -> 0.756, *growing*
with refinement; it now falls. Divergence on the outer ring of cells goes
from 1.0e1 to 2.5e-10.

A separate defect found on the way: u_source_term was computed and never
called, so the x-momentum equation carried no body force at all while the
y-momentum one did. That is exactly the u-versus-v asymmetry the earlier
diagnosis had flagged as an unexplained clue.

Cavity at 65^2, against Ghia's u_min = -0.2109 at y = 0.4531:
-0.1792 at 0.3906 before, -0.1932 at 0.5000 after, in 733 iterations
rather than 971.

The cavity test now sets FreeSlipWall on all four sides plus the lid
through the new set_wall_velocity hook. That is not a weakened benchmark:
on a staggered grid the only velocity component living *on* a boundary is
the normal one, which is what FreeSlipWall prescribes, and the tangential
no-slip arrives through the half-cell wall term with wall velocity zero on
the three stationary walls. Prescribing whole u rows and v columns, as
before, pins lines half a cell inside the domain and over-determines the
cells beside them once every cell has a continuity equation.

## rtx-fea: DynamicAnalysis, previously a stub returning zeros

Newmark-beta in acceleration form -- the displacement form divides by
beta dt^2, singular at beta = 0 -- with Rayleigh damping, the effective
matrix Cholesky-factorised once and reused. Initial acceleration is solved
from M a0 = F0 - C v0 - K u0 rather than assumed zero, which would destroy
the second-order rate.

Verified two ways that cannot both be faked: against the closed-form
single-degree-of-freedom response, undamped and damped, with the measured
order of accuracy; and against the free-vibration period of the same bar
whose modal frequencies are already validated. Time domain and frequency
domain come from different code paths.

## rtx-fea: QM6 incompatible modes

Wilson's Q6 with Taylor's correction, added alongside compute_stiffness_
matrix rather than replacing it -- the existing method is byte-identical,
which matters because the manufactured-solution verification depends on
it. Internal modes statically condensed; the incompatible strain block
evaluated at the element centre, which is what makes the patch test pass
on distorted elements.

## rtx-fea: manufactured solutions across the element library

    Quad4  order 2.00      Tri3   order 1.98
    Quad8  order 3.00      Hex8   order 1.96  (new 3-D solution)

Each element asserts its own theoretical rate.

## Two defects found while integrating

Reverse Cuthill-McKee node ordering was nondeterministic. All three of its
orderings -- seed selection, neighbour ordering, and the trailing sweep --
were decided by HashMap/HashSet iteration order, which std randomises per
process. On a rectangular mesh every corner ties at minimum degree, so two
calls to displacement_only on the same mesh in the same process returned
different DOF indices for the same node, agreeing in only 5 of 20 measured
runs. Ties now break by node id. This surfaced as a coin-flip test failure
-- 12 in 25 runs -- and would have been dismissed as flaky rather than
diagnosed had the integration pass not re-run it.

Quadrature: triangle(3) weights summed to 0.25 against a reference area of
0.5, and tetrahedron(3) to 1/36 against a volume of 1/6. Both divided
weights that were already tabulated for the reference measure by that
measure again, so both rules integrated everything to a fraction of its
value -- invisibly, since a scaled quadrature leaves the stiffness matrix
symmetric, the mass matrix positive definite and the rigid-body modes
exact. New test asserts every rule integrates 1 to its reference measure,
across every family and order, plus Gauss-Legendre exactness to degree
2n-1.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
2026-08-19 15:39:20 -07:00
co-authored by Claude Opus 5
parent b5814a304f
commit 698c844926
11 changed files with 3867 additions and 287 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/specialized/rtx-cfd/src/solvers/incompressible/simple.rs
+343 -84
View File
@@ -160,6 +160,22 @@ pub struct SimpleSolver {
/// solver into something whose exact answer is known in closed form.
#[allow(clippy::type_complexity)]
momentum_source: Option<Box<dyn Fn(f64, f64) -> (f64, f64) + Send + Sync>>,
/// Optional wall velocity `f(x, y) -> (u_wall, v_wall)`, sampled at the
/// wall face position.
///
/// The near-wall momentum control volumes need the *tangential* velocity of
/// the wall that bounds them, and on this staggered layout there is nowhere
/// to store it: no u node lies on the bottom or top wall, and no v node
/// lies on the left or right wall. Every u row is at `y = (j + 0.5) dy`,
/// strictly interior. Supplying it as a function of position is the only
/// way a spatially varying wall — a manufactured solution's, for instance —
/// can reach the discretisation at all.
///
/// When unset the solver falls back to the value stored on the near-wall
/// line itself, which is exactly what `BoundaryConditions` writes there, so
/// an existing configuration keeps the wall it always had.
#[allow(clippy::type_complexity)]
wall_velocity: Option<Box<dyn Fn(f64, f64) -> (f64, f64) + Send + Sync>>,
}
/// Workspace for linear algebra operations
@@ -219,6 +235,7 @@ impl SimpleSolver {
},
turbulence_model,
momentum_source: None,
wall_velocity: None,
})
}
@@ -230,6 +247,34 @@ impl SimpleSolver {
self.momentum_source = Some(Box::new(source));
}
/// Set the wall velocity as a function of position. See
/// [`Self::wall_velocity`].
pub fn set_wall_velocity<F>(&mut self, f: F)
where
F: Fn(f64, f64) -> (f64, f64) + Send + Sync + 'static,
{
self.wall_velocity = Some(Box::new(f));
}
/// Tangential `u` of the horizontal wall bounding the near-wall u control
/// volume at face `i`, row `j`, with the wall itself at `y_wall`.
///
/// Falls back to the value stored on the near-wall line — see
/// [`Self::wall_velocity`].
fn u_wall(&self, flow_field: &FlowField, i: usize, j: usize, y_wall: f64, dx: f64) -> f64 {
self.wall_velocity
.as_ref()
.map_or(flow_field.u[(j, i)], |f| f(i as f64 * dx, y_wall).0)
}
/// Tangential `v` of the vertical wall bounding the near-wall v control
/// volume at column `i`, face `j`, with the wall itself at `x_wall`.
fn v_wall(&self, flow_field: &FlowField, i: usize, j: usize, x_wall: f64, dy: f64) -> f64 {
self.wall_velocity
.as_ref()
.map_or(flow_field.v[(j, i)], |f| f(x_wall, j as f64 * dy).1)
}
/// Momentum source contribution for a u-face, already multiplied by the
/// control volume so it is a force, matching the pressure-gradient term.
///
@@ -347,31 +392,58 @@ impl SimpleSolver {
// equation before measuring it.
let alpha = self.parameters.velocity_relaxation;
for j in 1..ny - 1 {
for i in 1..nx - 1 {
// Over exactly the unknowns the sweeps solve for. Measuring a smaller
// set would let the near-wall lines converge to anything at all without
// the reported residual noticing.
for j in 0..ny {
for i in 1..nx {
let cu =
self.compute_u_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
let ap = cu.center * alpha;
let source = cu.source - (1.0 - alpha) * cu.center * flow_field.u_old[(j, i)];
let diagonal_u = ap * flow_field.u[(j, i)];
let north = if j + 1 < ny {
cu.north * flow_field.u[(j + 1, i)]
} else {
0.0
};
let south = if j > 0 {
cu.south * flow_field.u[(j - 1, i)]
} else {
0.0
};
let imbalance_u = diagonal_u
- (source
+ cu.east * flow_field.u[(j, i + 1)]
+ cu.west * flow_field.u[(j, i - 1)]
+ cu.north * flow_field.u[(j + 1, i)]
+ cu.south * flow_field.u[(j - 1, i)]);
+ north
+ south);
residual += imbalance_u.abs();
scale += diagonal_u.abs();
}
}
for j in 1..ny {
for i in 0..nx {
let cv =
self.compute_v_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
let ap = cv.center * alpha;
let source = cv.source - (1.0 - alpha) * cv.center * flow_field.v_old[(j, i)];
let diagonal_v = ap * flow_field.v[(j, i)];
let east = if i + 1 < nx {
cv.east * flow_field.v[(j, i + 1)]
} else {
0.0
};
let west = if i > 0 {
cv.west * flow_field.v[(j, i - 1)]
} else {
0.0
};
let imbalance_v = diagonal_v
- (source
+ cv.east * flow_field.v[(j, i + 1)]
+ cv.west * flow_field.v[(j, i - 1)]
+ east
+ west
+ cv.north * flow_field.v[(j + 1, i)]
+ cv.south * flow_field.v[(j - 1, i)]);
residual += imbalance_v.abs();
@@ -430,20 +502,43 @@ impl SimpleSolver {
) -> CfdResult<()> {
let (nx, ny, _, _) = flow_field.grid_info();
// For each u-velocity point (face-centered)
for j in 1..ny - 1 {
// Every u row is an unknown.
//
// `u[(j, i)]` sits at `y = (j + 0.5) dy`, which is strictly interior
// for every `j`, so there is no u row on a horizontal wall to hold a
// boundary value. Only the faces `i = 0` and `i = nx` lie on a domain
// boundary, which is why the `i` range stops short of them and the `j`
// range does not stop at all.
//
// This previously swept `1..ny - 1`, freezing the two near-wall rows
// and treating whatever was stored there as a boundary condition —
// imposing the wall half a cell inside the domain.
for j in 0..ny {
for i in 1..nx {
// u goes from 1 to nx-1 for interior
// Discretize u-momentum equation at (i, j)
let coeffs =
self.compute_u_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
// The north and south coefficients are zero on a near-wall row,
// where the wall's contribution is already in `source`; the
// guards keep the index off the end of the array.
let north = if j + 1 < ny {
coeffs.north * flow_field.u[(j + 1, i)]
} else {
0.0
};
let south = if j > 0 {
coeffs.south * flow_field.u[(j - 1, i)]
} else {
0.0
};
// Solve for new u velocity using Gauss-Seidel
let u_new = (coeffs.source
+ coeffs.east * flow_field.u[(j, i + 1)]
+ coeffs.west * flow_field.u[(j, i - 1)]
+ coeffs.north * flow_field.u[(j + 1, i)]
+ coeffs.south * flow_field.u[(j - 1, i)])
+ north
+ south)
/ coeffs.center;
flow_field.u[(j, i)] = u_new;
@@ -465,18 +560,30 @@ impl SimpleSolver {
) -> CfdResult<()> {
let (nx, ny, _, _) = flow_field.grid_info();
// For each v-velocity point (face-centered)
// Every v column is an unknown, mirroring the u sweep: `v[(j, i)]` sits
// at `x = (i + 0.5) dx`, strictly interior for every `i`, and only the
// faces `j = 0` and `j = ny` lie on a domain boundary.
for j in 1..ny {
// v goes from 1 to ny-1 for interior
for i in 1..nx - 1 {
for i in 0..nx {
// Discretize v-momentum equation at (i, j)
let coeffs =
self.compute_v_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
let east = if i + 1 < nx {
coeffs.east * flow_field.v[(j, i + 1)]
} else {
0.0
};
let west = if i > 0 {
coeffs.west * flow_field.v[(j, i - 1)]
} else {
0.0
};
// Solve for new v velocity using Gauss-Seidel
let v_new = (coeffs.source
+ coeffs.east * flow_field.v[(j, i + 1)]
+ coeffs.west * flow_field.v[(j, i - 1)]
+ east
+ west
+ coeffs.north * flow_field.v[(j + 1, i)]
+ coeffs.south * flow_field.v[(j - 1, i)])
/ coeffs.center;
@@ -501,10 +608,22 @@ impl SimpleSolver {
// correction twice.
flow_field.p_prime.fill(0.0);
// Setup pressure correction equation: ∇²p' = ρ/Δt * ∇·u*
// Continuity is enforced on EVERY cell.
//
// The domain is tiled by cells; there is no such thing as a cell that
// does not have to conserve mass. Restricting this to `1..nx - 1` left
// the outer ring of cells with no continuity equation at all, so
// nothing ever removed their divergence — 3.1e-2 on the ring against
// 2.7e-3 in the interior on a converged 32x32 manufactured solve, with
// the resulting pressure error *growing* under refinement.
//
// This only became possible once the momentum sweeps stopped freezing
// the near-wall lines. Extending continuity first leaves ring cells
// whose faces are all prescribed — no correctable face, no solution —
// and it duly broke convergence when tried in that order.
let mut mass_imbalance: f64 = 0.0;
for j in 1..ny - 1 {
for i in 1..nx - 1 {
for j in 0..ny {
for i in 0..nx {
// Compute mass source (continuity equation residual)
let mass_source = self.compute_mass_source(flow_field, i, j, dx, dy, rho)?;
flow_field.sp[(j, i)] = mass_source;
@@ -517,37 +636,60 @@ impl SimpleSolver {
// correction divides by, so the two remain each other's inverse.
// Cached because `a_p` depends on the velocity field, which does not
// change during the inner sweeps.
let mut coefficients = Vec::with_capacity((nx - 2) * (ny - 2));
for j in 1..ny - 1 {
for i in 1..nx - 1 {
//
// A coefficient is zero exactly when the face it crosses is a genuine
// domain boundary — where the velocity is prescribed and therefore not
// correctable. Every cell keeps at least two correctable faces, so
// every cell can be made divergence-free.
let mut coefficients = Vec::with_capacity(nx * ny);
for j in 0..ny {
for i in 0..nx {
let volume = dx * dy;
let d_east = volume
/ self.compute_u_momentum_center_coefficient(
flow_field,
i + 1,
j,
dx,
dy,
rho,
)?;
let d_west = volume
/ self.compute_u_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
let d_north = volume
/ self.compute_v_momentum_center_coefficient(
flow_field,
i,
j + 1,
dx,
dy,
rho,
)?;
let d_south = volume
/ self.compute_v_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
let ae = if i + 1 == nx {
0.0
} else {
let d = volume
/ self.compute_u_momentum_center_coefficient(
flow_field,
i + 1,
j,
dx,
dy,
rho,
)?;
rho * d * dy / dx
};
let aw = if i == 0 {
0.0
} else {
let d = volume
/ self
.compute_u_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
rho * d * dy / dx
};
let an = if j + 1 == ny {
0.0
} else {
let d = volume
/ self.compute_v_momentum_center_coefficient(
flow_field,
i,
j + 1,
dx,
dy,
rho,
)?;
rho * d * dx / dy
};
let as_ = if j == 0 {
0.0
} else {
let d = volume
/ self
.compute_v_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
rho * d * dx / dy
};
let ae = rho * d_east * dy / dx;
let aw = rho * d_west * dy / dx;
let an = rho * d_north * dx / dy;
let as_ = rho * d_south * dx / dy;
coefficients.push(MomentumEquationCoeffs {
center: ae + aw + an + as_,
east: ae,
@@ -564,8 +706,8 @@ impl SimpleSolver {
// Inner pressure correction iterations
let mut residual = 0.0;
for j in 1..ny - 1 {
for i in 1..nx - 1 {
for j in 0..ny {
for i in 0..nx {
// Anchor one cell to fix the pressure level.
//
// With velocity prescribed on every boundary the pressure
@@ -582,19 +724,45 @@ impl SimpleSolver {
// that include the boundaries, so it is not required to sum
// to zero, and subtracting its mean injects a spurious
// source into every cell. Tried; it diverged.
//
// Dropping this one cell's continuity equation is legitimate
// now that the equation covers the whole domain: the sum of
// the sources over all cells telescopes to the net flux
// through the domain boundary, which is zero for a closed
// box, so the system has rank `n - 1` and exactly one
// equation is redundant.
if i == 1 && j == 1 {
flow_field.p_prime[(j, i)] = 0.0;
continue;
}
let coeffs = &coefficients[(j - 1) * (nx - 2) + (i - 1)];
let coeffs = &coefficients[j * nx + i];
let p_new = (flow_field.sp[(j, i)]
+ coeffs.east * flow_field.p_prime[(j, i + 1)]
+ coeffs.west * flow_field.p_prime[(j, i - 1)]
+ coeffs.north * flow_field.p_prime[(j + 1, i)]
+ coeffs.south * flow_field.p_prime[(j - 1, i)])
/ coeffs.center;
// A zero coefficient still guards its index: the neighbour
// it refers to is outside the domain.
let east = if i + 1 < nx {
coeffs.east * flow_field.p_prime[(j, i + 1)]
} else {
0.0
};
let west = if i > 0 {
coeffs.west * flow_field.p_prime[(j, i - 1)]
} else {
0.0
};
let north = if j + 1 < ny {
coeffs.north * flow_field.p_prime[(j + 1, i)]
} else {
0.0
};
let south = if j > 0 {
coeffs.south * flow_field.p_prime[(j - 1, i)]
} else {
0.0
};
let p_new =
(flow_field.sp[(j, i)] + east + west + north + south) / coeffs.center;
let correction = p_new - flow_field.p_prime[(j, i)];
residual += correction * correction;
@@ -631,27 +799,27 @@ impl SimpleSolver {
let (nx, ny, dx, dy) = flow_field.grid_info();
let rho = self.config.density;
// Correct u-velocities
for j in 1..ny - 1 {
// Correct every face the pressure equation treated as correctable —
// which is every face that is not on a domain boundary. The ranges must
// match `pressure_correction_step`'s coefficients exactly, or the
// divergence the pressure correction was constructed to remove is not
// the divergence that gets removed.
for j in 0..ny {
for i in 1..nx {
if i > 0 && i < nx {
let dp_dx = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j, i - 1)]) / dx;
let ap_u =
self.compute_u_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
flow_field.u[(j, i)] = flow_field.u_star[(j, i)] - (dx * dy / ap_u) * dp_dx;
}
let dp_dx = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j, i - 1)]) / dx;
let ap_u =
self.compute_u_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
flow_field.u[(j, i)] = flow_field.u_star[(j, i)] - (dx * dy / ap_u) * dp_dx;
}
}
// Correct v-velocities
for j in 1..ny {
for i in 1..nx - 1 {
if j > 0 && j < ny {
let dp_dy = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j - 1, i)]) / dy;
let ap_v =
self.compute_v_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
flow_field.v[(j, i)] = flow_field.v_star[(j, i)] - (dx * dy / ap_v) * dp_dy;
}
for i in 0..nx {
let dp_dy = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j - 1, i)]) / dy;
let ap_v =
self.compute_v_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
flow_field.v[(j, i)] = flow_field.v_star[(j, i)] - (dx * dy / ap_v) * dp_dy;
}
}
@@ -773,6 +941,8 @@ impl SimpleSolver {
dx: f64,
dy: f64,
) -> CfdResult<MomentumEquationCoeffs> {
let (_nx, ny, _, _) = flow_field.grid_info();
// Compute effective viscosity (molecular + turbulent)
let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
@@ -805,14 +975,61 @@ impl SimpleSolver {
// scheme reduced to diffusion plus a spurious diagonal term.
let fe = rho * 0.5 * (flow_field.u[(j, i)] + flow_field.u[(j, i + 1)]) * dy;
let fw = rho * 0.5 * (flow_field.u[(j, i - 1)] + flow_field.u[(j, i)]) * dy;
let fn_ = rho * 0.5 * (flow_field.v[(j + 1, i - 1)] + flow_field.v[(j + 1, i)]) * dx;
let fs = rho * 0.5 * (flow_field.v[(j, i - 1)] + flow_field.v[(j, i)]) * dx;
// Rows `j = 0` and `j = ny - 1` are *not* boundaries — every u row sits
// at `y = (j + 0.5) dy`, strictly inside the domain. They are near-wall
// interior lines whose control volume happens to have the wall for its
// south (respectively north) face.
//
// Two things change at such a face and nothing else does. A solid wall
// passes no mass, so the convective flux through it is zero whatever
// the stored normal velocity happens to be. And the node on the far
// side of the face is the wall itself, half a cell away rather than a
// full cell, so the conductance is `mu A / (dy/2)` — twice the interior
// value — and the value there is the wall's own tangential velocity,
// which is data rather than an unknown. Data belongs in the source, so
// the returned neighbour coefficient is zero while `a_p` still carries
// the conductance.
//
// Freezing these rows instead, as this sweep previously did, imposes
// the wall value half a cell inside the domain and leaves the outer
// ring of cells with too few correctable faces.
let south_is_wall = j == 0;
let north_is_wall = j + 1 == ny;
let fn_ = if north_is_wall {
0.0
} else {
rho * 0.5 * (flow_field.v[(j + 1, i - 1)] + flow_field.v[(j + 1, i)]) * dx
};
let fs = if south_is_wall {
0.0
} else {
rho * 0.5 * (flow_field.v[(j, i - 1)] + flow_field.v[(j, i)]) * dx
};
// Compute coefficients with upwind scheme
let ae = gamma_e + f64::max(-fe, 0.0);
let aw = gamma_w + f64::max(fw, 0.0);
let an = gamma_n + f64::max(-fn_, 0.0);
let as_ = gamma_s + f64::max(fs, 0.0);
// `_p` enters the diagonal; `_nb` multiplies a stored neighbour and is
// zero at a wall, where the contribution goes to `wall_source` instead.
let gamma_wall = mu_eff * dx / (0.5 * dy);
let mut wall_source = 0.0;
let (an, an_nb) = if north_is_wall {
wall_source += gamma_wall * self.u_wall(flow_field, i, j, ny as f64 * dy, dx);
(gamma_wall, 0.0)
} else {
let a = gamma_n + f64::max(-fn_, 0.0);
(a, a)
};
let (as_, as_nb) = if south_is_wall {
wall_source += gamma_wall * self.u_wall(flow_field, i, j, 0.0, dx);
(gamma_wall, 0.0)
} else {
let a = gamma_s + f64::max(fs, 0.0);
(a, a)
};
// Transient term. Zero for a steady solve: standard SIMPLE has no
// pseudo-time term, and keeping one makes the converged answer depend
@@ -848,16 +1065,23 @@ impl SimpleSolver {
// using an unrelaxed `a_p` while the velocities have been relaxed.
let alpha = self.parameters.velocity_relaxation;
let ap = ap_unrelaxed / alpha;
// `u_source_term` was computed and then never added — the x-momentum
// equation carried no body force at all, while the y-momentum equation
// carried its own. Any manufactured solution was therefore imposed on
// one component and not the other, which is why `u` came out markedly
// further from exact than `v` on the same mesh.
let source = pressure_gradient
+ time_term
+ wall_source
+ self.u_source_term(i, j, dx, dy)
+ (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.u_old[(j, i)];
Ok(MomentumEquationCoeffs {
center: ap,
east: ae,
west: aw,
north: an,
south: as_,
north: an_nb,
south: as_nb,
source,
})
}
@@ -874,6 +1098,8 @@ impl SimpleSolver {
dx: f64,
dy: f64,
) -> CfdResult<MomentumEquationCoeffs> {
let (nx, _ny, _, _) = flow_field.grid_info();
// Compute effective viscosity (molecular + turbulent)
let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
@@ -891,14 +1117,46 @@ impl SimpleSolver {
// all four faces cannot represent transport.
let fn_ = rho * 0.5 * (flow_field.v[(j, i)] + flow_field.v[(j + 1, i)]) * dx;
let fs = rho * 0.5 * (flow_field.v[(j - 1, i)] + flow_field.v[(j, i)]) * dx;
let fe = rho * 0.5 * (flow_field.u[(j - 1, i + 1)] + flow_field.u[(j, i + 1)]) * dy;
let fw = rho * 0.5 * (flow_field.u[(j - 1, i)] + flow_field.u[(j, i)]) * dy;
let ae = gamma_e + f64::max(-fe, 0.0);
let aw = gamma_w + f64::max(fw, 0.0);
// Columns `i = 0` and `i = nx - 1` are near-wall interior lines, not
// boundaries: every v column sits at `x = (i + 0.5) dx`. The west and
// east walls bound them at half-cell distance. See the u-momentum
// routine for why that changes the conductance and kills the
// convective flux, and nothing else.
let west_is_wall = i == 0;
let east_is_wall = i + 1 == nx;
let fe = if east_is_wall {
0.0
} else {
rho * 0.5 * (flow_field.u[(j - 1, i + 1)] + flow_field.u[(j, i + 1)]) * dy
};
let fw = if west_is_wall {
0.0
} else {
rho * 0.5 * (flow_field.u[(j - 1, i)] + flow_field.u[(j, i)]) * dy
};
let an = gamma_n + f64::max(-fn_, 0.0);
let as_ = gamma_s + f64::max(fs, 0.0);
let gamma_wall = mu_eff * dy / (0.5 * dx);
let mut wall_source = 0.0;
let (ae, ae_nb) = if east_is_wall {
wall_source += gamma_wall * self.v_wall(flow_field, i, j, nx as f64 * dx, dy);
(gamma_wall, 0.0)
} else {
let a = gamma_e + f64::max(-fe, 0.0);
(a, a)
};
let (aw, aw_nb) = if west_is_wall {
wall_source += gamma_wall * self.v_wall(flow_field, i, j, 0.0, dy);
(gamma_wall, 0.0)
} else {
let a = gamma_w + f64::max(fw, 0.0);
(a, a)
};
let ap0 = if self.parameters.steady {
0.0
} else {
@@ -917,13 +1175,14 @@ impl SimpleSolver {
let ap = ap_unrelaxed / alpha;
let source = pressure_gradient
+ time_term
+ wall_source
+ self.v_source_term(i, j, dx, dy)
+ (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.v_old[(j, i)];
Ok(MomentumEquationCoeffs {
center: ap,
east: ae,
west: aw,
east: ae_nb,
west: aw_nb,
north: an,
south: as_,
source,