grillkol/FluidSim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

stable

Browse Source
This commit is contained in:
max
2026-06-09 17:49:11 +02:00
parent 1489f278dc
commit 21a62fb46e
6 changed files with 401 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Components/Cylinder.cs
View File

@@ -32,7 +32,7 @@ namespace FluidSim.Components
public float StoichiometricAFR = 14.7f;
public float FuelLowerHeatingValue = 44e6f;
public float EnergyVariationFraction = 0.05f;
public float MisfireProbability = 0.01f;
public float MisfireProbability = 0.0f;
public float CylinderWallArea = 0.02f;
public float HeatTransferCoefficient = 100f;
public float AmbientTemperature = 300f;

116
Core/BoundarySystem.cs
View File

@@ -165,7 +165,7 @@ namespace FluidSim.Core
continue;
}
// ---- Preliminary isentropic solution (used for face pressure if inertance is on) ----
// ---- Preliminary isentropic solution ----
float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
if (volP >= pipeP)
{
@@ -179,7 +179,7 @@ namespace FluidSim.Core
mdotEst = -mdotEst;
}
// ---- Compute ghost state ----
// ---- Compute final mass flow with limiters ----
float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
if (d.UseInertance)
@@ -189,11 +189,10 @@ namespace FluidSim.Core
float dp = volP - pipeP;
float Rlin = d.LossCoefficient;
// Forward Euler
float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
float mdotNew = d.CurrentMdot + dmdot_dt * dt;
// ---------- Symmetric flow limiters (existing) ----------
// Limit outflow from volume (if volume owner is Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
@@ -201,15 +200,19 @@ namespace FluidSim.Core
if (mdotNew < -maxOut) mdotNew = -maxOut;
}
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
// Limit inflow from pipe – pipe cell cannot be emptied in one step
{
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); // true cell area, not orifice
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
}
// ---------- Velocity clamp to Mach 0.9 ----------
// Velocity clamp to Mach 0.9
float rhoFacePrelim = mdotNew >= 0 ? volRho : pipeRho;
float uFacePrelim = MathF.Abs(mdotNew) / MathF.Max(rhoFacePrelim * area, 1e-12f);
float cUp = mdotNew >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
@@ -220,46 +223,64 @@ namespace FluidSim.Core
mdotNew = rhoFacePrelim * uFacePrelim * area * (mdotNew >= 0 ? 1f : -1f);
}
// NaN safety
if (float.IsNaN(mdotNew)) mdotNew = 0f;
d.CurrentMdot = mdotNew;
mdotFinal = mdotNew;
// Ghost state
rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
pFace = pFaceEst;
uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
}
else
{
// ---- Standard quasi‑steady orifice ----
// Standard quasi‑steady orifice
mdotFinal = mdotEst;
rhoFace = rhoFaceEst;
uFace = uFaceEst;
pFace = pFaceEst;
// Limit outflow from cavity
// Limit outflow from volume (if Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
if (mdotFinal > maxOut) mdotFinal = maxOut;
}
// ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
if (mdotFinal < 0)
{
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotFinal < -maxFromPipe)
mdotFinal = -maxFromPipe;
}
d.CurrentMdot = mdotFinal;
// Limit outflow from cylinder into pipe (positive mdot = volume → pipe)
if (mdotFinal > 0f && d.VolumePort?.Owner is Cylinder cyl)
{
float maxOut = cyl.Mass / dt;
if (mdotFinal > maxOut)
mdotFinal = maxOut;
}
}
// ---- Determine air fraction for ghost ----
// ---- Air fraction for ghost ----
if (mdotFinal >= 0)
{
airFracGhost = volAF;
}
else
{
airFracGhost = pipeAF;
if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
}
// ---- Apply sign convention for velocity ----
// ---- Sign convention for velocity ----
if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace;
else if (mdotFinal >= 0 && !d.IsLeftEnd) uFace = -uFace;
else if (mdotFinal < 0 && d.IsLeftEnd) uFace = -uFace;
@@ -271,18 +292,17 @@ namespace FluidSim.Core
else
_pipeSystem.SetGhostRight(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
// ---- Update volume port (mass flow: positive into volume) ----
// ---- Update volume port ----
if (d.VolumePort != null)
{
d.VolumePort.MassFlowRate = -mdotFinal;
// Set enthalpy of the stream entering the volume
if (-mdotFinal >= 0) // mass flowing into the volume (out of pipe)
if (-mdotFinal >= 0) // mass entering volume (out of pipe)
{
float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
d.VolumePort.SpecificEnthalpy = pipeH;
}
else // mass flowing out of the volume (into pipe)
else // mass leaving volume (into pipe)
{
d.VolumePort.SpecificEnthalpy = volH;
}
@@ -309,6 +329,7 @@ namespace FluidSim.Core
float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f));
float pAmb = d.AmbientPressure;
// Characteristic solution (isentropic expansion to ambient)
float Jplus = uInt + 2f * cInt / gm1;
float Jminus = uInt - 2f * cInt / gm1;
float s = pInt / MathF.Pow(rhoInt, gamma);
@@ -318,9 +339,14 @@ namespace FluidSim.Core
? (Jminus + 2f * cIso / gm1)
: (Jplus - 2f * cIso / gm1);
// Supersonic check
bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
float rhoGhost, uGhost, pGhost, afGhost;
if (!supersonic)
{
supersonic = d.IsLeftEnd ? (uIso <= -cIso) : (uIso >= cIso);
}
float rhoGhost, uGhost, pGhost, afGhost;
if (supersonic)
{
rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt;
@@ -332,15 +358,45 @@ namespace FluidSim.Core
afGhost = inflow ? 1f : afInt;
}
// ------- Mass flow limiter -------
int adjCell = d.IsLeftEnd
? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float cellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float area = d.PipeArea;
float mdotRaw = rhoGhost * uGhost * area; // positive out of pipe
if (d.IsLeftEnd) mdotRaw = -mdotRaw; // now positive = out of pipe
// Outflow limit
if (mdotRaw > 0 && mdotRaw * dt > cellMass)
{
mdotRaw = cellMass / dt;
}
// Inflow limit (allow up to 10× cell mass to avoid starving the pipe)
else if (mdotRaw < 0 && -mdotRaw * dt > 10f * cellMass)
{
mdotRaw = -10f * cellMass / dt;
}
// Recompute uGhost from the limited mdot, keeping rhoGhost, pGhost
float mdotMag = MathF.Abs(mdotRaw);
uGhost = mdotMag / MathF.Max(rhoGhost * area, 1e-12f);
if (d.IsLeftEnd)
uGhost = (mdotRaw >= 0f) ? -uGhost : uGhost;
else
uGhost = (mdotRaw >= 0f) ? uGhost : -uGhost;
// Apply ghost
if (d.IsLeftEnd)
_pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
else
_pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
float area = d.PipeArea;
float mdot = rhoGhost * uGhost * area;
if (d.IsLeftEnd) mdot = -mdot;
d.LastMassFlowRate = mdot;
d.LastMassFlowRate = mdotRaw;
d.LastFacePressure = pGhost;
}
}

428
Core/Pipesystem.cs
View File

@@ -16,23 +16,28 @@ namespace FluidSim.Core
private readonly int _allCells; // total allocated (padded to Vector<float>.Count)
private readonly int _pipeCount;
// Derived state – _p is kept for visualization, _c is gone
// Derived state – _p is kept for visualization
private float[] _p;
// Flux arrays (size = _allCells + 1)
// 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;
// Damping and relaxation (computed on‑the‑fly only if used)
// 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
// Ghost buffer (per‑pipe ghost states)
private readonly GhostBuffer _ghost;
// Wall mask – precomputed once
private readonly bool[] _isWallFace;
// 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;
@@ -102,6 +107,16 @@ namespace FluidSim.Core
_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;
@@ -110,18 +125,12 @@ namespace FluidSim.Core
_ghost = new GhostBuffer(_pipeCount);
_ambientEnergyRef = initialP * Gm1Inv;
// Pre‑compute wall face flags: each face that sits between two different pipes is a wall
_isWallFace = new bool[faceCount];
for (int f = 1; f < _totalCells; f++)
// Mark faces that coincide with a pipe boundary (start or end)
_isPipeBoundaryFace = new bool[faceCount];
for (int p = 0; p < _pipeCount; p++)
{
for (int p = 0; p < _pipeCount; p++)
{
if (f == _pipeEnd[p] && f < _totalCells)
{
_isWallFace[f] = true;
break;
}
}
_isPipeBoundaryFace[_pipeStart[p]] = true;
_isPipeBoundaryFace[_pipeEnd[p]] = true;
}
// Initialize uniform state
@@ -150,6 +159,7 @@ namespace FluidSim.Core
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];
@@ -215,13 +225,13 @@ namespace FluidSim.Core
}
}
// ---------- Flux computation: fuses primitive calculation and flux evaluation ----------
// ---------- Flux computation ----------
private void ComputeFluxes(float dt)
{
float fm, fp, fe;
int vecSize = Vector<float>.Count;
// ---- 1. Left ghost boundaries ----
// ---- 1. Left ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeStart[p];
@@ -239,22 +249,18 @@ namespace FluidSim.Core
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
float YR = _Y[idx];
// store pressure for cell idx
_p[idx] = pR;
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[idx] = fm; _fluxP[idx] = fp; _fluxE[idx] = 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);
_fluxY[idx] = fy;
_leftFluxY[p] = fy;
}
// ---- 2. Right ghost boundaries ----
// ---- 2. Right ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeEnd[p] - 1;
int face = idx + 1;
int ghostIdx = p * 2 + 1;
float rR = _ghost.Rho[ghostIdx];
float uR = _ghost.U[ghostIdx];
@@ -269,45 +275,35 @@ namespace FluidSim.Core
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
float YL = _Y[idx];
// store pressure for cell idx
_p[idx] = pL;
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = 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);
_fluxY[face] = fy;
_rightFluxY[p] = fy;
}
// ---- 3. Interior faces – vectorised SIMD ----
// ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
for (int face = 1; face < _totalCells; face++)
{
// Handle walls (rare) with scalar code
if (_isWallFace[face])
{
int iL = face - 1;
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);
_p[iL] = pL;
LaxFlux(rL, uL, pL, cL, rL, -uL, pL, cL, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
_fluxY[face] = 0f;
// Skip faces that belong to a pipe boundary (they are already handled)
if (_isPipeBoundaryFace[face])
continue;
}
// If the next vecSize faces contain a wall, fall back to scalar for this block
// Try to vectorize a block of contiguous non‑boundary faces
if (face + vecSize - 1 < _totalCells)
{
bool hasWall = false;
bool canVectorize = true;
for (int f = face; f < face + vecSize; f++)
if (_isWallFace[f]) { hasWall = true; break; }
{
if (_isPipeBoundaryFace[f])
{
canVectorize = false;
break;
}
}
if (!hasWall)
if (canVectorize)
{
// --- Vectorised block ---
var rhoL = new Vector<float>(_rho, face - 1);
@@ -330,11 +326,7 @@ namespace FluidSim.Core
var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
// Store pressures for visualisation (left cell of each face)
pL.CopyTo(_p, face - 1);
// Lax‑Friedrichs fluxes
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
var FmL = rhoL * uL;
@@ -362,50 +354,45 @@ namespace FluidSim.Core
feVec.CopyTo(_fluxE, face);
fyVec.CopyTo(_fluxY, face);
face += vecSize - 1; // loop increment will add 1, so we advance vecSize faces
face += vecSize - 1; // loop increment will add 1
continue;
}
}
// --- Scalar interior face (fallback) ---
// --- Scalar fallback for a single interior face ---
{
int iLf = face - 1, iRf = face;
float rLf = _rho[iLf], rhouLf = _rhou[iLf];
float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f));
float uLf = rhouLf * invRhoLf;
float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf);
float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf);
float YLf = _Y[iLf];
_p[iLf] = pLf;
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 rRf = _rho[iRf], rhouRf = _rhou[iRf];
float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f));
float uRf = rhouRf * invRhoRf;
float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf);
float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf);
float YRf = _Y[iRf];
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(rLf, uLf, pLf, cLf, rRf, uRf, pRf, cRf, out fm, out fp, out fe);
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(uLf) + cLf, MathF.Abs(uRf) + cRf);
ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy);
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;
}
}
// If damping/relaxation are active, compute the factors here (re-uses _dampingFactors/_relaxFactors arrays,
// but we no longer have a separate precompute pass). We compute them on demand in UpdateCells anyway?
// Actually UpdateCells multiplies by these factors; we can compute them there if needed.
}
// ---------- Cell update (unchanged core, but skips relaxation/damping when not needed) ----------
// ---------- Cell update (per pipe, using correct boundary fluxes) ----------
private void UpdateCells(float dt)
{
int vecSize = Vector<float>.Count;
float dtRelax = -_relaxRate * dt;
// Compute damping and relaxation factors if needed
// Precompute damping and relaxation factors globally
if (_applyDamping)
{
for (int i = 0; i < _totalCells; i++)
@@ -418,89 +405,217 @@ namespace FluidSim.Core
}
if (_applyRelax)
{
var relaxVal = MathF.Exp(dtRelax);
float relaxVal = MathF.Exp(dtRelax);
for (int i = 0; i < _totalCells; i++)
_relaxFactors[i] = relaxVal;
}
int iCell = 0;
for (; iCell <= _totalCells - vecSize; iCell += vecSize)
// Update each pipe separately
for (int p = 0; p < _pipeCount; p++)
{
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);
int start = _pipeStart[p];
int end = _pipeEnd[p]; // exclusive
int len = end - start;
if (len == 0) continue;
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)
// ------- Left boundary cell (i = start) ------
{
var ambRef = new Vector<float>(_ambientEnergyRef);
var relax = new Vector<float>(_relaxFactors, iCell);
ENew = ambRef + (ENew - ambRef) * relax;
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);
}
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);
// ------- Interior cells (i = start+1 to end-2) ------
if (len > 2)
{
int iCell = start + 1;
int iEnd = end - 1; // exclusive upper bound
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);
// 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);
}
}
// Scalar remainder (only a few cells)
for (; iCell < _totalCells; iCell++)
// Recompute pressure for all cells (for visualization)
for (int i = 0; i < _totalCells; i++)
{
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);
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 (used in boundaries and scalar fallback) ----------
// ---------- 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)
@@ -528,6 +643,23 @@ namespace FluidSim.Core
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()
{

3
Core/Solver.cs
View File

@@ -36,7 +36,8 @@ namespace FluidSim.Core
{
if (_pipeSystem == null || _boundarySystem == null) return;
int nSub = SubStepCount;
int nSub = _pipeSystem.GetRequiredSubSteps((float)_dt, 0.8f);
nSub = Math.Max(nSub, SubStepCount); // never go below fixed minimum
float dtSub = (float)(_dt / nSub);
for (int sub = 0; sub < nSub; sub++)

2
Program.cs
View File

@@ -50,7 +50,7 @@ public class Program
{
var window = CreateWindow();
LoadFont();
_scenario = new HelmholtzScenario();
_scenario = new SingleCylScenario();
_scenario.Initialize(SampleRate);
_lastThrottleUpdateTime = 0.0f;

42
Scenarios/SingleCylScenario.cs
View File

@@ -30,7 +30,7 @@ namespace FluidSim.Tests
private double dt;
private int stepCount;
public float MaxThrottleArea = 1e-4f; // 1 cm²
public float MaxThrottleArea = 100e-4f; // 1 cm²
// pipe area for open end calculations
private float pipeArea;
@@ -41,9 +41,9 @@ namespace FluidSim.Tests
// ---- Crankshaft ----
crankshaft = new Crankshaft(600);
crankshaft.Inertia = 0.2f;
crankshaft.Inertia = 0.05f;
crankshaft.FrictionConstant = 2f;
crankshaft.FrictionViscous = 0.04f;
crankshaft.FrictionViscous = 0.01f;
// ---- Cylinder ----
float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 9.2f;
@@ -77,12 +77,12 @@ namespace FluidSim.Tests
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
1.225f, 0f, 101325f);
pipeSystem.DampingMultiplier = 0.5f;
pipeSystem.EnergyRelaxationRate = 0f;
pipeSystem.DampingMultiplier = 1.0f;
pipeSystem.EnergyRelaxationRate = 0.5f;
pipeSystem.AmbientPressure = 101325f;
// ---- Volumes ----
intakePlenum = new Volume0D(5e-6f, 101325f, 300f); // 5 mL
intakePlenum = new Volume0D(100e-6f, 101325f, 300f); // 100 mL
plenumInlet = intakePlenum.CreatePort();
plenumOutlet = intakePlenum.CreatePort();
exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // 10 mL (unused but present)
@@ -130,8 +130,8 @@ namespace FluidSim.Tests
solver.AddComponent(exhaustCollector);
// ---- Sound ----
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 1f };
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 1f };
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
reverb = new OutdoorExhaustReverb(sampleRate);
stepCount = 0;
@@ -163,11 +163,31 @@ namespace FluidSim.Tests
if (stepCount % 1000 == 0)
{
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
Console.WriteLine($"Step {stepCount}, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar");
Console.WriteLine($"intake flow: {intakeFlow:F12}, exhaust flow: {exhaustFlow:F16}");
float crankDeg = crankshaft.CrankAngle; // public property on Cylinder
Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1} deg, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar");
Console.WriteLine($" intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
// Pipe 0 (intake before throttle)
var (r0L, u0L, p0L) = pipeSystem.GetInteriorStateLeft(0);
var (r0R, u0R, p0R) = pipeSystem.GetInteriorStateRight(0);
Console.WriteLine($" Pipe0 L: rho={r0L:F4} u={u0L:F3} p={p0L/1e5:F3}bar | R: rho={r0R:F4} u={u0R:F3} p={p0R/1e5:F3}bar");
// Pipe 1 (runner)
var (r1L, u1L, p1L) = pipeSystem.GetInteriorStateLeft(1);
var (r1R, u1R, p1R) = pipeSystem.GetInteriorStateRight(1);
Console.WriteLine($" Pipe1 L: rho={r1L:F4} u={u1L:F3} p={p1L/1e5:F3}bar | R: rho={r1R:F4} u={u1R:F3} p={p1R/1e5:F3}bar");
// Pipe 2 (exhaust)
var (r2L, u2L, p2L) = pipeSystem.GetInteriorStateLeft(2);
var (r2R, u2R, p2R) = pipeSystem.GetInteriorStateRight(2);
Console.WriteLine($" Pipe2 L: rho={r2L:F4} u={u2L:F3} p={p2L/1e5:F3}bar | R: rho={r2R:F4} u={u2R:F3} p={p2R/1e5:F3}bar");
// Plenum and cylinder mass
Console.WriteLine($" Plenum P={intakePlenum.Pressure/1e5:F3}bar, mass={intakePlenum.Mass:E4} kg");
Console.WriteLine($" Cyl mass={cylinder.Mass:E4} kg");
}
return reverb.Process(intakeDry);
return reverb.Process(intakeDry + exhaustDry);
}
public override void Draw(RenderWindow target)