diff --git a/Components/Cylinder.cs b/Components/Cylinder.cs index 3f14d2f..0ba1625 100644 --- a/Components/Cylinder.cs +++ b/Components/Cylinder.cs @@ -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; diff --git a/Core/BoundarySystem.cs b/Core/BoundarySystem.cs index 5264250..aa30a70 100644 --- a/Core/BoundarySystem.cs +++ b/Core/BoundarySystem.cs @@ -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; } } diff --git a/Core/Pipesystem.cs b/Core/Pipesystem.cs index 6a927e8..7113efb 100644 --- a/Core/Pipesystem.cs +++ b/Core/Pipesystem.cs @@ -16,23 +16,28 @@ namespace FluidSim.Core private readonly int _allCells; // total allocated (padded to Vector.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.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(_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.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(_rho, iCell); - var rhouOld = new Vector(_rhou, iCell); - var EOld = new Vector(_E, iCell); - var YOld = new Vector(_Y, iCell); + int start = _pipeStart[p]; + int end = _pipeEnd[p]; // exclusive + int len = end - start; + if (len == 0) continue; - var fluxM_L = new Vector(_fluxM, iCell); - var fluxP_L = new Vector(_fluxP, iCell); - var fluxE_L = new Vector(_fluxE, iCell); - var fluxY_L = new Vector(_fluxY, iCell); - - var fluxM_R = new Vector(_fluxM, iCell + 1); - var fluxP_R = new Vector(_fluxP, iCell + 1); - var fluxE_R = new Vector(_fluxE, iCell + 1); - var fluxY_R = new Vector(_fluxY, iCell + 1); - - var dtdx = new Vector(dt) / new Vector(_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(_dampingFactors, iCell); - if (_applyRelax) + // ------- Left boundary cell (i = start) ------ { - var ambRef = new Vector(_ambientEnergyRef); - var relax = new Vector(_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(1e-12f)); - var kinNew = 0.5f * rhouNew * rhouNew / rhoNew; - var eMin = new Vector(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.Zero), Vector.One); - yNew.CopyTo(_Y, iCell); + // Vectorised path for interior cells (if available) + for (; iCell <= iEnd - vecSize; iCell += vecSize) + { + var rhoOld = new Vector(_rho, iCell); + var rhouOld = new Vector(_rhou, iCell); + var EOld = new Vector(_E, iCell); + var YOld = new Vector(_Y, iCell); + + var fluxM_L = new Vector(_fluxM, iCell); + var fluxP_L = new Vector(_fluxP, iCell); + var fluxE_L = new Vector(_fluxE, iCell); + var fluxY_L = new Vector(_fluxY, iCell); + + var fluxM_R = new Vector(_fluxM, iCell + 1); + var fluxP_R = new Vector(_fluxP, iCell + 1); + var fluxE_R = new Vector(_fluxE, iCell + 1); + var fluxY_R = new Vector(_fluxY, iCell + 1); + + var dtdx = new Vector(dt) / new Vector(_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(_dampingFactors, iCell); + if (_applyRelax) + { + var ambRef = new Vector(_ambientEnergyRef); + var relax = new Vector(_relaxFactors, iCell); + ENew = ambRef + (ENew - ambRef) * relax; + } + + rhoNew = Vector.Max(rhoNew, new Vector(1e-12f)); + var kinNew = 0.5f * rhouNew * rhouNew / rhoNew; + var eMin = new Vector(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.Zero), Vector.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() { diff --git a/Core/Solver.cs b/Core/Solver.cs index 1559fa2..00d24e4 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -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++) diff --git a/Program.cs b/Program.cs index 784dd3b..d9469ae 100644 --- a/Program.cs +++ b/Program.cs @@ -50,7 +50,7 @@ public class Program { var window = CreateWindow(); LoadFont(); - _scenario = new HelmholtzScenario(); + _scenario = new SingleCylScenario(); _scenario.Initialize(SampleRate); _lastThrottleUpdateTime = 0.0f; diff --git a/Scenarios/SingleCylScenario.cs b/Scenarios/SingleCylScenario.cs index f149fac..89c20b0 100644 --- a/Scenarios/SingleCylScenario.cs +++ b/Scenarios/SingleCylScenario.cs @@ -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)