693 lines
29 KiB
C#
693 lines
29 KiB
C#
using System;
|
||
using System.Diagnostics;
|
||
using System.Numerics;
|
||
|
||
namespace FluidSim.Core
|
||
{
|
||
public class PipeSystem
|
||
{
|
||
// ---------- Master arrays ----------
|
||
private float[] _rho, _rhou, _E, _Y;
|
||
private readonly float[] _area;
|
||
private readonly float[] _dx;
|
||
private readonly int[] _pipeStart;
|
||
private readonly int[] _pipeEnd;
|
||
private readonly int _totalCells; // original cell count (visible)
|
||
private readonly int _allCells; // total allocated (padded to Vector<float>.Count)
|
||
private readonly int _pipeCount;
|
||
|
||
// Derived state – _p is kept for visualization
|
||
private float[] _p;
|
||
|
||
// Flux arrays for faces INTERNAL to a single pipe (size = _allCells + 1)
|
||
// Only valid for faces that are NOT pipe boundaries.
|
||
private float[] _fluxM, _fluxP, _fluxE, _fluxY;
|
||
|
||
// Per‑pipe boundary flux buffers (size = _pipeCount)
|
||
private float[] _leftFluxM, _leftFluxP, _leftFluxE, _leftFluxY;
|
||
private float[] _rightFluxM, _rightFluxP, _rightFluxE, _rightFluxY;
|
||
|
||
// Damping and relaxation
|
||
private float[] _dampingFactors;
|
||
private float[] _relaxFactors;
|
||
private bool _applyDamping;
|
||
private bool _applyRelax;
|
||
|
||
// Ghost buffer (per‑pipe ghost states)
|
||
private readonly GhostBuffer _ghost;
|
||
|
||
// Precomputed flag: true if a face is a pipe boundary (start or end)
|
||
private readonly bool[] _isPipeBoundaryFace;
|
||
|
||
// ---------- Physical constants ----------
|
||
private const float Gamma = 1.4f;
|
||
private const float Gm1 = 0.4f;
|
||
private const float Gm1Inv = 1f / Gm1; // 2.5
|
||
private const float GammaOverGm1 = Gamma / Gm1; // 3.5
|
||
private float _coeffBase;
|
||
private float _relaxRate;
|
||
private float _ambientPressure = 101325f;
|
||
private float _ambientEnergyRef;
|
||
|
||
public float DampingMultiplier
|
||
{
|
||
set
|
||
{
|
||
_coeffBase = 0.1f * value;
|
||
_applyDamping = _coeffBase != 0f;
|
||
}
|
||
}
|
||
public float EnergyRelaxationRate
|
||
{
|
||
set
|
||
{
|
||
_relaxRate = value;
|
||
_applyRelax = _relaxRate != 0f;
|
||
}
|
||
}
|
||
public float AmbientPressure
|
||
{
|
||
set
|
||
{
|
||
_ambientPressure = value;
|
||
_ambientEnergyRef = value * Gm1Inv;
|
||
}
|
||
}
|
||
|
||
// ---------- Profiling ----------
|
||
public bool EnableProfiling { get; set; }
|
||
private long _profFluxTicks;
|
||
private long _profUpdateTicks;
|
||
private long _profCallCount;
|
||
|
||
// ---------- Construction ----------
|
||
public PipeSystem(int totalCells, int[] pipeStart, int[] pipeEnd,
|
||
float[] area, float[] dx,
|
||
float initialRho, float initialU, float initialP)
|
||
{
|
||
_pipeStart = pipeStart;
|
||
_pipeEnd = pipeEnd;
|
||
_pipeCount = pipeStart.Length;
|
||
_totalCells = totalCells;
|
||
_area = area;
|
||
_dx = dx;
|
||
|
||
// Pad to SIMD width so all vectorized loops cover the whole data
|
||
int vecSize = Vector<float>.Count;
|
||
_allCells = totalCells % vecSize == 0 ? totalCells : totalCells + vecSize - (totalCells % vecSize);
|
||
|
||
_rho = new float[_allCells];
|
||
_rhou = new float[_allCells];
|
||
_E = new float[_allCells];
|
||
_Y = new float[_allCells];
|
||
_p = new float[_allCells]; // pressure for drawing
|
||
int faceCount = _allCells + 1;
|
||
_fluxM = new float[faceCount];
|
||
_fluxP = new float[faceCount];
|
||
_fluxE = new float[faceCount];
|
||
_fluxY = new float[faceCount];
|
||
|
||
// Per‑pipe boundary flux buffers
|
||
_leftFluxM = new float[_pipeCount];
|
||
_leftFluxP = new float[_pipeCount];
|
||
_leftFluxE = new float[_pipeCount];
|
||
_leftFluxY = new float[_pipeCount];
|
||
_rightFluxM = new float[_pipeCount];
|
||
_rightFluxP = new float[_pipeCount];
|
||
_rightFluxE = new float[_pipeCount];
|
||
_rightFluxY = new float[_pipeCount];
|
||
|
||
_dampingFactors = new float[_allCells];
|
||
_relaxFactors = new float[_allCells];
|
||
_applyDamping = _coeffBase != 0f;
|
||
_applyRelax = _relaxRate != 0f;
|
||
|
||
_ghost = new GhostBuffer(_pipeCount);
|
||
_ambientEnergyRef = initialP * Gm1Inv;
|
||
|
||
// Mark faces that coincide with a pipe boundary (start or end)
|
||
_isPipeBoundaryFace = new bool[faceCount];
|
||
for (int p = 0; p < _pipeCount; p++)
|
||
{
|
||
_isPipeBoundaryFace[_pipeStart[p]] = true;
|
||
_isPipeBoundaryFace[_pipeEnd[p]] = true;
|
||
}
|
||
|
||
// Initialize uniform state
|
||
float initE = initialP / (Gm1 * initialRho);
|
||
float rhoE = initialRho * initE + 0.5f * initialRho * initialU * initialU;
|
||
for (int i = 0; i < totalCells; i++)
|
||
{
|
||
_rho[i] = initialRho;
|
||
_rhou[i] = initialRho * initialU;
|
||
_E[i] = rhoE;
|
||
_Y[i] = 1f;
|
||
}
|
||
}
|
||
|
||
// ---------- Ghost setters (for BoundarySystem) ----------
|
||
public void SetGhostLeft(int pipeIndex, float rho, float u, float p, float y)
|
||
=> _ghost.Set(pipeIndex, true, rho, u, p, y);
|
||
public void SetGhostRight(int pipeIndex, float rho, float u, float p, float y)
|
||
=> _ghost.Set(pipeIndex, false, rho, u, p, y);
|
||
|
||
// ---------- Public read methods ----------
|
||
public int TotalCells => _totalCells;
|
||
public int PipeCount => _pipeCount;
|
||
public int GetPipeStart(int pipeIdx) => _pipeStart[pipeIdx];
|
||
public int GetPipeEnd(int pipeIdx) => _pipeEnd[pipeIdx];
|
||
public float GetCellPressure(int i) => _p[i];
|
||
public float GetCellDensity(int i) => _rho[i];
|
||
public float GetCellDx(int i) => _dx[i];
|
||
public float GetCellArea(int i) => _area[i];
|
||
public float GetCellVelocity(int i)
|
||
{
|
||
float rho = _rho[i];
|
||
return rho > 1e-12f ? _rhou[i] / rho : 0f;
|
||
}
|
||
public float GetCellAirFraction(int i) => _Y[i];
|
||
|
||
public (float rho, float u, float p) GetInteriorStateLeft(int pipeIdx)
|
||
{
|
||
int i = _pipeStart[pipeIdx];
|
||
float rho = _rho[i];
|
||
float rhou = _rhou[i];
|
||
float u = rhou / MathF.Max(rho, 1e-12f);
|
||
float p = Gm1 * (_E[i] - 0.5f * rhou * u);
|
||
return (rho, u, p);
|
||
}
|
||
public (float rho, float u, float p) GetInteriorStateRight(int pipeIdx)
|
||
{
|
||
int i = _pipeEnd[pipeIdx] - 1;
|
||
float rho = _rho[i];
|
||
float rhou = _rhou[i];
|
||
float u = rhou / MathF.Max(rho, 1e-12f);
|
||
float p = Gm1 * (_E[i] - 0.5f * rhou * u);
|
||
return (rho, u, p);
|
||
}
|
||
public float GetInteriorAirFractionLeft(int pipeIdx) => _Y[_pipeStart[pipeIdx]];
|
||
public float GetInteriorAirFractionRight(int pipeIdx) => _Y[_pipeEnd[pipeIdx] - 1];
|
||
|
||
public void SetCellState(int i, float rho, float u, float p, float y = 1f)
|
||
{
|
||
if (i < 0 || i >= _totalCells) return;
|
||
_rho[i] = rho;
|
||
_rhou[i] = rho * u;
|
||
_E[i] = p * Gm1Inv + 0.5f * rho * u * u;
|
||
_Y[i] = y;
|
||
}
|
||
|
||
// ---------- Main step ----------
|
||
public void SimulateStep(float dt)
|
||
{
|
||
long t0 = 0, t1 = 0;
|
||
if (EnableProfiling)
|
||
{
|
||
_profCallCount++;
|
||
t0 = Stopwatch.GetTimestamp();
|
||
}
|
||
|
||
ComputeFluxes(dt);
|
||
|
||
if (EnableProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profFluxTicks += (t1 - t0);
|
||
t0 = t1;
|
||
}
|
||
|
||
UpdateCells(dt);
|
||
|
||
if (EnableProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profUpdateTicks += (t1 - t0);
|
||
}
|
||
}
|
||
|
||
// ---------- Flux computation ----------
|
||
private void ComputeFluxes(float dt)
|
||
{
|
||
float fm, fp, fe;
|
||
int vecSize = Vector<float>.Count;
|
||
|
||
// ---- 1. Left ghost boundaries → per‑pipe buffers ----
|
||
for (int p = 0; p < _pipeCount; p++)
|
||
{
|
||
int idx = _pipeStart[p];
|
||
int ghostIdx = p * 2;
|
||
float rL = _ghost.Rho[ghostIdx];
|
||
float uL = _ghost.U[ghostIdx];
|
||
float pL = _ghost.P[ghostIdx];
|
||
float YL = _ghost.Y[ghostIdx];
|
||
float cL = MathF.Sqrt(Gamma * pL / MathF.Max(rL, 1e-12f));
|
||
|
||
float rR = _rho[idx], rhouR = _rhou[idx];
|
||
float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
|
||
float uR = rhouR * invRhoR;
|
||
float pR = Gm1 * (_E[idx] - 0.5f * rhouR * uR);
|
||
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
|
||
float YR = _Y[idx];
|
||
|
||
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
|
||
_leftFluxM[p] = fm; _leftFluxP[p] = fp; _leftFluxE[p] = fe;
|
||
|
||
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
|
||
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
|
||
_leftFluxY[p] = fy;
|
||
}
|
||
|
||
// ---- 2. Right ghost boundaries → per‑pipe buffers ----
|
||
for (int p = 0; p < _pipeCount; p++)
|
||
{
|
||
int idx = _pipeEnd[p] - 1;
|
||
int ghostIdx = p * 2 + 1;
|
||
float rR = _ghost.Rho[ghostIdx];
|
||
float uR = _ghost.U[ghostIdx];
|
||
float pR = _ghost.P[ghostIdx];
|
||
float YR = _ghost.Y[ghostIdx];
|
||
float cR = MathF.Sqrt(Gamma * pR / MathF.Max(rR, 1e-12f));
|
||
|
||
float rL = _rho[idx], rhouL = _rhou[idx];
|
||
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
|
||
float uL = rhouL * invRhoL;
|
||
float pL = Gm1 * (_E[idx] - 0.5f * rhouL * uL);
|
||
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
|
||
float YL = _Y[idx];
|
||
|
||
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
|
||
_rightFluxM[p] = fm; _rightFluxP[p] = fp; _rightFluxE[p] = fe;
|
||
|
||
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
|
||
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
|
||
_rightFluxY[p] = fy;
|
||
}
|
||
|
||
// ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
|
||
for (int face = 1; face < _totalCells; face++)
|
||
{
|
||
// Skip faces that belong to a pipe boundary (they are already handled)
|
||
if (_isPipeBoundaryFace[face])
|
||
continue;
|
||
|
||
// Try to vectorize a block of contiguous non‑boundary faces
|
||
if (face + vecSize - 1 < _totalCells)
|
||
{
|
||
bool canVectorize = true;
|
||
for (int f = face; f < face + vecSize; f++)
|
||
{
|
||
if (_isPipeBoundaryFace[f])
|
||
{
|
||
canVectorize = false;
|
||
break;
|
||
}
|
||
}
|
||
|
||
if (canVectorize)
|
||
{
|
||
// --- Vectorised block ---
|
||
var rhoL = new Vector<float>(_rho, face - 1);
|
||
var rhouL = new Vector<float>(_rhou, face - 1);
|
||
var EL = new Vector<float>(_E, face - 1);
|
||
var YL = new Vector<float>(_Y, face - 1);
|
||
var rhoR = new Vector<float>(_rho, face);
|
||
var rhouR = new Vector<float>(_rhou, face);
|
||
var ER = new Vector<float>(_E, face);
|
||
var YR = new Vector<float>(_Y, face);
|
||
|
||
var invRhoL = Vector<float>.One / Vector.Max(rhoL, new Vector<float>(1e-12f));
|
||
var invRhoR = Vector<float>.One / Vector.Max(rhoR, new Vector<float>(1e-12f));
|
||
var uL = rhouL * invRhoL;
|
||
var uR = rhouR * invRhoR;
|
||
var kinL = 0.5f * rhouL * uL;
|
||
var kinR = 0.5f * rhouR * uR;
|
||
var pL = Gm1 * (EL - kinL);
|
||
var pR = Gm1 * (ER - kinR);
|
||
var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
|
||
var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
|
||
|
||
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
|
||
var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
|
||
|
||
var FmL = rhoL * uL;
|
||
var FpL = rhoL * uL * uL + pL;
|
||
var FeL = (rhoL * ELs + pL) * uL;
|
||
|
||
var FmR = rhoR * uR;
|
||
var FpR = rhoR * uR * uR + pR;
|
||
var FeR = (rhoR * ERs + pR) * uR;
|
||
|
||
var absUL = Vector.Abs(uL);
|
||
var absUR = Vector.Abs(uR);
|
||
var alpha = Vector.Max(absUL + cL, absUR + cR);
|
||
|
||
var fmVec = 0.5f * (FmL + FmR) - 0.5f * alpha * (rhoR - rhoL);
|
||
var fpVec = 0.5f * (FpL + FpR) - 0.5f * alpha * (rhouR - rhouL);
|
||
var feVec = 0.5f * (FeL + FeR) - 0.5f * alpha * (rhoR * ERs - rhoL * ELs);
|
||
|
||
var fyL = FmL * YL;
|
||
var fyR = FmR * YR;
|
||
var fyVec = 0.5f * (fyL + fyR) - 0.5f * alpha * (rhoR * YR - rhoL * YL);
|
||
|
||
fmVec.CopyTo(_fluxM, face);
|
||
fpVec.CopyTo(_fluxP, face);
|
||
feVec.CopyTo(_fluxE, face);
|
||
fyVec.CopyTo(_fluxY, face);
|
||
|
||
face += vecSize - 1; // loop increment will add 1
|
||
continue;
|
||
}
|
||
}
|
||
|
||
// --- Scalar fallback for a single interior face ---
|
||
{
|
||
int iL = face - 1, iR = face;
|
||
float rL = _rho[iL], rhouL = _rhou[iL];
|
||
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
|
||
float uL = rhouL * invRhoL;
|
||
float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
|
||
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
|
||
float YL = _Y[iL];
|
||
|
||
float rR = _rho[iR], rhouR = _rhou[iR];
|
||
float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
|
||
float uR = rhouR * invRhoR;
|
||
float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
|
||
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
|
||
float YR = _Y[iR];
|
||
|
||
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
|
||
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
|
||
|
||
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
|
||
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
|
||
_fluxY[face] = fy;
|
||
}
|
||
}
|
||
}
|
||
|
||
// ---------- Cell update (per pipe, using correct boundary fluxes) ----------
|
||
private void UpdateCells(float dt)
|
||
{
|
||
int vecSize = Vector<float>.Count;
|
||
float dtRelax = -_relaxRate * dt;
|
||
|
||
// Precompute damping and relaxation factors globally
|
||
if (_applyDamping)
|
||
{
|
||
for (int i = 0; i < _totalCells; i++)
|
||
{
|
||
float rho = _rho[i];
|
||
_dampingFactors[i] = rho > 1e-12f
|
||
? MathF.Exp(-_coeffBase * dt / rho)
|
||
: 1f;
|
||
}
|
||
}
|
||
if (_applyRelax)
|
||
{
|
||
float relaxVal = MathF.Exp(dtRelax);
|
||
for (int i = 0; i < _totalCells; i++)
|
||
_relaxFactors[i] = relaxVal;
|
||
}
|
||
|
||
// Update each pipe separately
|
||
for (int p = 0; p < _pipeCount; p++)
|
||
{
|
||
int start = _pipeStart[p];
|
||
int end = _pipeEnd[p]; // exclusive
|
||
int len = end - start;
|
||
if (len == 0) continue;
|
||
|
||
// ------- Left boundary cell (i = start) ------
|
||
{
|
||
int i = start;
|
||
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
|
||
|
||
// left face: always the pipe's left boundary flux
|
||
float fluxM_L = _leftFluxM[p];
|
||
float fluxP_L = _leftFluxP[p];
|
||
float fluxE_L = _leftFluxE[p];
|
||
float fluxY_L = _leftFluxY[p];
|
||
|
||
// right face: depends on pipe length
|
||
float fluxM_R, fluxP_R, fluxE_R, fluxY_R;
|
||
if (len == 1)
|
||
{
|
||
// Only one cell: right face is the pipe's right boundary flux
|
||
fluxM_R = _rightFluxM[p];
|
||
fluxP_R = _rightFluxP[p];
|
||
fluxE_R = _rightFluxE[p];
|
||
fluxY_R = _rightFluxY[p];
|
||
}
|
||
else
|
||
{
|
||
// interior face (global flux at index i+1)
|
||
fluxM_R = _fluxM[i + 1];
|
||
fluxP_R = _fluxP[i + 1];
|
||
fluxE_R = _fluxE[i + 1];
|
||
fluxY_R = _fluxY[i + 1];
|
||
}
|
||
|
||
float dtdx = dt / _dx[i];
|
||
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
|
||
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
|
||
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
|
||
float rhoYOld = rhoOld * YOld;
|
||
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
|
||
|
||
if (_applyDamping) rhouNew *= _dampingFactors[i];
|
||
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
|
||
|
||
rhoNew = MathF.Max(rhoNew, 1e-12f);
|
||
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
|
||
float eMin = 100f * Gm1Inv + kin;
|
||
ENew = MathF.Max(ENew, eMin);
|
||
|
||
_rho[i] = rhoNew;
|
||
_rhou[i] = rhouNew;
|
||
_E[i] = ENew;
|
||
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
|
||
}
|
||
|
||
// ------- Interior cells (i = start+1 to end-2) ------
|
||
if (len > 2)
|
||
{
|
||
int iCell = start + 1;
|
||
int iEnd = end - 1; // exclusive upper bound
|
||
|
||
// Vectorised path for interior cells (if available)
|
||
for (; iCell <= iEnd - vecSize; iCell += vecSize)
|
||
{
|
||
var rhoOld = new Vector<float>(_rho, iCell);
|
||
var rhouOld = new Vector<float>(_rhou, iCell);
|
||
var EOld = new Vector<float>(_E, iCell);
|
||
var YOld = new Vector<float>(_Y, iCell);
|
||
|
||
var fluxM_L = new Vector<float>(_fluxM, iCell);
|
||
var fluxP_L = new Vector<float>(_fluxP, iCell);
|
||
var fluxE_L = new Vector<float>(_fluxE, iCell);
|
||
var fluxY_L = new Vector<float>(_fluxY, iCell);
|
||
|
||
var fluxM_R = new Vector<float>(_fluxM, iCell + 1);
|
||
var fluxP_R = new Vector<float>(_fluxP, iCell + 1);
|
||
var fluxE_R = new Vector<float>(_fluxE, iCell + 1);
|
||
var fluxY_R = new Vector<float>(_fluxY, iCell + 1);
|
||
|
||
var dtdx = new Vector<float>(dt) / new Vector<float>(_dx, iCell);
|
||
|
||
var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
|
||
var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
|
||
var ENew = EOld - dtdx * (fluxE_R - fluxE_L);
|
||
var rhoYOld = rhoOld * YOld;
|
||
var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
|
||
|
||
if (_applyDamping)
|
||
rhouNew *= new Vector<float>(_dampingFactors, iCell);
|
||
if (_applyRelax)
|
||
{
|
||
var ambRef = new Vector<float>(_ambientEnergyRef);
|
||
var relax = new Vector<float>(_relaxFactors, iCell);
|
||
ENew = ambRef + (ENew - ambRef) * relax;
|
||
}
|
||
|
||
rhoNew = Vector.Max(rhoNew, new Vector<float>(1e-12f));
|
||
var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
|
||
var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
|
||
ENew = Vector.Max(ENew, eMin);
|
||
|
||
rhoNew.CopyTo(_rho, iCell);
|
||
rhouNew.CopyTo(_rhou, iCell);
|
||
ENew.CopyTo(_E, iCell);
|
||
var yNew = rhoYNew / rhoNew;
|
||
yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
|
||
yNew.CopyTo(_Y, iCell);
|
||
}
|
||
|
||
// Scalar remainder for interior cells
|
||
for (; iCell < iEnd; iCell++)
|
||
{
|
||
float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
|
||
float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
|
||
float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
|
||
float dtdx = dt / _dx[iCell];
|
||
|
||
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
|
||
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
|
||
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
|
||
float rhoYOld = rhoOld * YOld;
|
||
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
|
||
|
||
if (_applyDamping) rhouNew *= _dampingFactors[iCell];
|
||
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell];
|
||
|
||
rhoNew = MathF.Max(rhoNew, 1e-12f);
|
||
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
|
||
float eMin = 100f * Gm1Inv + kin;
|
||
ENew = MathF.Max(ENew, eMin);
|
||
|
||
_rho[iCell] = rhoNew;
|
||
_rhou[iCell] = rhouNew;
|
||
_E[iCell] = ENew;
|
||
_Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
|
||
}
|
||
}
|
||
|
||
// ------- Right boundary cell (i = end-1, if len > 1) ------
|
||
if (len > 1)
|
||
{
|
||
int i = end - 1;
|
||
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
|
||
|
||
// left face
|
||
float fluxM_L, fluxP_L, fluxE_L, fluxY_L;
|
||
if (len == 2)
|
||
{
|
||
// Only two cells: left face is the pipe's left boundary flux
|
||
fluxM_L = _leftFluxM[p];
|
||
fluxP_L = _leftFluxP[p];
|
||
fluxE_L = _leftFluxE[p];
|
||
fluxY_L = _leftFluxY[p];
|
||
}
|
||
else
|
||
{
|
||
// interior face (global flux at i)
|
||
fluxM_L = _fluxM[i];
|
||
fluxP_L = _fluxP[i];
|
||
fluxE_L = _fluxE[i];
|
||
fluxY_L = _fluxY[i];
|
||
}
|
||
|
||
// right face: always the pipe's right boundary flux
|
||
float fluxM_R = _rightFluxM[p];
|
||
float fluxP_R = _rightFluxP[p];
|
||
float fluxE_R = _rightFluxE[p];
|
||
float fluxY_R = _rightFluxY[p];
|
||
|
||
float dtdx = dt / _dx[i];
|
||
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
|
||
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
|
||
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
|
||
float rhoYOld = rhoOld * YOld;
|
||
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
|
||
|
||
if (_applyDamping) rhouNew *= _dampingFactors[i];
|
||
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
|
||
|
||
rhoNew = MathF.Max(rhoNew, 1e-12f);
|
||
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
|
||
float eMin = 100f * Gm1Inv + kin;
|
||
ENew = MathF.Max(ENew, eMin);
|
||
|
||
_rho[i] = rhoNew;
|
||
_rhou[i] = rhouNew;
|
||
_E[i] = ENew;
|
||
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
|
||
}
|
||
}
|
||
|
||
// Recompute pressure for all cells (for visualization)
|
||
for (int i = 0; i < _totalCells; i++)
|
||
{
|
||
float rho = _rho[i];
|
||
float rhou = _rhou[i];
|
||
float u = rhou / MathF.Max(rho, 1e-12f);
|
||
_p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
|
||
}
|
||
}
|
||
|
||
// ---------- Scalar flux helpers ----------
|
||
private static void LaxFlux(float rL, float uL, float pL, float cL,
|
||
float rR, float uR, float pR, float cR,
|
||
out float fm, out float fp, out float fe)
|
||
{
|
||
float EL = pL * Gm1Inv / rL + 0.5f * uL * uL;
|
||
float ER = pR * Gm1Inv / rR + 0.5f * uR * uR;
|
||
float FmL = rL * uL;
|
||
float FpL = rL * uL * uL + pL;
|
||
float FeL = (rL * EL + pL) * uL;
|
||
float FmR = rR * uR;
|
||
float FpR = rR * uR * uR + pR;
|
||
float FeR = (rR * ER + pR) * uR;
|
||
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
|
||
fm = 0.5f * (FmL + FmR) - 0.5f * alpha * (rR - rL);
|
||
fp = 0.5f * (FpL + FpR) - 0.5f * alpha * (rR * uR - rL * uL);
|
||
fe = 0.5f * (FeL + FeR) - 0.5f * alpha * (rR * ER - rL * EL);
|
||
}
|
||
|
||
private static void ScalarFlux(float rL, float uL, float YL,
|
||
float rR, float uR, float YR,
|
||
float alpha, out float fy)
|
||
{
|
||
float FyL = rL * uL * YL;
|
||
float FyR = rR * uR * YR;
|
||
fy = 0.5f * (FyL + FyR) - 0.5f * alpha * (rR * YR - rL * YL);
|
||
}
|
||
|
||
public int GetRequiredSubSteps(float dtGlobal, float cflTarget = 0.8f)
|
||
{
|
||
float maxW = 0f;
|
||
for (int i = 0; i < _totalCells; i++)
|
||
{
|
||
float rho = MathF.Max(_rho[i], 1e-12f);
|
||
float u = MathF.Abs(_rhou[i] / rho);
|
||
float p = Gm1 * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
|
||
float c = MathF.Sqrt(Gamma * p / rho);
|
||
float w = u + c;
|
||
if (w > maxW) maxW = w;
|
||
}
|
||
maxW = MathF.Max(maxW, 1e-8f);
|
||
float minDx = _dx.Min(); // need using System.Linq;
|
||
return Math.Max(1, (int)MathF.Ceiling(dtGlobal * maxW / (cflTarget * minDx)));
|
||
}
|
||
|
||
// ---------- Profiling report ----------
|
||
public string GetProfileReport()
|
||
{
|
||
if (!EnableProfiling || _profCallCount == 0)
|
||
return "Pipe profiling disabled or no data.";
|
||
|
||
double freq = Stopwatch.Frequency;
|
||
long totalTicks = _profFluxTicks + _profUpdateTicks;
|
||
if (totalTicks == 0) return "No pipe profile data collected.";
|
||
|
||
double totalMs = totalTicks * 1000.0 / freq;
|
||
double avgCallUs = totalMs * 1000.0 / _profCallCount;
|
||
|
||
double fluxMs = _profFluxTicks * 1000.0 / freq;
|
||
double updateMs = _profUpdateTicks * 1000.0 / freq;
|
||
|
||
double fluxAvgUs = fluxMs * 1000.0 / _profCallCount;
|
||
double updateAvgUs = updateMs * 1000.0 / _profCallCount;
|
||
|
||
string report = $" Pipe kernel (over {_profCallCount} calls, total {totalMs:F2} ms, avg {avgCallUs:F2} µs/call):\n";
|
||
report += $" Fluxes (incl. primitives): {fluxMs:F2} ms ({_profFluxTicks * 100.0 / totalTicks:F1}%), avg {fluxAvgUs:F2} µs/call\n";
|
||
report += $" Update cells: {updateMs:F2} ms ({_profUpdateTicks * 100.0 / totalTicks:F1}%), avg {updateAvgUs:F2} µs/call\n";
|
||
|
||
_profFluxTicks = 0;
|
||
_profUpdateTicks = 0;
|
||
_profCallCount = 0;
|
||
|
||
return report;
|
||
}
|
||
}
|
||
} |