config tuning

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;

Core/BoundarySystem.cs

@@ -25,8 +25,7 @@ namespace FluidSim.Core
             public float EffectiveLength;
             public float CurrentMdot;       // kg/s, positive = volume → pipe
 
-            // --- Loss coefficient (linear resistance) – inertance only ---
+            // --- Loss coefficient (linear resistance) ---
-            // If 0 when UseInertance is true, a stable default is auto‑computed at runtime.
             public float LossCoefficient;   // N·s/m⁵  or  kg/(m⁴·s)
         }
 
@@ -58,10 +57,9 @@ namespace FluidSim.Core
         public int OpenEndCount { get; private set; }
 
         // ---------- Add orifice (no inertance) ----------
-        // Simple isentropic nozzle – no built‑in loss. For dissipation use pipe damping
-        // or the inertance model if you need a damped resonator.
         public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
-                               int areaIndex, float dischargeCoeff = 1f)
+                               int areaIndex, float dischargeCoeff = 1f,
+                               float lossCoefficient = 0f)
         {
             _orifices[OrificeCount] = new OrificeDesc
             {
@@ -73,24 +71,22 @@ namespace FluidSim.Core
                 UseInertance = false,
                 EffectiveLength = 0f,
                 CurrentMdot = 0f,
-                LossCoefficient = 0f
+                LossCoefficient = lossCoefficient
             };
             OrificeCount++;
         }
 
         // ---------- Add orifice with inertance ----------
-        // effectiveLength – length of the inertial slug (m), typically the physical neck length.
-        // lossCoefficient – linear resistance (N·s/m⁵). If 0 (or omitted) an automatic stable
-        //                   value will be computed from the pipe's characteristic impedance.
         public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
                                             int areaIndex, float dischargeCoeff,
                                             float effectiveLength, float lossCoefficient = 0f)
         {
-            AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff);
+            // Reuse the base AddOrifice and then override fields
+            AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient);
             ref var d = ref _orifices[OrificeCount - 1];
             d.UseInertance = true;
             d.EffectiveLength = effectiveLength;
-            d.LossCoefficient = lossCoefficient;
+            d.LossCoefficient = lossCoefficient;   // store the linear resistance
         }
 
         public void AddOpenEnd(int pipeIndex, bool isLeftEnd,
@@ -150,7 +146,7 @@ namespace FluidSim.Core
                     ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
                     : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
 
-                // ---- Handle closed orifice as a wall ----
+                // ---- Handle closed orifice (area ≈ 0) as a wall ----
                 if (area < 1e-12f || d.VolumePort == null)
                 {
                     var (rInt, uInt, pInt) = d.IsLeftEnd
@@ -169,7 +165,7 @@ namespace FluidSim.Core
                     continue;
                 }
 
-                // ---- Preliminary isentropic solution (for reference) ----
+                // ---- Preliminary isentropic solution ----
                 float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
                 if (volP >= pipeP)
                 {
@@ -183,31 +179,20 @@ namespace FluidSim.Core
                     mdotEst = -mdotEst;
                 }
 
-                // ---- Compute ghost state ----
+                // ---- Compute final mass flow with limiters ----
                 float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
 
                 if (d.UseInertance)
                 {
-                    // ---- Inertance ODE with (possibly automatic) linear loss ----
                     float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
                     float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
                     float dp = volP - pipeP;
-
-                    // If loss coefficient not provided, use a tiny fraction of the pipe's characteristic impedance
                     float Rlin = d.LossCoefficient;
-                    if (Rlin <= 0f)
-                    {
-                        // Auto‑sized linear drag: 0.5 % of Z_char
-                        float rhoRef = d.CurrentMdot >= 0 ? volRho : pipeRho;
-                        float cRef = d.CurrentMdot >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
-                        float Z_char = rhoRef * cRef / MathF.Max(area, 1e-12f);
-                        Rlin = 0.005f * Z_char;
-                    }
 
                     float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
                     float mdotNew = d.CurrentMdot + dmdot_dt * dt;
 
-                    // Symmetric flow limiters
+                    // Limit outflow from volume (if volume owner is Volume0D)
                     if (d.VolumePort.Owner is Volume0D vol0)
                     {
                         float maxOut = vol0.Mass / dt;
@@ -215,15 +200,19 @@ namespace FluidSim.Core
                         if (mdotNew < -maxOut) mdotNew = -maxOut;
                     }
 
-                    int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
+                    // Limit inflow from pipe – pipe cell cannot be emptied in one step
-                                              : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
+                    {
-                    float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
+                        int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
-                    float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
+                                                : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
-                    float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
+                        float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
-                    float maxFromPipe = pipeCellMass / dt;
+                        float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);  // true cell area, not orifice
-                    if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
+                        float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
+                        float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
+                        float maxFromPipe = pipeCellMass / dt;
+                        if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
+                    }
 
-                    // Velocity clamp 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);
@@ -238,51 +227,60 @@ namespace FluidSim.Core
 
                     d.CurrentMdot = mdotNew;
                     mdotFinal = mdotNew;
-
                     rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
                     pFace = pFaceEst;
                     uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
                 }
                 else
                 {
-                    // ---- Standard quasi‑steady orifice (purely isentropic) ----
+                    // 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;
                     }
 
-                    // Safety velocity clamp (Mach 0.9)
+                    // ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
-                    float cLocal = mdotFinal >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
+                    if (mdotFinal < 0)
-                    float maxULocal = 0.9f * cLocal;
-                    float uCheck = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
-                    if (uCheck > maxULocal)
                     {
-                        uFace = maxULocal;
+                        int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
-                        mdotFinal = rhoFace * uFace * area * (mdotFinal >= 0 ? 1f : -1f);
+                                                : _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;
@@ -299,12 +297,12 @@ namespace FluidSim.Core
                 {
                     d.VolumePort.MassFlowRate = -mdotFinal;
 
-                    if (-mdotFinal >= 0)   // mass flowing into the volume
+                    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
+                    else   // mass leaving volume (into pipe)
                     {
                         d.VolumePort.SpecificEnthalpy = volH;
                     }
@@ -331,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);
@@ -340,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;
@@ -354,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;
+                d.LastMassFlowRate = mdotRaw;
-                float mdot = rhoGhost * uGhost * area;
-                if (d.IsLeftEnd) mdot = -mdot;
-                d.LastMassFlowRate = mdot;
                 d.LastFacePressure = pGhost;
             }
         }

Core/Pipesystem.cs

@@ -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
+        // Precomputed flag: true if a face is a pipe boundary (start or end)
-        private readonly bool[] _isWallFace;
+        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
+            // Mark faces that coincide with a pipe boundary (start or end)
-            _isWallFace = new bool[faceCount];
+            _isPipeBoundaryFace = new bool[faceCount];
-            for (int f = 1; f < _totalCells; f++)
+            for (int p = 0; p < _pipeCount; p++)
             {
-                for (int p = 0; p < _pipeCount; p++)
+                _isPipeBoundaryFace[_pipeStart[p]] = true;
-                {
+                _isPipeBoundaryFace[_pipeEnd[p]] = true;
-                    if (f == _pipeEnd[p] && f < _totalCells)
-                    {
-                        _isWallFace[f] = true;
-                        break;
-                    }
-                }
             }
 
             // 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
+                // Skip faces that belong to a pipe boundary (they are already handled)
-                if (_isWallFace[face])
+                if (_isPipeBoundaryFace[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;
                     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)
+                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
-                        pL.CopyTo(_p, face - 1);
-
-                        // Lax‑Friedrichs fluxes
-                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass
                         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;
+                    int iL = face - 1, iR = face;
-                    float rLf = _rho[iLf], rhouLf = _rhou[iLf];
+                    float rL = _rho[iL], rhouL = _rhou[iL];
-                    float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f));
+                    float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
-                    float uLf = rhouLf * invRhoLf;
+                    float uL = rhouL * invRhoL;
-                    float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf);
+                    float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
-                    float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf);
+                    float cL = MathF.Sqrt(Gamma * pL * invRhoL);
-                    float YLf = _Y[iLf];
+                    float YL = _Y[iL];
-                    _p[iLf] = pLf;
 
-                    float rRf = _rho[iRf], rhouRf = _rhou[iRf];
+                    float rR = _rho[iR], rhouR = _rhou[iR];
-                    float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f));
+                    float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
-                    float uRf = rhouRf * invRhoRf;
+                    float uR = rhouR * invRhoR;
-                    float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf);
+                    float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
-                    float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf);
+                    float cR = MathF.Sqrt(Gamma * pR * invRhoR);
-                    float YRf = _Y[iRf];
+                    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);
+                    float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
-                    ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy);
+                    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;
+            // Update each pipe separately
-            for (; iCell <= _totalCells - vecSize; iCell += vecSize)
+            for (int p = 0; p < _pipeCount; p++)
             {
-                var rhoOld = new Vector<float>(_rho, iCell);
+                int start = _pipeStart[p];
-                var rhouOld = new Vector<float>(_rhou, iCell);
+                int end = _pipeEnd[p];    // exclusive
-                var EOld = new Vector<float>(_E, iCell);
+                int len = end - start;
-                var YOld = new Vector<float>(_Y, iCell);
+                if (len == 0) continue;
 
-                var fluxM_L = new Vector<float>(_fluxM, iCell);
+                // ------- Left boundary cell (i = start) ------
-                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);
+                    int i = start;
-                    var relax = new Vector<float>(_relaxFactors, iCell);
+                    float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
-                    ENew = ambRef + (ENew - ambRef) * relax;
+
+                    // 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));
+                // ------- Interior cells (i = start+1 to end-2) ------
-                var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
+                if (len > 2)
-                var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
+                {
-                ENew = Vector.Max(ENew, eMin);
+                    int iCell = start + 1;
+                    int iEnd = end - 1; // exclusive upper bound
 
-                rhoNew.CopyTo(_rho, iCell);
+                    // Vectorised path for interior cells (if available)
-                rhouNew.CopyTo(_rhou, iCell);
+                    for (; iCell <= iEnd - vecSize; iCell += vecSize)
-                ENew.CopyTo(_E, iCell);
+                    {
-                var yNew = rhoYNew / rhoNew;
+                        var rhoOld = new Vector<float>(_rho, iCell);
-                yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
+                        var rhouOld = new Vector<float>(_rhou, iCell);
-                yNew.CopyTo(_Y, 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)
+            // Recompute pressure for all cells (for visualization)
-            for (; iCell < _totalCells; iCell++)
+            for (int i = 0; i < _totalCells; i++)
             {
-                float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
+                float rho = _rho[i];
-                float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
+                float rhou = _rhou[i];
-                float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
+                float u = rhou / MathF.Max(rho, 1e-12f);
-                float dtdx = dt / _dx[iCell];
+                _p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
-
-                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);
             }
         }
 
-        // ---------- 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()
         {

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++)

Scenarios/SingleCylScenario.cs

@@ -1,6 +1,7 @@
 using FluidSim.Components;
 using FluidSim.Core;
 using FluidSim.Interfaces;
+using FluidSim.Utils;
 using SFML.Graphics;
 using SFML.System;
 using System;
@@ -9,208 +10,190 @@ namespace FluidSim.Tests
 {
     public class SingleCylScenario : Scenario
     {
-        // ---------- Engine components ----------
         private Crankshaft crankshaft;
         private Cylinder cylinder;
 
-        // ---------- Fluid network ----------
         private PipeSystem pipeSystem;
         private BoundarySystem boundaries;
         private Solver solver;
 
-        // Volumes
         private Volume0D intakePlenum;
-
-        // Ports
         private Port plenumInlet, plenumOutlet;
+        private Volume0D exhaustCollector;
+        private Port colIn, colOut;
 
-        // Orifice / open‑end indices
+        private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx;
-        private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
-        private int intakeOpenIdx, exhaustOpenIdx;
         private float[] orificeAreas;
+        private int intakeOpenIdx, exhaustOpenIdx;
 
-        // Sound
         private SoundProcessor exhaustSound, intakeSound;
         private OutdoorExhaustReverb reverb;
 
-        // ---------- Simulation state ----------
         private double dt;
         private int stepCount;
 
-        // ---------- Geometry (Lifan YX140) ----------
+        // Use a private field for the maximum throttle area, avoiding any base‑class conflicts
-        // Bore 56 mm, Stroke 57 mm, CR 9.5
+        private float _maxThrottleArea;
-        private const float Bore = 0.056f;
-        private const float Stroke = 0.057f;
-        private const float ConRod = 0.110f;           // typical for 57 mm stroke
-        private const float CompressionRatio = 9.5f;
 
-        // Valve diameters (intake 27 mm, exhaust 23 mm)
+        // pipe area for open end calculations
-        private const float IntakeValveDiam = 0.027f;
+        private float pipeArea;
-        private const float ExhaustValveDiam = 0.023f;
-        private const float ValveLift = 0.006f;         // 6 mm peak lift
-
-        // Valve timings (degrees, 720° four‑stroke)
-        // Intake: 15° BTDC → 45° ABDC
-        private const float IVO = 345f;   // 15° BTDC
-        private const float IVC = 585f;   // 45° ABDC (180°+45°)
-        // Exhaust: 45° BBDC → 15° ATDC
-        private const float EVO = 135f;   // 45° BBDC (180°-45°)
-        private const float EVC = 375f;   // 15° ATDC (360°+15°)
-
-        // Spark advance: 30° BTDC
-        private const float SparkAdv = 30f;
-
-        // Pipe / plenum sizes
-        private const float PipeDiam = 0.025f;          // 25 mm intake / exhaust
-        private const float PipeArea = 0.00049087f;      // π*D²/4
-        private const float PlenumVolume = 0.0005f;      // 500 mL
-        private const float MaxThrottleArea = 1e-4f;     // ~1 cm² (fully open)
-
-        // Pipe lengths and cell counts
-        private const float IntakeLenBefore = 0.15f;     // 15 cm before throttle
-        private const float RunnerLen = 0.25f;           // 25 cm runner
-        private const float ExhaustLen = 0.60f;          // 60 cm exhaust
-        private const int CellsBefore = 6;
-        private const int CellsRunner = 10;
-        private const int CellsExhaust = 24;
 
         public override void Initialize(int sampleRate)
         {
             dt = 1.0 / sampleRate;
 
-            // ---------- Crankshaft ----------
+            // Maximum throttle area – independent of base class
-            crankshaft = new Crankshaft(600);            // start at ~600 RPM
+            _maxThrottleArea = (float)Units.AreaFromDiameter(3 * Units.cm);   // 1 cm²
-            crankshaft.Inertia = 0.2f;
-            crankshaft.FrictionConstant = 2.0f;
-            crankshaft.FrictionViscous = 0.04f;
 
-            // ---------- Cylinder ----------
+            // ---- Crankshaft ----
-            cylinder = new Cylinder(Bore, Stroke, ConRod, CompressionRatio,
+            crankshaft = new Crankshaft(2000);
-                                    IVO, IVC, EVO, EVC, crankshaft)
+            crankshaft.Inertia = 0.01f;
+            crankshaft.FrictionConstant = 2f;
+            crankshaft.FrictionViscous = 0.0f;
+
+            // ---- Cylinder ----
+            float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 11f;
+            float ivo = 350f, ivc = 580f, evo = 120f, evc = 370f;
+            cylinder = new Cylinder(bore, stroke, conRod, compRatio,
+                                    ivo, ivc, evo, evc, crankshaft)
             {
-                IntakeValveDiameter = IntakeValveDiam,
+                IntakeValveDiameter = 0.03f,
-                ExhaustValveDiameter = ExhaustValveDiam,
+                IntakeValveLift = 0.005f,
-                IntakeValveLift = ValveLift,
+                ExhaustValveDiameter = 0.028f,
-                ExhaustValveLift = ValveLift,
+                ExhaustValveLift = 0.005f
-                SparkAdvance = SparkAdv,
-                EnergyVariationFraction = 0.03f,   // small cycle‑to‑cycle variation
-                MisfireProbability = 0.0f
             };
 
-            // ---------- Pipe system ----------
+            // ---- Pipe system ----
-            int totalCells = CellsBefore + CellsRunner + CellsExhaust;
+            int[] pipeStart = { 0, 10, 20 };
-            int[] pipeStart = { 0, CellsBefore, CellsBefore + CellsRunner };
+            int[] pipeEnd   = { 10, 20, 70 };
-            int[] pipeEnd   = { CellsBefore, CellsBefore + CellsRunner, totalCells };
+            int totalCells = pipeEnd[^1];   // automatically 70, stays in sync
-
+            float[] area = new float[totalCells];
-            float[] areas = new float[totalCells];
+            float[] dx = new float[totalCells];
-            float[] dxs   = new float[totalCells];
+            float pipeDiameter = 0.02f;          // 2 cm
-            float dxBefore = IntakeLenBefore / CellsBefore;
+            pipeArea = MathF.PI * 0.25f * pipeDiameter * pipeDiameter;
-            float dxRunner = RunnerLen / CellsRunner;
+            float areaVal = pipeArea;
-            float dxExh    = ExhaustLen / CellsExhaust;
+            float intakeLenBefore = 0.2f, intakeLenRunner = 0.2f, exhaustLen = 0.4f;
-
             for (int i = 0; i < totalCells; i++)
             {
-                areas[i] = PipeArea;
+                area[i] = areaVal;
-                if (i < CellsBefore)
+                if (i < 10) dx[i] = intakeLenBefore / 10f;
-                    dxs[i] = dxBefore;
+                else if (i < 20) dx[i] = intakeLenRunner / 10f;
-                else if (i < CellsBefore + CellsRunner)
+                else dx[i] = exhaustLen / 50f;
-                    dxs[i] = dxRunner;
-                else
-                    dxs[i] = dxExh;
             }
 
-            float rho0 = 101325f / (287f * 300f);
+            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
-            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, areas, dxs,
+                                        1.225f, 0f, 101325f);
-                                        rho0, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 1.0f;
-            pipeSystem.DampingMultiplier = 0.5f;
+            pipeSystem.EnergyRelaxationRate = 0.5f;
-            pipeSystem.EnergyRelaxationRate = 0f;        // adiabatic pipes
             pipeSystem.AmbientPressure = 101325f;
 
-            // ---------- Volumes ----------
+            // ---- Volumes ----
-            intakePlenum = new Volume0D(PlenumVolume, 101325f, 300f);
+            intakePlenum = new Volume0D(100e-6f, 101325f, 300f);  // 100 mL
-            plenumInlet  = intakePlenum.CreatePort();
+            plenumInlet = intakePlenum.CreatePort();
             plenumOutlet = intakePlenum.CreatePort();
+            exhaustCollector = new Volume0D(10e-6f, 101325f, 800f);  // 10 mL (unused but present)
+            colIn = exhaustCollector.CreatePort();
+            colOut = exhaustCollector.CreatePort();
 
-            // ---------- Boundary system ----------
+            // ---- Boundary system ----
             boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
 
-            throttleAreaIdx  = 0;
+            throttleAreaIdx = 0;
-            plenumRunnerIdx  = 1;
+            plenumRunnerAreaIdx = 1;
-            intakeValveIdx   = 2;
+            intakeValveIdx = 2;
-            exhaustValveIdx  = 3;
+            exhaustValveIdx = 3;
 
-            // Open ends
+            // Intake open end (pipe0 left)
-            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, PipeArea);
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, pipeArea);
             intakeOpenIdx = 0;
-            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, PipeArea);
+
+            // Throttle orifice (plenum inlet to pipe0 right)
+            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.2f);
+
+            // Plenum to runner (plenum outlet to pipe1 left)
+            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1f);
+
+            // Intake valve (cylinder intake to pipe1 right)
+            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1f);
+
+            // Exhaust valve (cylinder exhaust to pipe2 left)
+            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1f);
+
+            // Exhaust open end (pipe2 right)
+            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, pipeArea);
             exhaustOpenIdx = 1;
 
-            // Orifices
-            // throttle – variable area, low discharge for restriction
-            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false,
-                                  throttleAreaIdx, dischargeCoeff: 0.8f);
-            // plenum → runner
-            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true,
-                                  plenumRunnerIdx, dischargeCoeff: 1.0f);
-            // intake valve
-            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false,
-                                  intakeValveIdx, dischargeCoeff: 1.0f);
-            // exhaust valve
-            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true,
-                                  exhaustValveIdx, dischargeCoeff: 1.0f);
-
             orificeAreas = new float[4];
-            orificeAreas[plenumRunnerIdx] = PipeArea;   // fixed full‑bore
+            orificeAreas[plenumRunnerAreaIdx] = areaVal; // fixed plenum->runner area
 
-            // ---------- Solver ----------
+            // ---- Solver ----
-            solver = new Solver { SubStepCount = 5, EnableProfiling = false };
+            solver = new Solver { SubStepCount = 4, EnableProfiling = false };
             solver.SetTimeStep(dt);
             solver.SetPipeSystem(pipeSystem);
             solver.SetBoundarySystem(boundaries);
             solver.AddComponent(cylinder);
             solver.AddComponent(intakePlenum);
+            solver.AddComponent(exhaustCollector);
 
-            // ---------- Sound ----------
+            // ---- Sound ----
-            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
+            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
-            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
+            intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
             reverb = new OutdoorExhaustReverb(sampleRate);
 
             stepCount = 0;
-            Console.WriteLine("Single‑cylinder engine (YX140) ready.");
+            Console.WriteLine("TestScenario ready.");
         }
 
         public override float Process()
         {
-            // ---- Crank and cylinder pre‑step ----
             crankshaft.Step((float)dt);
             cylinder.PreStep((float)dt);
 
-            // ---- Update variable areas ----
+            // Update variable orifice areas – use the private _maxThrottleArea
-            float throttledArea = MaxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1.0f);
+            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1f);
             orificeAreas[throttleAreaIdx] = throttledArea;
-            orificeAreas[intakeValveIdx]  = cylinder.IntakeValveArea;
+            orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
             orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
             boundaries.SetOrificeAreas(orificeAreas);
 
-            // ---- Fluids step ----
             solver.Step();
             stepCount++;
 
-            // ---- Sound ----
+            // Retrieve open‑end mass flows for sound synthesis
             float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
-            float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
+            float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
 
             float exhaustDry = exhaustSound.Process(exhaustFlow);
-            float intakeDry  = intakeSound.Process(intakeFlow);
+            float intakeDry = intakeSound.Process(intakeFlow);
 
-            if (stepCount % 2000 == 0)
+            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, " +
+                float crankDeg = crankshaft.CrankAngle;   // degrees (0–720)
-                                  $"Throttle={Throttle * 100:F0}%");
+                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(exhaustDry + intakeDry);
+            return reverb.Process(intakeDry + exhaustDry);
         }
 
         public override void Draw(RenderWindow target)
@@ -220,53 +203,44 @@ namespace FluidSim.Tests
 
             float intakeY = winH / 2f - 40f;
             float exhaustY = winH / 2f + 80f;
-            float leftX = 40f;
+            float openEndX = 40f;
 
-            // Intake open end marker
+            // Intake pipe before throttle (pipe 0)
-            var om = new CircleShape(5f) { FillColor = Color.Cyan };
+            float pipe1StartX = openEndX;
-            om.Position = new Vector2f(leftX - 5f, intakeY - 5f);
+            float pipe1EndX = pipe1StartX + 120f;
-            target.Draw(om);
+            DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
-
-            // Pipe 0 – before throttle
-            float p0EndX = leftX + 80f;
-            DrawPipe(target, pipeSystem, 0, intakeY, leftX, p0EndX);
 
             // Throttle symbol
-            float thrX = p0EndX + 5f;
+            float throttleX = pipe1EndX + 5f;
-            var thr = new RectangleShape(new Vector2f(8f, 30f))
+            var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
             {
                 FillColor = Color.Yellow,
-                Position = new Vector2f(thrX, intakeY - 15f)
+                Position = new Vector2f(throttleX, intakeY - 15f)
             };
-            target.Draw(thr);
+            target.Draw(throttleRect);
 
-            // Plenum volume
+            // Plenum
-            float plenW = 60f, plenH = 50f;
+            float plenW = 60f, plenH = 80f;
-            float plenLeftX = thrX + 12f;
+            float plenLeftX = throttleX + 10f;
             float plenCenterX = plenLeftX + plenW / 2f;
             float plenTopY = intakeY - plenH / 2f;
             DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
 
-            // Pipe 1 – runner
+            // Runner pipe (pipe 1)
-            float rStartX = plenLeftX + plenW + 10f;
+            float runnerStartX = plenLeftX + plenW + 5f;
-            float rEndX = rStartX + 100f;
+            float runnerEndX = runnerStartX + 100f;
-            DrawPipe(target, pipeSystem, 1, intakeY, rStartX, rEndX);
+            DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
 
             // Cylinder
-            float cylCX = rEndX + 50f;
+            float cylCX = runnerEndX + 50f;
             float cylTopY = intakeY - 120f;
             float cylW = 80f, cylMaxH = 240f;
             DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
 
-            // Pipe 2 – exhaust
+            // Exhaust pipe (pipe 2)
             float exhStartX = cylCX + cylW / 2f + 20f;
             float exhEndX = winW - 60f;
             DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
-
-            // Exhaust open end
-            var em = new CircleShape(5f) { FillColor = Color.Magenta };
-            em.Position = new Vector2f(exhEndX - 5f, exhaustY - 5f);
-            target.Draw(em);
         }
     }
 }