"better" two stroke engine

Components/Crankshaft.cs

@@ -4,25 +4,51 @@ namespace FluidSim.Components
 {
     public class Crankshaft
     {
-        public float AngularVelocity;   // rad/s
-        public float CrankAngle;        // rad, 0 … 4π
+        public float AngularVelocity;
+        public float CrankAngle;
         public float PreviousAngle;
 
         public float Inertia = 0.2f;
-        public float FrictionConstant;   // N·m
-        public float FrictionViscous;    // N·m per rad/s
+        public float FrictionConstant;
+        public float FrictionViscous;
+        public float LastNetTorque { get; private set; }
+        public float AveragePower { get; private set; }
+        public float AverageTorque { get; private set; }
 
         private float externalTorque;
+        private float _loadTorque;
+
+        private readonly float[] _powerBuffer;
+        private int _powerBufIdx, _powerBufCount;
+        private float _powerBufSum;
+        private readonly float[] _torqueBuffer;
+        private int _torqueBufIdx, _torqueBufCount;
+        private float _torqueBufSum;
+
+        /// <summary>Engine cycle length in radians. 4π = four‑stroke, 2π = two‑stroke.</summary>
+        public float CycleLength { get; set; } = 4f * MathF.PI;
 
         public Crankshaft(float initialRPM = 400f)
         {
             AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
             CrankAngle = 0f;
             PreviousAngle = 0f;
+
+            _powerBuffer = new float[16384];
+            _torqueBuffer = new float[16384];
         }
 
         public void AddTorque(float torque) => externalTorque += torque;
 
+        public void SetLoadTorque(float torque) => _loadTorque = Math.Max(torque, 0f);
+
+        private float _effectiveInertia;   // if >0, overrides Inertia
+
+        public void SetEffectiveInertia(float inertia)
+        {
+            _effectiveInertia = inertia;
+        }
+
         public void Step(float dt)
         {
             if (float.IsNaN(AngularVelocity) || float.IsInfinity(AngularVelocity))
@@ -34,17 +60,42 @@ namespace FluidSim.Components
 
             float friction = FrictionConstant * MathF.Sign(AngularVelocity)
                            + FrictionViscous * AngularVelocity;
-            float netTorque = externalTorque - friction;
-            float alpha = netTorque / Inertia;
-            AngularVelocity += alpha * dt;
 
+            float netTorque = externalTorque - friction;
+            LastNetTorque = netTorque;
+
+            float totalNetTorque = netTorque - _loadTorque;
+            float currentInertia = _effectiveInertia > 0f ? _effectiveInertia : Inertia;
+            float alpha = totalNetTorque / currentInertia;
+            AngularVelocity += alpha * dt;
             if (AngularVelocity < 0f) AngularVelocity = 0f;
 
             CrankAngle += AngularVelocity * dt;
-            if (CrankAngle >= 4f * MathF.PI)
-                CrankAngle -= 4f * MathF.PI;
+            if (CrankAngle >= CycleLength)
+                CrankAngle -= CycleLength;
             else if (CrankAngle < 0f)
-                CrankAngle += 4f * MathF.PI;
+                CrankAngle += CycleLength;
+
+            // Power averaging
+            float instantPower = netTorque * AngularVelocity;
+            if (_powerBufCount == _powerBuffer.Length)
+                _powerBufSum -= _powerBuffer[_powerBufIdx];
+            else
+                _powerBufCount++;
+            _powerBuffer[_powerBufIdx] = instantPower;
+            _powerBufSum += instantPower;
+            _powerBufIdx = (_powerBufIdx + 1) % _powerBuffer.Length;
+            AveragePower = _powerBufSum / _powerBufCount;
+
+            // Torque averaging
+            if (_torqueBufCount == _torqueBuffer.Length)
+                _torqueBufSum -= _torqueBuffer[_torqueBufIdx];
+            else
+                _torqueBufCount++;
+            _torqueBuffer[_torqueBufIdx] = netTorque;
+            _torqueBufSum += netTorque;
+            _torqueBufIdx = (_torqueBufIdx + 1) % _torqueBuffer.Length;
+            AverageTorque = _torqueBufSum / _torqueBufCount;
 
             externalTorque = 0f;
         }

Components/Cylinder.cs

@@ -1,99 +1,25 @@
 using System;
-using System.Collections.Generic;
-using FluidSim.Interfaces;
+using FluidSim.Components;   // if needed
 
 namespace FluidSim.Components
 {
-    public class Cylinder : IComponent
+    public class Cylinder : EngineCylinder
     {
-        public Port IntakePort { get; }
-        public Port ExhaustPort { get; }
-        public Crankshaft Crankshaft { get; }
+        public float IVO, IVC, EVO, EVC;   // degrees in a 720° cycle
 
-        private readonly Port[] _ports;
-        IReadOnlyList<Port> IComponent.Ports => _ports;
+        protected override float CycleLengthRad => 4f * MathF.PI;
+        protected override float MaxCycleDeg => 720f;
 
-        public float Bore { get; }
-        public float Stroke { get; }
-        public float ConRodLength { get; }
-        public float CompressionRatio { get; }
-
-        public float IVO, IVC, EVO, EVC;     // degrees
-        public float IntakeValveDiameter = 0.03f;
-        public float ExhaustValveDiameter = 0.028f;
-        public float IntakeValveLift = 0.005f;
-        public float ExhaustValveLift = 0.005f;
-
-        public float IntakeValveMaxArea => MathF.PI * IntakeValveDiameter * IntakeValveLift;
-        public float ExhaustValveMaxArea => MathF.PI * ExhaustValveDiameter * ExhaustValveLift;
-
-        public float SparkAdvance = 20f;
-        public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f;
-        public float StoichiometricAFR = 14.7f;
-        public float FuelLowerHeatingValue = 44e6f;
-        public float EnergyVariationFraction = 0.05f;
-        public float MisfireProbability = 0.01f;
-        public float CylinderWallArea = 0.02f;
-        public float HeatTransferCoefficient = 100f;
-        public float AmbientTemperature = 300f;
-
-        public float PhaseOffset;   // rad
-
-        public float Volume => cylinderVolume;
-        public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
-        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
-        public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
-        public float Mass => _airMass + _exhaustMass;
-        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
-        public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
-
-        private float cylinderVolume, cylinderEnergy;
-        private float _airMass, _exhaustMass;
-        private float trappedAirMass, fuelMass, burnFraction;
-        private bool combustionActive, fuelInjected;
-        private float _energyFactor = 1f;
-        private readonly Random _random = new Random();
-
-        private const float Gamma = 1.4f;
-        private const float GasConstant = 287f;
-        private const float MaxPressurePa = 200e5f;
-        private const float MaxTemperatureK = 3500f;
+        public override float IntakeValveArea =>
+            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
+        public override float ExhaustValveArea =>
+            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
 
         public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio,
                         float ivo, float ivc, float evo, float evc, Crankshaft crankshaft)
+            : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
         {
-            Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
-            CompressionRatio = compressionRatio;
             IVO = ivo; IVC = ivc; EVO = evo; EVC = evc;
-            Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
-
-            cylinderVolume = clearanceVolume;
-            float initRho = 1.225f;
-            _airMass = initRho * clearanceVolume;
-            _exhaustMass = 0f;
-            cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
-
-            IntakePort = new Port { Owner = this };
-            ExhaustPort = new Port { Owner = this };
-            _ports = new[] { IntakePort, ExhaustPort };
-        }
-
-        private float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
-        private float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
-        private float CrankRadius => Stroke * 0.5f;
-        private float Obliquity => CrankRadius / ConRodLength;
-
-        private float CrankDeg =>
-            ((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f;
-
-        public float ComputeVolume(float thetaRad)
-        {
-            float r = CrankRadius, l = ConRodLength;
-            float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
-            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
-            float x = r * (1f - cosTh) + l * (1f - term);
-            float area = MathF.PI * 0.25f * Bore * Bore;
-            return clearanceVolume + area * x;
         }
 
         private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
@@ -101,19 +27,12 @@ namespace FluidSim.Components
             float deg = thetaDeg % 720f;
             if (deg < 0f) deg += 720f;
 
-            float duration;
             float effectiveOpen = opens;
             float effectiveClose = closes;
-
-            if (closes < opens)
-            {
-                // Wrap‑around case (e.g., exhaust: opens near 480°, closes near 30°)
-                effectiveClose += 720f;
-            }
-            duration = effectiveClose - effectiveOpen;
+            if (closes < opens) effectiveClose += 720f;
+            float duration = effectiveClose - effectiveOpen;
             if (duration <= 0f) return 0f;
 
-            // Map the angle into the [opens, opens+duration] window
             float mapped = deg;
             if (mapped < opens) mapped += 720f;
             if (mapped < opens || mapped > effectiveClose) return 0f;
@@ -138,39 +57,9 @@ namespace FluidSim.Components
             return 0f;
         }
 
-        public float IntakeValveArea =>
-            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-        public float ExhaustValveArea =>
-            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-
-        private float Wiebe(float angleSinceSpark)
+        protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
         {
-            if (angleSinceSpark < WiebeStart) return 0f;
-            float phi = (angleSinceSpark - WiebeStart) / WiebeDuration;
-            if (phi <= 0f) return 0f;
-            return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f));
-        }
-
-        public void PreStep(float dt)
-        {
-            float prevVolume = cylinderVolume;
-            float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
-            cylinderVolume = ComputeVolume(crankAngleRad);
-
-            float dV = cylinderVolume - prevVolume;
-            float pRel = Pressure - 101325f;
-            float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
-            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
-            float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
-            float pistonArea = MathF.PI * 0.25f * Bore * Bore;
-            Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
-
-            cylinderEnergy -= Pressure * dV;
-
-            float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % 720f;
-            float currDeg = crankAngleRad * 180f / MathF.PI % 720f;
-
-            // Intake closing
+            // Intake closing → fuel injection
             if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
             {
                 trappedAirMass = _airMass;
@@ -178,11 +67,14 @@ namespace FluidSim.Components
                 fuelInjected = true;
             }
 
-            // Spark
-            float sparkAngle = 0f - SparkAdvance;
-            if (sparkAngle < 0f) sparkAngle += 720f;
-            bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) ||
-                                (prevDeg > sparkAngle + 360f && currDeg < sparkAngle);
+            // Spark – occurs at TDC (0°) minus advance, every 720°
+            float sparkAngle = (0f - SparkAdvance + 720f) % 720f;
+            bool crossedSpark = false;
+            if (prevDeg < sparkAngle && currDeg >= sparkAngle)
+                crossedSpark = true;
+            else if (prevDeg > sparkAngle && currDeg < sparkAngle)
+                crossedSpark = true;
+
             if (crossedSpark && !combustionActive && fuelInjected)
             {
                 if (_random.NextDouble() < MisfireProbability)
@@ -197,7 +89,7 @@ namespace FluidSim.Components
                 }
             }
 
-            // Combustion
+            // Combustion progression
             if (combustionActive)
             {
                 float angleSinceSpark = currDeg - sparkAngle;
@@ -220,62 +112,6 @@ namespace FluidSim.Components
                     burnFraction = newFraction;
                 }
             }
-
-            // Heat loss
-            float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
-                            (Temperature - AmbientTemperature) * dt;
-            cylinderEnergy -= dQ_loss;
-
-            // Update port states
-            float p = Pressure, rho = Density, T = Temperature;
-            float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
-            float af = AirFraction;
-            IntakePort.Pressure = p; IntakePort.Density = rho;
-            IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
-            ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
-            ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
-        }
-
-        public void UpdateState(float dt)
-        {
-            float dmAir = 0f, dmExhaust = 0f, dE = 0f;
-            foreach (var port in _ports)
-            {
-                float mdot = port.MassFlowRate;
-                float af = mdot >= 0f ? port.AirFraction : AirFraction;
-                dmAir += mdot * af * dt;
-                dmExhaust += mdot * (1f - af) * dt;
-                dE += mdot * port.SpecificEnthalpy * dt;
-            }
-
-            _airMass += dmAir; _exhaustMass += dmExhaust;
-            cylinderEnergy += dE;
-
-            float V = MathF.Max(cylinderVolume, 1e-12f);
-            float currentP = (Gamma - 1f) * cylinderEnergy / V;
-            if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f);
-
-            float currentRho = (_airMass + _exhaustMass) / V;
-            float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f);
-            if (currentT > MaxTemperatureK)
-            {
-                float pAtTlimit = currentRho * GasConstant * MaxTemperatureK;
-                cylinderEnergy = pAtTlimit * V / (Gamma - 1f);
-            }
-
-            float totalMass = _airMass + _exhaustMass;
-            if (totalMass < 1e-9f)
-            {
-                _airMass = 1e-9f; _exhaustMass = 0f;
-                cylinderEnergy = 101325f * V / (Gamma - 1f);
-            }
-            else if (cylinderEnergy < 0f)
-            {
-                cylinderEnergy = 101325f * V / (Gamma - 1f);
-            }
-
-            if (_airMass < 0f) _airMass = 0f;
-            if (_exhaustMass < 0f) _exhaustMass = 0f;
         }
     }
 }

Components/EngineCylinder.cs

@@ -0,0 +1,203 @@
+using System;
+using System.Collections.Generic;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    /// <summary>Common base for all reciprocating engine cylinders.</summary>
+    public abstract class EngineCylinder : IComponent
+    {
+        public Port IntakePort { get; }
+        public Port ExhaustPort { get; }
+        public Crankshaft Crankshaft { get; }
+
+        private readonly Port[] _ports;
+        IReadOnlyList<Port> IComponent.Ports => _ports;
+
+        // ----- Geometry -----
+        public float Bore { get; }
+        public float Stroke { get; }
+        public float ConRodLength { get; }
+        public float CompressionRatio { get; }
+
+        // ----- Valve / port sizes (used for curtain area) -----
+        public float IntakeValveDiameter = 0.03f;
+        public float ExhaustValveDiameter = 0.028f;
+        public float IntakeValveLift = 0.005f;
+        public float ExhaustValveLift = 0.005f;
+
+        // ----- Combustion -----
+        public float SparkAdvance = 20f;
+        public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f;
+        public float StoichiometricAFR = 14.7f;
+        public float FuelLowerHeatingValue = 44e6f;
+        public float EnergyVariationFraction = 0.05f;
+        public float MisfireProbability = 0f;
+        public float CylinderWallArea = 0.02f;
+        public float HeatTransferCoefficient = 100f;
+        public float AmbientTemperature = 300f;
+
+        public float PhaseOffset;   // radians
+
+        // ----- State (public, used by drawing) -----
+        public float Volume => cylinderVolume;
+        public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
+        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
+        public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
+        public float Mass => _airMass + _exhaustMass;
+        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
+        public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
+
+        protected float cylinderVolume, cylinderEnergy;
+        protected float _airMass, _exhaustMass;
+        protected float trappedAirMass, fuelMass, burnFraction;
+        protected bool combustionActive, fuelInjected;
+        protected float _energyFactor = 1f;
+        protected readonly Random _random = new Random();
+
+        protected const float Gamma = 1.4f;
+        protected const float GasConstant = 287f;
+        protected const float MaxPressurePa = 200e5f;
+        protected const float MaxTemperatureK = 3500f;
+
+        // ----- Derived geometry (cycle‑independent) -----
+        protected float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
+        protected float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
+        protected float CrankRadius => Stroke * 0.5f;
+        protected float Obliquity => CrankRadius / ConRodLength;
+
+        // ----- Abstract members (cycle‑specific) -----
+        protected abstract float CycleLengthRad { get; }      // 4π or 2π
+        protected abstract float MaxCycleDeg { get; }         // 720 or 360
+        public abstract float IntakeValveArea { get; }
+        public abstract float ExhaustValveArea { get; }
+        protected abstract void HandleCycleEvents(float prevDeg, float currDeg, float dt);
+
+        protected EngineCylinder(float bore, float stroke, float conRodLength,
+                                 float compressionRatio, Crankshaft crankshaft)
+        {
+            Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
+            CompressionRatio = compressionRatio;
+            Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
+
+            cylinderVolume = clearanceVolume;
+            float initRho = 1.225f;
+            _airMass = initRho * clearanceVolume;
+            _exhaustMass = 0f;
+            cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
+
+            IntakePort = new Port { Owner = this };
+            ExhaustPort = new Port { Owner = this };
+            _ports = new[] { IntakePort, ExhaustPort };
+
+            // Set crankshaft cycle length
+            crankshaft.CycleLength = CycleLengthRad;
+        }
+
+        public float ComputeVolume(float thetaRad)
+        {
+            float r = CrankRadius, l = ConRodLength;
+            float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
+            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
+            float x = r * (1f - cosTh) + l * (1f - term);
+            float area = MathF.PI * 0.25f * Bore * Bore;
+            return clearanceVolume + area * x;
+        }
+
+        protected float CrankDeg =>
+            ((Crankshaft.CrankAngle + PhaseOffset) % CycleLengthRad) * 180f / MathF.PI;
+
+        protected float Wiebe(float angleSinceSpark)
+        {
+            if (angleSinceSpark < WiebeStart) return 0f;
+            float phi = (angleSinceSpark - WiebeStart) / WiebeDuration;
+            return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f));
+        }
+
+        // ----- Main update called before flow solver -----
+        public void PreStep(float dt)
+        {
+            // Speed‑dependent spark advance
+            float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f);
+
+            float prevVolume = cylinderVolume;
+            float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
+            cylinderVolume = ComputeVolume(crankAngleRad);
+
+            // Piston work
+            float dV = cylinderVolume - prevVolume;
+            float pRel = Pressure - 101325f;
+            float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
+            float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
+            float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
+            float pistonArea = MathF.PI * 0.25f * Bore * Bore;
+            Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
+
+            cylinderEnergy -= Pressure * dV;
+
+            float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % MaxCycleDeg;
+            float currDeg = crankAngleRad * 180f / MathF.PI % MaxCycleDeg;
+
+            // Let derived class handle valve events, spark, fuel
+            HandleCycleEvents(prevDeg, currDeg, dt);
+
+            // Heat loss
+            float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
+                            (Temperature - AmbientTemperature) * dt;
+            cylinderEnergy -= dQ_loss;
+
+            // Update port states
+            float p = Pressure, rho = Density, T = Temperature;
+            float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
+            float af = AirFraction;
+            IntakePort.Pressure = p; IntakePort.Density = rho;
+            IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
+            ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
+            ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
+        }
+
+        // ----- State update (mass/energy balance) -----
+        public void UpdateState(float dt)
+        {
+            float dmAir = 0f, dmExhaust = 0f, dE = 0f;
+            foreach (var port in _ports)
+            {
+                float mdot = port.MassFlowRate;
+                float af = mdot >= 0f ? port.AirFraction : AirFraction;
+                dmAir += mdot * af * dt;
+                dmExhaust += mdot * (1f - af) * dt;
+                dE += mdot * port.SpecificEnthalpy * dt;
+            }
+
+            _airMass += dmAir; _exhaustMass += dmExhaust;
+            cylinderEnergy += dE;
+
+            float V = MathF.Max(cylinderVolume, 1e-12f);
+            float currentP = (Gamma - 1f) * cylinderEnergy / V;
+            if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f);
+
+            float currentRho = (_airMass + _exhaustMass) / V;
+            float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f);
+            if (currentT > MaxTemperatureK)
+            {
+                float pAtTlimit = currentRho * GasConstant * MaxTemperatureK;
+                cylinderEnergy = pAtTlimit * V / (Gamma - 1f);
+            }
+
+            float totalMass = _airMass + _exhaustMass;
+            if (totalMass < 1e-9f)
+            {
+                _airMass = 1e-9f; _exhaustMass = 0f;
+                cylinderEnergy = 101325f * V / (Gamma - 1f);
+            }
+            else if (cylinderEnergy < 0f)
+            {
+                cylinderEnergy = 101325f * V / (Gamma - 1f);
+            }
+
+            if (_airMass < 0f) _airMass = 0f;
+            if (_exhaustMass < 0f) _exhaustMass = 0f;
+        }
+    }
+}

Components/TwoStrokeCylinder.cs

@@ -0,0 +1,183 @@
+using System;
+
+namespace FluidSim.Components
+{
+    /// <summary>
+    /// Two-stroke cylinder with symmetrical port timings centred on BDC (180°).
+    ///
+    /// Changes vs. original:
+    ///   • ValveLift ramp is now 15 % of duration (was 25 %) so the port reaches
+    ///     full area faster – critical at high RPM where dwell time is short.
+    ///   • Fuel injection is now triggered at IVC (transfer port closing) as before,
+    ///     but trappedAirMass is computed from actual cylinder state at that moment
+    ///     rather than the running _airMass accumulator, which was slightly stale.
+    ///   • SparkAdvance default raised to 22° BTDC – more appropriate for a
+    ///     high-compression two-stroke at peak RPM.  The scenario can still override it.
+    /// </summary>
+    public class TwoStrokeCylinder : EngineCylinder
+    {
+        // ── Port timing read-outs (degrees, 0 = TDC) ───────────────────────────
+        public float IVO => 180f - TransferDuration / 2f;   // transfer opens
+        public float IVC => 180f + TransferDuration / 2f;   // transfer closes
+        public float EVO => 180f - ExhaustDuration  / 2f;   // exhaust opens
+        public float EVC => 180f + ExhaustDuration  / 2f;   // exhaust closes
+
+        // ── Configurable durations ──────────────────────────────────────────────
+        public float TransferDuration { get; }   // default: 155°
+        public float ExhaustDuration  { get; }   // default: 195°
+
+        // Fraction of port-open duration used for ramp-up / ramp-down.
+        // 0.15 → port at full area for the middle 70 % of open time.
+        private const float RampFraction = 0.15f;
+
+        protected override float CycleLengthRad => 2f * MathF.PI;
+        protected override float MaxCycleDeg    => 360f;
+
+        public override float IntakeValveArea =>
+            MathF.PI * IntakeValveDiameter
+                     * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
+
+        public override float ExhaustValveArea =>
+            MathF.PI * ExhaustValveDiameter
+                     * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
+
+        // ── Constructor ─────────────────────────────────────────────────────────
+        public TwoStrokeCylinder(float bore, float stroke, float conRodLength,
+                                 float compressionRatio,
+                                 float transferDuration, float exhaustDuration,
+                                 Crankshaft crankshaft)
+            : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
+        {
+            TransferDuration = transferDuration;
+            ExhaustDuration  = exhaustDuration;
+
+            if (EVO >= IVO)
+                throw new ArgumentException(
+                    $"Exhaust must open before transfer port. " +
+                    $"EVO={EVO:F1}° must be less than IVO={IVO:F1}°. " +
+                    $"Increase exhaustDuration or decrease transferDuration.");
+        }
+
+        // ── Valve lift profile ──────────────────────────────────────────────────
+        /// <summary>
+        /// Smooth trapezoidal lift: fast ramp (15 % of duration), flat top (70 %),
+        /// fast ramp-down (15 %).  Ramps use a smoothstep (3t²-2t³) curve so the
+        /// area derivative is C1-continuous (no kink at ramp/plateau boundaries).
+        /// </summary>
+        private static float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
+        {
+            // Normalise to [0, 360)
+            float deg = thetaDeg % 360f;
+            if (deg < 0f) deg += 360f;
+
+            // Handle wrap-around (e.g. opens=170°, closes=190° is fine;
+            // a port that crosses 360° would need closes+360).
+            float effectiveClose = closes < opens ? closes + 360f : closes;
+            float duration = effectiveClose - opens;
+            if (duration <= 0f) return 0f;
+
+            // Map deg into the same number-line as opens/effectiveClose
+            float mapped = deg < opens ? deg + 360f : deg;
+            if (mapped < opens || mapped > effectiveClose) return 0f;
+
+            float rampDur = duration * RampFraction;
+            float holdEnd = effectiveClose - rampDur;
+
+            if (mapped < opens + rampDur)
+            {
+                // Opening ramp: smoothstep
+                float t = (mapped - opens) / rampDur;
+                return peakLift * t * t * (3f - 2f * t);
+            }
+            else if (mapped <= holdEnd)
+            {
+                // Flat top – full area
+                return peakLift;
+            }
+            else
+            {
+                // Closing ramp: smoothstep reversed
+                float t = (mapped - holdEnd) / rampDur;
+                return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
+            }
+        }
+
+        // ── Cycle event handler ─────────────────────────────────────────────────
+        protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
+        {
+            // ── Fuel injection at transfer-port closing (IVC) ──────────────────
+            // At IVC the cylinder is sealed; whatever air is trapped is what we burn.
+            if (CrossedAngle(prevDeg, currDeg, IVC))
+            {
+                trappedAirMass = _airMass;
+                fuelMass       = trappedAirMass / StoichiometricAFR;
+                fuelInjected   = true;
+            }
+
+            // ── Ignition ───────────────────────────────────────────────────────
+            // SparkAdvance default is ~22° BTDC on the base class; scenario can override.
+            float sparkAngle = (360f - SparkAdvance) % 360f;
+
+            if (CrossedAngle(prevDeg, currDeg, sparkAngle) && !combustionActive && fuelInjected)
+            {
+                if (_random.NextDouble() < MisfireProbability)
+                {
+                    combustionActive = false;
+                }
+                else
+                {
+                    combustionActive = true;
+                    burnFraction     = 0f;
+                    float range = EnergyVariationFraction;
+                    _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
+                }
+            }
+
+            // ── Combustion heat release (Wiebe) ────────────────────────────────
+            if (combustionActive)
+            {
+                float angleSinceSpark = currDeg - sparkAngle;
+                if (angleSinceSpark < 0f) angleSinceSpark += 360f;
+
+                float newFraction = Wiebe(angleSinceSpark);
+                bool burnComplete = newFraction >= 1f
+                                 || angleSinceSpark > WiebeDuration + WiebeStart + SparkAdvance;
+
+                if (burnComplete)
+                {
+                    newFraction      = 1f;
+                    combustionActive = false;
+                    fuelInjected     = false;
+                    float totalMass  = _airMass + _exhaustMass;
+                    _airMass         = 0f;
+                    _exhaustMass     = totalMass;
+                }
+
+                float dFraction = newFraction - burnFraction;
+                if (dFraction > 0f)
+                {
+                    float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction;
+                    cylinderEnergy += dQ;
+                    _exhaustMass   += fuelMass * dFraction;
+                    burnFraction    = newFraction;
+                }
+            }
+        }
+
+        // ── Helper: did the crank cross a target angle this step? ───────────────
+        /// <summary>
+        /// Returns true if the crank swept through <paramref name="target"/> going
+        /// from <paramref name="prev"/> to <paramref name="curr"/> in a single step.
+        /// Handles wrap-around at 360°.
+        /// </summary>
+        private static bool CrossedAngle(float prev, float curr, float target)
+        {
+            // Normal case (no wrap)
+            if (curr >= prev)
+                return prev < target && target <= curr;
+
+            // Wrapped past 360° → two intervals to check
+            return prev < target || target <= curr;
+        }
+    }
+}

Components/Vehicle.cs

@@ -0,0 +1,166 @@
+using System;
+
+namespace FluidSim.Components
+{
+    public class Vehicle
+    {
+        // ---- Gearbox ----
+        public int CurrentGear { get; private set; } = 0;
+        public readonly float[] GearRatios = { 2.5f, 1.8f, 1.4f, 1.1f, 0.9f, 0.75f };
+        public float FinalDriveRatio = 3.0f;
+        public float PrimaryReduction = 2.5f;
+
+        // ---- Clutch ----
+        public float ClutchInput { get; set; }
+        public float ClutchDisengageTime = 0.15f;
+        private float _clutchTimer;
+        private float _currentEngagement = 0f;
+
+        /// <summary>Time constant for clutch engagement smoothing (seconds).</summary>
+        public float EngagementSmoothTime = 0.5f;   // longer, gentler bite
+
+        private float TargetEngagement
+        {
+            get
+            {
+                if (ClutchInput > 0.01f) return 1f - ClutchInput;
+                if (CurrentGear == 0 || _clutchTimer > 0f) return 0f;
+                return 1f;
+            }
+        }
+
+        public float Engagement => _currentEngagement;
+
+        // ---- Clutch torque model ----
+        /// <summary>Peak clutch friction torque (Nm) when fully engaged at high RPM.</summary>
+        public float BaseMaxTorque = 80f;            // much lower than before
+
+        /// <summary>Stiffness when slipping (Nm per rad/s). Lower = softer engagement.</summary>
+        public float ClutchStiffness = 50f;          // very soft
+
+        /// <summary>Below this engine RPM, the clutch torque is progressively reduced to prevent stalling.</summary>
+        public float IdleRpm = 1200f;
+        public float StallPreventionRamp = 300f;     // RPM band above idle where torque ramps up
+
+        // ---- Physical constants ----
+        public float Mass = 160f;
+        public float WheelRadius = 0.32f;
+        public float DragCoefficient = 0.35f;
+        public float FrontalArea = 0.8f;
+        public float AirDensity = 1.225f;
+        public float RollingFrictionCoeff = 0.01f;
+        public float Gravity = 9.81f;
+
+        // ---- State ----
+        public float Speed { get; private set; }
+
+        public (float clutchTorqueOnEngine, float effectiveEngineInertia) Update(float engineRpm, float engineInertia, float dt)
+        {
+            if (_clutchTimer > 0f)
+            {
+                _clutchTimer -= dt;
+                if (_clutchTimer < 0f) _clutchTimer = 0f;
+            }
+
+            float target = TargetEngagement;
+            float smoothing = 1f - MathF.Exp(-dt / Math.Max(EngagementSmoothTime, 0.001f));
+            _currentEngagement += (target - _currentEngagement) * smoothing;
+            if (MathF.Abs(_currentEngagement - target) < 0.001f)
+                _currentEngagement = target;
+
+            float engagement = _currentEngagement;
+
+            float totalGear = 1f;
+            if (CurrentGear > 0)
+                totalGear = GearRatios[CurrentGear - 1] * FinalDriveRatio * PrimaryReduction;
+
+            float engineRadPerSec = engineRpm * 2f * MathF.PI / 60f;
+
+            float v = MathF.Max(Speed, 0f);
+            float drag = 0.5f * AirDensity * DragCoefficient * FrontalArea * v * v;
+            float rolling = RollingFrictionCoeff * Mass * Gravity;
+            float resistanceForce = drag + rolling;
+
+            float clutchTorque = 0f;
+            float effectiveInertia = engineInertia;
+
+            if (engagement > 0f && CurrentGear > 0)
+            {
+                float vehicleReflectedRadPerSec = (Speed / WheelRadius) * totalGear;
+                float slip = engineRadPerSec - vehicleReflectedRadPerSec;
+
+                // Stall prevention: reduce max torque when engine RPM is near idle
+                float torqueLimit = BaseMaxTorque * engagement;
+                if (engineRpm < IdleRpm + StallPreventionRamp)
+                {
+                    float factor = Math.Clamp((engineRpm - IdleRpm) / StallPreventionRamp, 0f, 1f);
+                    torqueLimit *= factor;
+                }
+
+                float stiffnessTorque = ClutchStiffness * engagement * slip;
+                clutchTorque = Math.Clamp(stiffnessTorque, -torqueLimit, torqueLimit);
+
+                // Lock if slip negligible and engagement high
+                if (engagement >= 0.99f && MathF.Abs(slip) < 1.0f)
+                {
+                    float vehicleInertia = Mass * WheelRadius * WheelRadius;
+                    float reflectedVehicleInertia = vehicleInertia / (totalGear * totalGear);
+                    effectiveInertia = engineInertia + reflectedVehicleInertia;
+
+                    Speed = engineRadPerSec * WheelRadius / totalGear;
+                    float loadTorque = resistanceForce * WheelRadius / totalGear;
+                    return (loadTorque, effectiveInertia);
+                }
+            }
+
+            float driveTorqueAtWheel = clutchTorque * totalGear;
+            float driveForce = driveTorqueAtWheel / WheelRadius;
+            float netForce = driveForce - resistanceForce;
+            float acceleration = netForce / Mass;
+            Speed += acceleration * dt;
+            if (Speed < 0f) Speed = 0f;
+
+            return (clutchTorque, engineInertia);
+        }
+
+        public void ShiftUp()
+        {
+            if (CurrentGear < GearRatios.Length)
+            {
+                CurrentGear++;
+                AutoDisengageClutch();
+            }
+        }
+
+        public void ShiftDown()
+        {
+            if (CurrentGear > 1)
+            {
+                CurrentGear--;
+                AutoDisengageClutch();
+            }
+        }
+
+        public void SetNeutral()
+        {
+            CurrentGear = 0;
+            _clutchTimer = 0f;
+        }
+
+        public void SetFirstGear()
+        {
+            if (CurrentGear == 0)
+            {
+                CurrentGear = 1;
+                AutoDisengageClutch();
+            }
+        }
+
+        private void AutoDisengageClutch()
+        {
+            _clutchTimer = ClutchDisengageTime;
+        }
+
+        public float SpeedKmh => Speed * 3.6f;
+    }
+}

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 ---
-            // If 0 when UseInertance is true, a stable default is auto‑computed at runtime.
+            // --- Loss coefficient (linear resistance) ---
             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;
                     }
 
+                    // 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 * area * pipeDxAdj;
+                        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)
-                    float cLocal = mdotFinal >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
-                    float maxULocal = 0.9f * cLocal;
-                    float uCheck = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
-                    if (uCheck > maxULocal)
+                    // ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
+                    if (mdotFinal < 0)
                     {
-                        uFace = maxULocal;
-                        mdotFinal = rhoFace * uFace * area * (mdotFinal >= 0 ? 1f : -1f);
+                        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;
@@ -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;
-                float mdot = rhoGhost * uGhost * area;
-                if (d.IsLeftEnd) mdot = -mdot;
-                d.LastMassFlowRate = mdot;
+                d.LastMassFlowRate = mdotRaw;
                 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
-        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++)
             {
-                    if (f == _pipeEnd[p] && f < _totalCells)
-                    {
-                        _isWallFace[f] = true;
-                        break;
-                    }
-                }
+                _isPipeBoundaryFace[_pipeStart[p]] = true;
+                _isPipeBoundaryFace[_pipeEnd[p]] = true;
             }
 
             // Initialize uniform state
@@ -150,6 +159,7 @@ namespace FluidSim.Core
         public float GetCellPressure(int i) => _p[i];
         public float GetCellDensity(int i) => _rho[i];
         public float GetCellDx(int i) => _dx[i];
+        public float GetCellArea(int i) => _area[i];
         public float GetCellVelocity(int i)
         {
             float rho = _rho[i];
@@ -215,13 +225,13 @@ namespace FluidSim.Core
             }
         }
 
-        // ---------- Flux computation: fuses primitive calculation and flux evaluation ----------
+        // ---------- Flux computation ----------
         private void ComputeFluxes(float dt)
         {
             float fm, fp, fe;
             int vecSize = Vector<float>.Count;
 
-            // ---- 1. Left ghost boundaries ----
+            // ---- 1. Left ghost boundaries → per‑pipe buffers ----
             for (int p = 0; p < _pipeCount; p++)
             {
                 int idx = _pipeStart[p];
@@ -239,22 +249,18 @@ namespace FluidSim.Core
                 float cR = MathF.Sqrt(Gamma * pR * invRhoR);
                 float YR = _Y[idx];
 
-                // store pressure for cell idx
-                _p[idx] = pR;
-
                 LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
-                _fluxM[idx] = fm; _fluxP[idx] = fp; _fluxE[idx] = fe;
+                _leftFluxM[p] = fm; _leftFluxP[p] = fp; _leftFluxE[p] = fe;
 
                 float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
                 ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
-                _fluxY[idx] = fy;
+                _leftFluxY[p] = fy;
             }
 
-            // ---- 2. Right ghost boundaries ----
+            // ---- 2. Right ghost boundaries → per‑pipe buffers ----
             for (int p = 0; p < _pipeCount; p++)
             {
                 int idx = _pipeEnd[p] - 1;
-                int face = idx + 1;
                 int ghostIdx = p * 2 + 1;
                 float rR = _ghost.Rho[ghostIdx];
                 float uR = _ghost.U[ghostIdx];
@@ -269,45 +275,35 @@ namespace FluidSim.Core
                 float cL = MathF.Sqrt(Gamma * pL * invRhoL);
                 float YL = _Y[idx];
 
-                // store pressure for cell idx
-                _p[idx] = pL;
-
                 LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
-                _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
+                _rightFluxM[p] = fm; _rightFluxP[p] = fp; _rightFluxE[p] = fe;
 
                 float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
                 ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
-                _fluxY[face] = fy;
+                _rightFluxY[p] = fy;
             }
 
-            // ---- 3. Interior faces – vectorised SIMD ----
+            // ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
             for (int face = 1; face < _totalCells; face++)
             {
-                // Handle walls (rare) with scalar code
-                if (_isWallFace[face])
-                {
-                    int iL = face - 1;
-                    float rL = _rho[iL], rhouL = _rhou[iL];
-                    float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
-                    float uL = rhouL * invRhoL;
-                    float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
-                    float cL = MathF.Sqrt(Gamma * pL * invRhoL);
-                    _p[iL] = pL;
-
-                    LaxFlux(rL, uL, pL, cL, rL, -uL, pL, cL, out fm, out fp, out fe);
-                    _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
-                    _fluxY[face] = 0f;
+                // Skip faces that belong to a pipe boundary (they are already handled)
+                if (_isPipeBoundaryFace[face])
                     continue;
-                }
 
-                // If the next vecSize faces contain a wall, fall back to scalar for this block
+                // Try to vectorize a block of contiguous non‑boundary faces
                 if (face + vecSize - 1 < _totalCells)
                 {
-                    bool hasWall = false;
+                    bool canVectorize = true;
                     for (int f = face; f < face + vecSize; f++)
-                        if (_isWallFace[f]) { hasWall = true; break; }
+                    {
+                        if (_isPipeBoundaryFace[f])
+                        {
+                            canVectorize = false;
+                            break;
+                        }
+                    }
 
-                    if (!hasWall)
+                    if (canVectorize)
                     {
                         // --- Vectorised block ---
                         var rhoL = new Vector<float>(_rho, face - 1);
@@ -330,11 +326,7 @@ namespace FluidSim.Core
                         var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
                         var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
 
-                        // Store pressures for visualisation (left cell of each face)
-                        pL.CopyTo(_p, face - 1);
-
-                        // Lax‑Friedrichs fluxes
-                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass
+                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
                         var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
 
                         var FmL = rhoL * uL;
@@ -362,50 +354,45 @@ namespace FluidSim.Core
                         feVec.CopyTo(_fluxE, face);
                         fyVec.CopyTo(_fluxY, face);
 
-                        face += vecSize - 1; // loop increment will add 1, so we advance vecSize faces
+                        face += vecSize - 1; // loop increment will add 1
                         continue;
                     }
                 }
 
-                // --- Scalar interior face (fallback) ---
+                // --- Scalar fallback for a single interior face ---
                 {
-                    int iLf = face - 1, iRf = face;
-                    float rLf = _rho[iLf], rhouLf = _rhou[iLf];
-                    float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f));
-                    float uLf = rhouLf * invRhoLf;
-                    float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf);
-                    float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf);
-                    float YLf = _Y[iLf];
-                    _p[iLf] = pLf;
+                    int iL = face - 1, iR = face;
+                    float rL = _rho[iL], rhouL = _rhou[iL];
+                    float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
+                    float uL = rhouL * invRhoL;
+                    float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
+                    float cL = MathF.Sqrt(Gamma * pL * invRhoL);
+                    float YL = _Y[iL];
 
-                    float rRf = _rho[iRf], rhouRf = _rhou[iRf];
-                    float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f));
-                    float uRf = rhouRf * invRhoRf;
-                    float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf);
-                    float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf);
-                    float YRf = _Y[iRf];
+                    float rR = _rho[iR], rhouR = _rhou[iR];
+                    float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
+                    float uR = rhouR * invRhoR;
+                    float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
+                    float cR = MathF.Sqrt(Gamma * pR * invRhoR);
+                    float YR = _Y[iR];
 
-                    LaxFlux(rLf, uLf, pLf, cLf, rRf, uRf, pRf, cRf, out fm, out fp, out fe);
+                    LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
                     _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
 
-                    float alpha = MathF.Max(MathF.Abs(uLf) + cLf, MathF.Abs(uRf) + cRf);
-                    ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy);
+                    float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
+                    ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
                     _fluxY[face] = fy;
                 }
             }
-
-            // If damping/relaxation are active, compute the factors here (re-uses _dampingFactors/_relaxFactors arrays,
-            // but we no longer have a separate precompute pass). We compute them on demand in UpdateCells anyway?
-            // Actually UpdateCells multiplies by these factors; we can compute them there if needed.
         }
 
-        // ---------- Cell update (unchanged core, but skips relaxation/damping when not needed) ----------
+        // ---------- Cell update (per pipe, using correct boundary fluxes) ----------
         private void UpdateCells(float dt)
         {
             int vecSize = Vector<float>.Count;
             float dtRelax = -_relaxRate * dt;
 
-            // Compute damping and relaxation factors if needed
+            // Precompute damping and relaxation factors globally
             if (_applyDamping)
             {
                 for (int i = 0; i < _totalCells; i++)
@@ -418,13 +405,78 @@ 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++)
+            {
+                int start = _pipeStart[p];
+                int end = _pipeEnd[p];    // exclusive
+                int len = end - start;
+                if (len == 0) continue;
+
+                // ------- Left boundary cell (i = start) ------
+                {
+                    int i = start;
+                    float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
+
+                    // left face: always the pipe's left boundary flux
+                    float fluxM_L = _leftFluxM[p];
+                    float fluxP_L = _leftFluxP[p];
+                    float fluxE_L = _leftFluxE[p];
+                    float fluxY_L = _leftFluxY[p];
+
+                    // right face: depends on pipe length
+                    float fluxM_R, fluxP_R, fluxE_R, fluxY_R;
+                    if (len == 1)
+                    {
+                        // Only one cell: right face is the pipe's right boundary flux
+                        fluxM_R = _rightFluxM[p];
+                        fluxP_R = _rightFluxP[p];
+                        fluxE_R = _rightFluxE[p];
+                        fluxY_R = _rightFluxY[p];
+                    }
+                    else
+                    {
+                        // interior face (global flux at index i+1)
+                        fluxM_R = _fluxM[i + 1];
+                        fluxP_R = _fluxP[i + 1];
+                        fluxE_R = _fluxE[i + 1];
+                        fluxY_R = _fluxY[i + 1];
+                    }
+
+                    float dtdx = dt / _dx[i];
+                    float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                    float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                    float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                    float rhoYOld = rhoOld * YOld;
+                    float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                    if (_applyDamping) rhouNew *= _dampingFactors[i];
+                    if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
+
+                    rhoNew = MathF.Max(rhoNew, 1e-12f);
+                    float kin = 0.5f * rhouNew * rhouNew / rhoNew;
+                    float eMin = 100f * Gm1Inv + kin;
+                    ENew = MathF.Max(ENew, eMin);
+
+                    _rho[i] = rhoNew;
+                    _rhou[i] = rhouNew;
+                    _E[i] = ENew;
+                    _Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
+                }
+
+                // ------- Interior cells (i = start+1 to end-2) ------
+                if (len > 2)
+                {
+                    int iCell = start + 1;
+                    int iEnd = end - 1; // exclusive upper bound
+
+                    // Vectorised path for interior cells (if available)
+                    for (; iCell <= iEnd - vecSize; iCell += vecSize)
                     {
                         var rhoOld = new Vector<float>(_rho, iCell);
                         var rhouOld = new Vector<float>(_rhou, iCell);
@@ -471,8 +523,8 @@ namespace FluidSim.Core
                         yNew.CopyTo(_Y, iCell);
                     }
 
-            // Scalar remainder (only a few cells)
-            for (; iCell < _totalCells; 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];
@@ -500,7 +552,70 @@ namespace FluidSim.Core
                     }
                 }
 
-        // ---------- Scalar flux helpers (used in boundaries and scalar fallback) ----------
+                // ------- Right boundary cell (i = end-1, if len > 1) ------
+                if (len > 1)
+                {
+                    int i = end - 1;
+                    float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
+
+                    // left face
+                    float fluxM_L, fluxP_L, fluxE_L, fluxY_L;
+                    if (len == 2)
+                    {
+                        // Only two cells: left face is the pipe's left boundary flux
+                        fluxM_L = _leftFluxM[p];
+                        fluxP_L = _leftFluxP[p];
+                        fluxE_L = _leftFluxE[p];
+                        fluxY_L = _leftFluxY[p];
+                    }
+                    else
+                    {
+                        // interior face (global flux at i)
+                        fluxM_L = _fluxM[i];
+                        fluxP_L = _fluxP[i];
+                        fluxE_L = _fluxE[i];
+                        fluxY_L = _fluxY[i];
+                    }
+
+                    // right face: always the pipe's right boundary flux
+                    float fluxM_R = _rightFluxM[p];
+                    float fluxP_R = _rightFluxP[p];
+                    float fluxE_R = _rightFluxE[p];
+                    float fluxY_R = _rightFluxY[p];
+
+                    float dtdx = dt / _dx[i];
+                    float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                    float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                    float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                    float rhoYOld = rhoOld * YOld;
+                    float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                    if (_applyDamping) rhouNew *= _dampingFactors[i];
+                    if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
+
+                    rhoNew = MathF.Max(rhoNew, 1e-12f);
+                    float kin = 0.5f * rhouNew * rhouNew / rhoNew;
+                    float eMin = 100f * Gm1Inv + kin;
+                    ENew = MathF.Max(ENew, eMin);
+
+                    _rho[i] = rhoNew;
+                    _rhou[i] = rhouNew;
+                    _E[i] = ENew;
+                    _Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
+                }
+            }
+
+            // Recompute pressure for all cells (for visualization)
+            for (int i = 0; i < _totalCells; i++)
+            {
+                float rho = _rho[i];
+                float rhou = _rhou[i];
+                float u = rhou / MathF.Max(rho, 1e-12f);
+                _p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
+            }
+        }
+
+        // ---------- Scalar flux helpers ----------
         private static void LaxFlux(float rL, float uL, float pL, float cL,
                                     float rR, float uR, float pR, float cR,
                                     out float fm, out float fp, out float fe)
@@ -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++)

Program.cs

@@ -33,24 +33,33 @@ public class Program
     // Audio & simulation
     private static SimulationRingBuffer _simRingBuffer = null!;
     private static SoundEngine _soundEngine = null!;
-    private static Scenario _scenario = null!;   // cast to access ThrottleArea
+    private static Scenario _scenario = null!;
     private static Font? _overlayFont;
     private static Text? _overlayText;
 
     // Throttle control
-    private static float _throttleTarget = 1.0f;          // 0‑1, set by arrow keys
-    private static float _throttleCurrent = 0.0f;         // actual current fraction (lerped)
-    private const float ThrottleLerpRate = 10.0f;          // times per second (speed of movement)
+    private static float _throttleTarget = 1.0f;
+    private static float _throttleCurrent = 0.0f;
+    private const float ThrottleLerpRate = 10.0f;
     private static bool _wKeyHeld = false;
     private static float _lastThrottleUpdateTime;
 
+    // Load
+    private static float _loadTarget = 0.0f;          // 0‑1
+    private static float _loadCurrent = 0.0f;
+
+    private static float _clutchTarget = 0f;
+    private static float _clutchCurrent = 0f;
+    private static bool _cKeyHeld = false;
+
     private const int TargetMaxFill = (int)(SampleRate * 0.2);
 
     public static void Main()
     {
         var window = CreateWindow();
         LoadFont();
-        _scenario = new SingleCylScenario();
+        _scenario = new TwoStrokeScenario();
+        _scenario.Font = _overlayFont;
         _scenario.Initialize(SampleRate);
         _lastThrottleUpdateTime = 0.0f;
 
@@ -76,14 +85,12 @@ public class Program
                 (1.0 - Math.Exp(-8.0 * (now - lastDrawTime)));
             _soundEngine.Speed = _currentDisplaySpeed;
 
-            // ---- Throttle update ----
+            // ---- Throttle & Load update (shared dt) ----
             float dtThrottle = (float)now - _lastThrottleUpdateTime;
             _lastThrottleUpdateTime = (float)now;
 
             float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
-
-            // Snap to zero instantly when target is zero (key released)
-            if (throttleDesiredFraction == 0.0)
+            if (throttleDesiredFraction == 0.0f)
             {
                 _throttleCurrent = 0.0f;
             }
@@ -93,8 +100,18 @@ public class Program
                 _throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
             }
 
+            float loadSmoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
+            _loadCurrent += (_loadTarget - _loadCurrent) * loadSmoothing;
+            _scenario.Load = _loadCurrent;
+
             _scenario.Throttle = _throttleCurrent;
 
+            float clutchDesired = _cKeyHeld ? 1f : 0f;
+            float clutchSmoothing = 1f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
+            _clutchCurrent += (clutchDesired - _clutchCurrent) * clutchSmoothing;
+            _scenario.Clutch = _clutchCurrent;
+
+
             // ---- Drawing ----
             if (now - lastDrawTime >= 1.0 / DrawFrequency)
             {
@@ -103,7 +120,8 @@ public class Program
                     string toggleHint = _isRealTime ? "[Space] slow mo" : "[Space] real time";
                     _overlayText.DisplayedString =
                         $"{toggleHint}  Speed: {_currentDisplaySpeed:F3}x  RT: {(_currentDisplaySpeed * 100.0):F1}%  Sim load: {_loadTracker.LoadPercent:F0}%\n" +
-                        $"Throttle: {_throttleCurrent * 100:F0}%  Target: {_throttleTarget * 100:F0}%  [W] {(_wKeyHeld ? "BLIP" : "---")}";
+                        $"Clutch: {_clutchCurrent*100:F0}%  [C]" +
+                        $"Load: {_loadCurrent*100:F0}%  [←][→]  Throttle: {_throttleCurrent * 100:F0}%  Target: {_throttleTarget * 100:F0}%  [W] {(_wKeyHeld ? "BLIP" : "---")}";
                 }
 
                 window.Clear(Color.Black);
@@ -205,6 +223,25 @@ public class Program
             case Keyboard.Key.Down:
                 _throttleTarget = MathF.Max(0.0f, _throttleTarget - 0.05f);
                 break;
+
+            case Keyboard.Key.Left:
+                _loadTarget = MathF.Max(0.0f, _loadTarget - 0.05f);
+                break;
+
+            case Keyboard.Key.Right:
+                _loadTarget = MathF.Min(1.0f, _loadTarget + 0.05f);
+                break;
+
+            case Keyboard.Key.E:
+                _scenario.ShiftUp();
+                break;
+            case Keyboard.Key.Q:
+                _scenario.ShiftDown();
+                break;
+
+            case Keyboard.Key.C:
+                _cKeyHeld = true;
+                break;
         }
     }
 
@@ -212,5 +249,8 @@ public class Program
     {
         if (e.Code == Keyboard.Key.W)
             _wKeyHeld = false;
+
+        if (e.Code == Keyboard.Key.C)
+            _cKeyHeld = false;
     }
 }

Scenarios/Scenario.cs

@@ -2,6 +2,8 @@
 using SFML.System;
 using FluidSim.Core;
 using FluidSim.Components;
+using System;
+using System.Collections.Generic;
 
 namespace FluidSim.Tests
 {
@@ -10,11 +12,204 @@ namespace FluidSim.Tests
         protected const float AmbientPressure = 101325f;
         protected const float AmbientTemperature = 300f;
         public float Throttle { get; set; }
+        public float Load { get; set; }
+        public float Clutch { get; set; }   // 0 = engaged, 1 = fully disengaged (manual lever)
+        public Font? Font { get; set; }
 
         public abstract void Initialize(int sampleRate);
         public abstract float Process();
         public abstract void Draw(RenderWindow target);
 
+        public virtual void ShiftUp() { }
+        public virtual void ShiftDown() { }     
+
+        // ---- Dyno curve graph ----
+        private const float RpmBinSize = 50f;
+        private readonly List<(float powerKw, float torqueNm)> _dynoBins = new();
+        private int _lastDynoBin = -1;
+
+        public void ResetDynoCurve()
+        {
+            _dynoBins.Clear();
+            _lastDynoBin = -1;
+        }
+
+        protected void UpdateDynoCurve(float rpm, float powerKw, float torqueNm)
+        {
+            if (rpm <= 0) return;
+            int bin = (int)(rpm / RpmBinSize);
+
+            while (_dynoBins.Count <= bin)
+                _dynoBins.Add((0f, 0f));
+
+            if (_lastDynoBin >= 0 && bin > _lastDynoBin + 1)
+            {
+                float lastPower = _dynoBins[_lastDynoBin].powerKw > 0 ? _dynoBins[_lastDynoBin].powerKw : 0f;
+                float lastTorque = _dynoBins[_lastDynoBin].torqueNm > 0 ? _dynoBins[_lastDynoBin].torqueNm : 0f;
+                for (int b = _lastDynoBin + 1; b < bin; b++)
+                {
+                    float t = (b - _lastDynoBin) / (float)(bin - _lastDynoBin);
+                    float interpPower = lastPower + (powerKw - lastPower) * t;
+                    float interpTorque = lastTorque + (torqueNm - lastTorque) * t;
+                    if (interpPower > _dynoBins[b].powerKw || _dynoBins[b].powerKw <= 0)
+                        _dynoBins[b] = (interpPower, _dynoBins[b].torqueNm);
+                    if (interpTorque > _dynoBins[b].torqueNm || _dynoBins[b].torqueNm <= 0)
+                        _dynoBins[b] = (_dynoBins[b].powerKw, interpTorque);
+                }
+            }
+
+            var current = _dynoBins[bin];
+            if (powerKw > current.powerKw || current.powerKw <= 0)
+                current.powerKw = powerKw;
+            if (torqueNm > current.torqueNm || current.torqueNm <= 0)
+                current.torqueNm = torqueNm;
+            _dynoBins[bin] = current;
+
+            _lastDynoBin = bin;
+        }
+
+        protected void DrawDynoCurve(RenderWindow target,
+                                     float graphX, float graphY, float graphWidth, float graphHeight,
+                                     float currentRpm, float currentPowerKw)
+        {
+            if (_dynoBins.Count == 0) return;
+
+            float maxPowerKw = 0.01f, maxTorqueNm = 0.01f, maxRpm = 1000f;
+            for (int b = 0; b < _dynoBins.Count; b++)
+            {
+                var bin = _dynoBins[b];
+                if (bin.powerKw > 0 || bin.torqueNm > 0)
+                {
+                    float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
+                    if (bin.powerKw > maxPowerKw) maxPowerKw = bin.powerKw;
+                    if (bin.torqueNm > maxTorqueNm) maxTorqueNm = bin.torqueNm;
+                    if (rpmBin > maxRpm) maxRpm = rpmBin;
+                }
+            }
+            maxPowerKw *= 1.1f;
+            maxTorqueNm *= 1.1f;
+            maxRpm = MathF.Max(maxRpm * 1.05f, 1000f);
+
+            var bg = new RectangleShape(new Vector2f(graphWidth, graphHeight))
+            {
+                FillColor = new Color(20, 20, 20, 200),
+                Position = new Vector2f(graphX, graphY)
+            };
+            target.Draw(bg);
+
+            const float leftMargin = 50f, rightMargin = 50f, topMargin = 20f, bottomMargin = 35f;
+            float plotX = graphX + leftMargin;
+            float plotY = graphY + topMargin;
+            float plotW = graphWidth - leftMargin - rightMargin;
+            float plotH = graphHeight - topMargin - bottomMargin;
+
+            float xMin = 0f, xMax = maxRpm;
+            float yLeftMin = 0f, yLeftMax = maxPowerKw;
+            float yRightMin = 0f, yRightMax = maxTorqueNm;
+
+            var powerColor = new Color(0xFF, 0x1B, 0x1B);
+            var torqueColor = new Color(0x09, 0x09, 0xFF);
+            var gridColor = new Color(50, 50, 50);
+
+            for (int i = 0; i <= 9; i++)
+            {
+                float t = i / 9f;
+                float x = plotX + t * plotW;
+                var vLine = new VertexArray(PrimitiveType.Lines, 2);
+                vLine[0] = new Vertex(new Vector2f(x, plotY), gridColor);
+                vLine[1] = new Vertex(new Vector2f(x, plotY + plotH), gridColor);
+                target.Draw(vLine);
+            }
+            for (int i = 0; i <= 5; i++)
+            {
+                float t = i / 5f;
+                float y = plotY + (1 - t) * plotH;
+                var hLine = new VertexArray(PrimitiveType.Lines, 2);
+                hLine[0] = new Vertex(new Vector2f(plotX, y), gridColor);
+                hLine[1] = new Vertex(new Vector2f(plotX + plotW, y), gridColor);
+                target.Draw(hLine);
+            }
+
+            DrawLabel(target, "RPM", new Vector2f(graphX + graphWidth / 2 - 12, graphY + graphHeight - 15), Color.White, 12);
+            DrawLabel(target, "kW", new Vector2f(graphX + 5, graphY + 2), Color.White, 11);
+            DrawLabel(target, "Nm", new Vector2f(graphX + graphWidth - 25, graphY + 2), Color.White, 11);
+
+            for (int i = 0; i <= 5; i++)
+            {
+                float leftValue = yLeftMin + (yLeftMax - yLeftMin) * i / 5f;
+                float rightValue = yRightMin + (yRightMax - yRightMin) * i / 5f;
+                float y = plotY + (1 - i / 5f) * plotH;
+                DrawLabel(target, $"{leftValue:F1}", new Vector2f(graphX + 2, y - 6), Color.White, 9);
+                DrawLabel(target, $"{rightValue:F1}", new Vector2f(graphX + graphWidth - 40, y - 6), Color.White, 9);
+            }
+
+            for (int i = 0; i <= 9; i++)
+            {
+                float value = xMin + (xMax - xMin) * i / 9f;
+                float x = plotX + i / 9f * plotW;
+                DrawLabel(target, $"{value / 1000f:F1}k", new Vector2f(x - 15, graphY + graphHeight - bottomMargin + 5), Color.White, 9);
+            }
+
+            var powerLine = new VertexArray(PrimitiveType.LineStrip);
+            bool firstPower = true;
+            for (int b = 0; b < _dynoBins.Count; b++)
+            {
+                float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
+                if (rpmBin > xMax) break;
+                var bin = _dynoBins[b];
+                if (bin.powerKw > 0)
+                {
+                    float sx = plotX + (rpmBin - xMin) / (xMax - xMin) * plotW;
+                    float sy = plotY + (1 - (bin.powerKw - yLeftMin) / (yLeftMax - yLeftMin)) * plotH;
+                    if (firstPower) { powerLine.Clear(); firstPower = false; }
+                    powerLine.Append(new Vertex(new Vector2f(sx, sy), powerColor));
+                }
+                else if (!firstPower)
+                {
+                    target.Draw(powerLine);
+                    powerLine.Clear();
+                    firstPower = true;
+                }
+            }
+            if (!firstPower) target.Draw(powerLine);
+
+            var torqueLine = new VertexArray(PrimitiveType.LineStrip);
+            bool firstTorque = true;
+            for (int b = 0; b < _dynoBins.Count; b++)
+            {
+                float rpmBin = b * RpmBinSize + RpmBinSize / 2f;
+                if (rpmBin > xMax) break;
+                var bin = _dynoBins[b];
+                if (bin.torqueNm > 0)
+                {
+                    float sx = plotX + (rpmBin - xMin) / (xMax - xMin) * plotW;
+                    float sy = plotY + (1 - (bin.torqueNm - yRightMin) / (yRightMax - yRightMin)) * plotH;
+                    if (firstTorque) { torqueLine.Clear(); firstTorque = false; }
+                    torqueLine.Append(new Vertex(new Vector2f(sx, sy), torqueColor));
+                }
+                else if (!firstTorque)
+                {
+                    target.Draw(torqueLine);
+                    torqueLine.Clear();
+                    firstTorque = true;
+                }
+            }
+            if (!firstTorque) target.Draw(torqueLine);
+
+            if (currentRpm > 0 && currentRpm <= xMax && currentPowerKw > 0)
+            {
+                float sx = plotX + (currentRpm - xMin) / (xMax - xMin) * plotW;
+                float sy = plotY + (1 - (currentPowerKw - yLeftMin) / (yLeftMax - yLeftMin)) * plotH;
+                var dot = new CircleShape(2.5f)
+                {
+                    FillColor = Color.White,
+                    Position = new Vector2f(sx - 2.5f, sy - 2.5f)
+                };
+                target.Draw(dot);
+            }
+        }
+
+        // ---- Drawing helpers ----
         protected Color PressureColor(float pressurePa)
         {
             float bar = pressurePa / 1e5f;
@@ -68,7 +263,7 @@ namespace FluidSim.Tests
             target.Draw(border);
         }
 
-        protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
+        protected void DrawCylinder(RenderWindow target, EngineCylinder cylinder,
                                     float centerX, float topY, float width, float maxHeight)
         {
             float fraction = cylinder.PistonFraction;
@@ -107,7 +302,8 @@ namespace FluidSim.Tests
         }
 
         protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex,
-                                float pipeCenterY, float pipeStartX, float pipeEndX)
+                        float pipeCenterY, float pipeStartX, float pipeEndX,
+                        float areaScale = 0f)
         {
             int start = pipeSystem.GetPipeStart(pipeIndex);
             int end = pipeSystem.GetPipeEnd(pipeIndex);
@@ -116,20 +312,34 @@ namespace FluidSim.Tests
 
             float pipeLen = pipeEndX - pipeStartX;
             float dx = pipeLen / (n - 1);
-            float baseRadius = 25f;
 
             var centers = new float[n];
             var radii = new float[n];
             var temps = new float[n];
+
             for (int i = 0; i < n; i++)
             {
                 int cell = start + i;
                 float p = pipeSystem.GetCellPressure(cell);
                 float rho = pipeSystem.GetCellDensity(cell);
                 temps[i] = p / MathF.Max(rho * 287f, 1e-12f);
+
+                if (areaScale > 0f)
+                {
+                    // Use actual cell area to determine visual radius
+                    float area = pipeSystem.GetCellArea(cell);
+                    radii[i] = MathF.Sqrt(area / MathF.PI) * areaScale;
+                    if (radii[i] < 1f) radii[i] = 1f;
+                }
+                else
+                {
+                    // Original pressure‑based radius
                     float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
+                    float baseRadius = 25f;   // default visual radius for constant‑area pipes
                     radii[i] = baseRadius * (1f + dev * 2f);
                     if (radii[i] < 2f) radii[i] = 2f;
+                }
+
                 centers[i] = pipeStartX + i * dx;
             }
 
@@ -157,5 +367,18 @@ namespace FluidSim.Tests
             }
             target.Draw(va);
         }
+
+        protected void DrawLabel(RenderWindow target, string text, Vector2f position, Color fillColor, uint characterSize = 14)
+        {
+            if (Font == null) return;
+            var txt = new Text(Font)
+            {
+                DisplayedString = text,
+                Position = position,
+                FillColor = fillColor,
+                CharacterSize = characterSize
+            };
+            target.Draw(txt);
+        }
     }
 }

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,200 @@ 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, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
-        private int intakeOpenIdx, exhaustOpenIdx;
+        private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx;
         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) ----------
-        // Bore 56 mm, Stroke 57 mm, CR 9.5
-        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)
-        private const float IntakeValveDiam = 0.027f;
-        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;
+        private float _maxThrottleArea;
+        private float intakePipeArea, exhaustPipeArea;
+        private const float MaxBrakeTorque = 30.0f;   // Nm at full load
 
         public override void Initialize(int sampleRate)
         {
             dt = 1.0 / sampleRate;
 
-            // ---------- Crankshaft ----------
-            crankshaft = new Crankshaft(600);            // start at ~600 RPM
-            crankshaft.Inertia = 0.2f;
-            crankshaft.FrictionConstant = 2.0f;
-            crankshaft.FrictionViscous = 0.04f;
+            // Throttle body diameter 44mm (typical for 250cc MX)
+            _maxThrottleArea = (float)Units.AreaFromDiameter(44 * Units.mm);
 
-            // ---------- Cylinder ----------
-            cylinder = new Cylinder(Bore, Stroke, ConRod, CompressionRatio,
-                                    IVO, IVC, EVO, EVC, crankshaft)
+            // ---- Crankshaft ----
+            crankshaft = new Crankshaft(2000);
+            crankshaft.Inertia = 0.02f;          // kg·m² (crank + flywheel)
+            crankshaft.FrictionConstant = 3.0f;  // Nm – bearings, rings, seals
+            crankshaft.FrictionViscous = 0.002f; // Nm/(rad/s) – oil windage
+
+            // ---- Cylinder (CRF250R) ----
+            float bore = 0.078f;          // 78 mm
+            float stroke = 0.0522f;       // 52.2 mm → 249.4 cc
+            float conRod = 0.1044f;       // 2× stroke
+            float compRatio = 13.5f;      // typical
+
+            // Valve events (high‑performance MX cam)
+            float ivo = 340f, ivc = 600f; // intake opens 20° BTDC (overlap), closes 60° ABDC
+            float evo = 120f, evc = 380f; // exhaust opens 60° BBDC, closes 20° ATDC
+
+            cylinder = new Cylinder(bore, stroke, conRod, compRatio,
+                                    ivo, ivc, evo, evc, crankshaft)
             {
-                IntakeValveDiameter = IntakeValveDiam,
-                ExhaustValveDiameter = ExhaustValveDiam,
-                IntakeValveLift = ValveLift,
-                ExhaustValveLift = ValveLift,
-                SparkAdvance = SparkAdv,
-                EnergyVariationFraction = 0.03f,   // small cycle‑to‑cycle variation
-                MisfireProbability = 0.0f
+                IntakeValveDiameter = 0.036f,   // 36 mm
+                IntakeValveLift = 0.0095f,      // 9.5 mm
+                ExhaustValveDiameter = 0.030f,  // 30 mm
+                ExhaustValveLift = 0.0085f      // 8.5 mm
             };
 
-            // ---------- Pipe system ----------
-            int totalCells = CellsBefore + CellsRunner + CellsExhaust;
-            int[] pipeStart = { 0, CellsBefore, CellsBefore + CellsRunner };
-            int[] pipeEnd   = { CellsBefore, CellsBefore + CellsRunner, totalCells };
+            // ---- Pipe system ----
+            int[] pipeStart = { 0, 10, 20 };
+            int[] pipeEnd   = { 10, 20, 70 };
+            int totalCells = pipeEnd[^1];
+            float[] area = new float[totalCells];
+            float[] dx   = new float[totalCells];
 
-            float[] areas = new float[totalCells];
-            float[] dxs   = new float[totalCells];
-            float dxBefore = IntakeLenBefore / CellsBefore;
-            float dxRunner = RunnerLen / CellsRunner;
-            float dxExh    = ExhaustLen / CellsExhaust;
+            float intakeDia = 0.040f;     // 40 mm intake runner
+            float exhaustDia = 0.038f;    // 38 mm exhaust primary
+            intakePipeArea  = MathF.PI * 0.25f * intakeDia * intakeDia;
+            exhaustPipeArea = MathF.PI * 0.25f * exhaustDia * exhaustDia;
+
+            float intakeLenBefore = 0.15f;  // throttle body to plenum
+            float intakeLenRunner = 0.25f;  // plenum to valve
+            float exhaustLen      = 0.50f;  // exhaust length
 
             for (int i = 0; i < totalCells; i++)
             {
-                areas[i] = PipeArea;
-                if (i < CellsBefore)
-                    dxs[i] = dxBefore;
-                else if (i < CellsBefore + CellsRunner)
-                    dxs[i] = dxRunner;
+                if (i < 10)
+                {
+                    area[i] = intakePipeArea; dx[i] = intakeLenBefore / 10f;
+                }
+                else if (i < 20)
+                {
+                    area[i] = intakePipeArea; dx[i] = intakeLenRunner / 10f;
+                }
                 else
-                    dxs[i] = dxExh;
+                {
+                    area[i] = exhaustPipeArea; dx[i] = exhaustLen / 50f;
+                }
             }
 
-            float rho0 = 101325f / (287f * 300f);
-            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, areas, dxs,
-                                        rho0, 0f, 101325f);
-            pipeSystem.DampingMultiplier = 0.5f;
-            pipeSystem.EnergyRelaxationRate = 0f;        // adiabatic pipes
+            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
+                                        1.225f, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 1.0f;
+            pipeSystem.EnergyRelaxationRate = 0.5f;
             pipeSystem.AmbientPressure = 101325f;
 
-            // ---------- Volumes ----------
-            intakePlenum = new Volume0D(PlenumVolume, 101325f, 300f);
+            // ---- Volumes ----
+            intakePlenum = new Volume0D(1.0e-3f, 101325f, 300f);   // 1 litre airbox
             plenumInlet = intakePlenum.CreatePort();
             plenumOutlet = intakePlenum.CreatePort();
+            exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // unused
+            colIn = exhaustCollector.CreatePort();
+            colOut = exhaustCollector.CreatePort();
 
-            // ---------- Boundary system ----------
+            // ---- Boundary system ----
             boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
 
             throttleAreaIdx = 0;
-            plenumRunnerIdx  = 1;
+            plenumRunnerAreaIdx = 1;
             intakeValveIdx = 2;
             exhaustValveIdx = 3;
 
-            // Open ends
-            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, PipeArea);
+            // Open ends (pipe area = pipe cross‑section)
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, intakePipeArea);
             intakeOpenIdx = 0;
-            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, PipeArea);
+            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, exhaustPipeArea);
             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);
+            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.7f);      // throttle
+            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1.0f); // plenum→runner
+            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1.0f); // intake valve
+            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1.0f); // exhaust valve
 
             orificeAreas = new float[4];
-            orificeAreas[plenumRunnerIdx] = PipeArea;   // fixed full‑bore
+            orificeAreas[plenumRunnerAreaIdx] = intakePipeArea; // runner cross‑section (fixed)
 
-            // ---------- Solver ----------
-            solver = new Solver { SubStepCount = 5, EnableProfiling = false };
+            // ---- Solver ----
+            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 ----------
-            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
-            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
+            // ---- Sound ----
+            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
+            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
             reverb = new OutdoorExhaustReverb(sampleRate);
 
             stepCount = 0;
-            Console.WriteLine("Single‑cylinder engine (YX140) ready.");
+            Console.WriteLine("CRF250R engine ready.");
         }
 
         public override float Process()
         {   
-            // ---- Crank and cylinder pre‑step ----
+            // Manual brake torque (0..30 Nm)
+            float loadTorque = Load * MaxBrakeTorque;
+            crankshaft.SetLoadTorque(loadTorque);
+
             crankshaft.Step((float)dt);
             cylinder.PreStep((float)dt);
 
-            // ---- Update variable areas ----
-            float throttledArea = MaxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1.0f);
+            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
             orificeAreas[throttleAreaIdx] = throttledArea;
+
             orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
             orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
             boundaries.SetOrificeAreas(orificeAreas);
 
-            // ---- Fluids step ----
             solver.Step();
             stepCount++;
 
-            // ---- Sound ----
             float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
             float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
 
             float exhaustDry = exhaustSound.Process(exhaustFlow);
             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, " +
-                                  $"Throttle={Throttle * 100:F0}%");
+                float crankDeg = (crankshaft.CrankAngle + cylinder.PhaseOffset) * 180f / MathF.PI % 720f;
+                Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1}°, RPM={rpm:F0}, CylP={cylinder.Pressure/1e5f:F2} bar");
+                Console.WriteLine($"  intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
+
+                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");
+
+                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");
+
+                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");
+
+                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) * 0.5f);
         }
 
         public override void Draw(RenderWindow target)
@@ -220,53 +213,54 @@ namespace FluidSim.Tests
 
             float intakeY = winH / 2f - 40f;
             float exhaustY = winH / 2f + 80f;
-            float leftX = 40f;
+            float openEndX = 40f;
 
-            // Intake open end marker
-            var om = new CircleShape(5f) { FillColor = Color.Cyan };
-            om.Position = new Vector2f(leftX - 5f, intakeY - 5f);
-            target.Draw(om);
+            float pipe1StartX = openEndX;
+            float pipe1EndX = pipe1StartX + 120f;
+            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;
-            var thr = new RectangleShape(new Vector2f(8f, 30f))
+            float throttleX = pipe1EndX + 5f;
+            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
-            float plenW = 60f, plenH = 50f;
-            float plenLeftX = thrX + 12f;
+            float plenW = 60f, plenH = 80f;
+            float plenLeftX = throttleX + 10f;
             float plenCenterX = plenLeftX + plenW / 2f;
             float plenTopY = intakeY - plenH / 2f;
             DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
 
-            // Pipe 1 – runner
-            float rStartX = plenLeftX + plenW + 10f;
-            float rEndX = rStartX + 100f;
-            DrawPipe(target, pipeSystem, 1, intakeY, rStartX, rEndX);
+            float runnerStartX = plenLeftX + plenW + 5f;
+            float runnerEndX = runnerStartX + 100f;
+            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
             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);
+            // --- RPM & Power labels ---
+            float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            float powerKw = crankshaft.AveragePower * 1e-3f;
+            DrawLabel(target, $"RPM: {rpm:F0}", new Vector2f(20, 90), Color.White, 24);
+            DrawLabel(target, $"Power: {powerKw:F2} kW", new Vector2f(20, 115), Color.White, 24);
+
+            // --- Dyno curve ---
+            float torqueNm = crankshaft.AverageTorque;
+            UpdateDynoCurve(rpm, powerKw, torqueNm);
+
+            float graphX = winW - 410f;
+            float graphY = winH - 260f;
+            float graphW = 400f;
+            float graphH = 250f;
+            DrawDynoCurve(target, graphX, graphY, graphW, graphH, rpm, powerKw);
         }
     }
 }

Scenarios/TwoStrokeScenario.cs

@@ -0,0 +1,350 @@
+using FluidSim.Components;
+using FluidSim.Core;
+using FluidSim.Interfaces;
+using FluidSim.Utils;
+using SFML.Graphics;
+using SFML.System;
+using System;
+
+namespace FluidSim.Tests
+{
+    public class TwoStrokeScenario : Scenario
+    {
+        private Crankshaft crankshaft;
+        private TwoStrokeCylinder cylinder;
+
+        private PipeSystem pipeSystem;
+        private BoundarySystem boundaries;
+        private Solver solver;
+
+        private Volume0D intakePlenum;
+        private Port plenumInlet, plenumOutlet;
+        private Volume0D exhaustMuffler;
+        private Port mufflerIn, mufflerOut;
+
+        private Vehicle vehicle;
+
+        private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
+        private float[] orificeAreas;
+        private int intakeOpenIdx, exhaustOpenIdx;
+
+        private SoundProcessor exhaustSound, intakeSound;
+        private OutdoorExhaustReverb reverb;
+
+        private double dt;
+        private int stepCount;
+
+        private float _maxThrottleArea;
+        private float intakePipeArea, exhaustHeaderArea;
+
+        public override void ShiftUp()   => vehicle.ShiftUp();
+        public override void ShiftDown() => vehicle.ShiftDown();
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            // ── Vehicle ──────────────────────────────────────────────────────────
+            vehicle = new Vehicle();
+
+            // ── Throttle body: 42 mm – wider to reduce high-RPM intake restriction ──
+            _maxThrottleArea = (float)Units.AreaFromDiameter(42 * Units.mm);
+
+            // ── Crankshaft ───────────────────────────────────────────────────────
+            // Lighter flywheel for quicker revving; friction tuned to ~0.5 kW loss at idle
+            crankshaft = new Crankshaft(2000);
+            crankshaft.CycleLength      = 2f * MathF.PI;   // two-stroke: fire every rev
+            crankshaft.Inertia          = 0.06f;            // lighter flywheel
+            crankshaft.FrictionConstant = 0.4f;             // ~0.4 Nm constant drag
+            crankshaft.FrictionViscous  = 0.0004f;          // ~2.5 Nm at 10 000 RPM
+
+            // ── Cylinder: 125 cc, motocross-style two-stroke ─────────────────────
+            // Bore × stroke = 54 × 54.5 mm → 124.9 cc
+            float bore         = 0.054f;
+            float stroke       = 0.0545f;
+            float conRod       = 0.110f;   // ~2× stroke
+            float compRatio    = 7.2f;     // geometric CR; effective CR after port closure is ~12:1
+
+            // Port timings: exhaust 195°, transfer 155° – competitive MX 125
+            float transferDuration = 155f;
+            float exhaustDuration  = 195f;
+
+            cylinder = new TwoStrokeCylinder(bore, stroke, conRod, compRatio,
+                                             transferDuration, exhaustDuration,
+                                             crankshaft)
+            {
+                IntakeValveDiameter  = 0.042f,   // matched to intake pipe
+                IntakeValveLift      = 0.015f,
+                ExhaustValveDiameter = 0.040f,
+                ExhaustValveLift     = 0.013f
+            };
+
+            // ── Pipe geometry ────────────────────────────────────────────────────
+            //
+            // Layout (all lengths in mm):
+            //   Intake path:  airbox stub 100 mm  |  runner 180 mm
+            //   Exhaust path: expansion chamber tuned to ~9 000 RPM power peak
+            //     header      170 mm  Ø 40 mm
+            //     diffuser    280 mm  Ø 40 → 72 mm
+            //     belly       200 mm  Ø 72 mm
+            //     convergent  130 mm  Ø 72 → 28 mm
+            //     stinger      70 mm  Ø 28 mm
+            //     total       850 mm
+            //
+            // Cell sizing: ~14 mm/cell.
+            // CFL: c_sound ≈ 550 m/s, dx=0.014 m → dt_max ≈ 25 µs
+            //      at 44100 Hz dt = 22.7 µs → SubStepCount=4 keeps CFL safely ≤ 1
+
+            // --- Cell counts ---
+            int intakeCells  = 7;    // 100 mm stub  → ~14 mm/cell
+            int runnerCells  = 13;   // 180 mm runner → ~14 mm/cell
+            int exhaustCells = 60;   // 850 mm total  → ~14 mm/cell
+
+            int totalCells = intakeCells + runnerCells + exhaustCells;
+            int[] pipeStart = { 0,           intakeCells,              intakeCells + runnerCells };
+            int[] pipeEnd   = { intakeCells, intakeCells + runnerCells, totalCells };
+
+            float[] area = new float[totalCells];
+            float[] dx   = new float[totalCells];
+
+            // --- Intake ---
+            float intakeDia      = 0.042f;   // matches throttle body
+            float intakeStubLen  = 0.100f;
+            float intakeRunnerLen= 0.160f;   // shorter runner → less pumping loss
+            intakePipeArea = MathF.PI * 0.25f * intakeDia * intakeDia;
+
+            for (int i = 0; i < intakeCells; i++)
+            { area[i] = intakePipeArea; dx[i] = intakeStubLen / intakeCells; }
+
+            for (int i = intakeCells; i < intakeCells + runnerCells; i++)
+            { area[i] = intakePipeArea; dx[i] = intakeRunnerLen / runnerCells; }
+
+            // Expansion chamber tuned for ~8 500 RPM power peak.
+            // Return-pulse travel distance = 0.5 × c_avg × (60 / RPM_target)
+            // c_avg ≈ 480 m/s → distance = 0.5 × 480 × (60/8500) ≈ 1.69 m round-trip
+            // → one-way pipe length ≈ 0.84 m  (matches total below)
+            float headerDia    = 0.040f;  float headerLen    = 0.130f;   // shorter header → earlier pulse
+            float diffEndDia   = 0.070f;  float diffuserLen  = 0.250f;   // slightly narrower belly
+            float bellyDia     = 0.070f;  float bellyLen     = 0.220f;
+            float convEndDia   = 0.028f;  float convergentLen= 0.160f;   // longer convergent → stronger return pulse
+            float stingerDia   = 0.028f;  float stingerLen   = 0.080f;
+            // total = 0.13+0.25+0.22+0.16+0.08 = 0.84 m
+
+            exhaustHeaderArea  = MathF.PI * 0.25f * headerDia  * headerDia;
+            float bellyArea    = MathF.PI * 0.25f * bellyDia   * bellyDia;
+            float stingerArea  = MathF.PI * 0.25f * stingerDia * stingerDia;
+
+            // Distribute cells proportionally by section length
+            int headerCells     = Math.Max(1, (int)MathF.Round(exhaustCells * headerLen    / 0.84f));
+            int diffuserCells   = Math.Max(1, (int)MathF.Round(exhaustCells * diffuserLen  / 0.84f));
+            int bellyCells      = Math.Max(1, (int)MathF.Round(exhaustCells * bellyLen     / 0.84f));
+            int convergentCells = Math.Max(1, (int)MathF.Round(exhaustCells * convergentLen/ 0.84f));
+            int stingerCells    = exhaustCells - headerCells - diffuserCells
+                                               - bellyCells  - convergentCells;
+            if (stingerCells < 1) stingerCells = 1;
+
+            int exhBase = intakeCells + runnerCells;
+            int idx = 0;
+            for (int i = exhBase; i < totalCells; i++, idx++)
+            {
+                if (idx < headerCells)
+                {
+                    area[i] = exhaustHeaderArea;
+                    dx[i]   = headerLen / headerCells;
+                }
+                else if (idx < headerCells + diffuserCells)
+                {
+                    float t  = (idx - headerCells) / (float)(diffuserCells - 1);
+                    // Smooth cosine taper instead of linear for better wave reflection
+                    float ct = 0.5f * (1f - MathF.Cos(MathF.PI * t));
+                    float dia = headerDia + (diffEndDia - headerDia) * ct;
+                    area[i] = MathF.PI * 0.25f * dia * dia;
+                    dx[i]   = diffuserLen / diffuserCells;
+                }
+                else if (idx < headerCells + diffuserCells + bellyCells)
+                {
+                    area[i] = bellyArea;
+                    dx[i]   = bellyLen / bellyCells;
+                }
+                else if (idx < headerCells + diffuserCells + bellyCells + convergentCells)
+                {
+                    float t  = (idx - headerCells - diffuserCells - bellyCells)
+                               / (float)(convergentCells - 1);
+                    // Steeper cosine convergent for a sharper return pulse
+                    float ct = 0.5f * (1f - MathF.Cos(MathF.PI * t));
+                    float dia = bellyDia + (convEndDia - bellyDia) * ct;
+                    area[i] = MathF.PI * 0.25f * dia * dia;
+                    dx[i]   = convergentLen / convergentCells;
+                }
+                else
+                {
+                    area[i] = stingerArea;
+                    dx[i]   = stingerLen / stingerCells;
+                }
+            }
+
+            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
+                                        1.225f, 0f, 101325f);
+            pipeSystem.DampingMultiplier    = 0.8f;   // slightly less damping → stronger pulses
+            pipeSystem.EnergyRelaxationRate = 0.4f;
+            pipeSystem.AmbientPressure      = 101325f;
+
+            // ── 0-D Volumes ──────────────────────────────────────────────────────
+            // Intake plenum: acts as a small airbox resonator (8 cc)
+            intakePlenum = new Volume0D(8e-3f, 101325f, 300f);
+            plenumInlet  = intakePlenum.CreatePort();
+            plenumOutlet = intakePlenum.CreatePort();
+
+            // Exhaust silencer volume: 600 cc is realistic for a small-bore muffler
+            exhaustMuffler = new Volume0D(600e-6f, 101325f, 650f);
+            mufflerIn  = exhaustMuffler.CreatePort();
+            mufflerOut = exhaustMuffler.CreatePort();
+
+            // ── Boundary system ───────────────────────────────────────────────────
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
+            throttleAreaIdx  = 0;
+            plenumRunnerIdx  = 1;
+            intakeValveIdx   = 2;
+            exhaustValveIdx  = 3;
+
+            // Open ends: atmosphere at both extremes
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true,  101325f, intakePipeArea);
+            intakeOpenIdx  = 0;
+            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, stingerArea);
+            exhaustOpenIdx = 1;
+
+            // Orifices: throttle → plenum → runner → cylinder → exhaust pipe
+            boundaries.AddOrifice(plenumInlet,          0, false, throttleAreaIdx, 0.72f);
+            boundaries.AddOrifice(plenumOutlet,         1, true,  plenumRunnerIdx, 1.00f);
+            boundaries.AddOrifice(cylinder.IntakePort,  1, false, intakeValveIdx,  0.68f);
+            boundaries.AddOrifice(cylinder.ExhaustPort, 2, true,  exhaustValveIdx, 0.70f);
+
+            orificeAreas = new float[4];
+            orificeAreas[plenumRunnerIdx] = intakePipeArea;  // runner always fully open
+
+            // ── Solver ────────────────────────────────────────────────────────────
+            // SubStepCount = 4 keeps CFL ≤ 1 for 5 mm cells at 44 100 Hz
+            solver = new Solver { SubStepCount = 4, EnableProfiling = false };
+            solver.SetTimeStep(dt);
+            solver.SetPipeSystem(pipeSystem);
+            solver.SetBoundarySystem(boundaries);
+            solver.AddComponent(cylinder);
+            solver.AddComponent(intakePlenum);
+            solver.AddComponent(exhaustMuffler);
+
+            // ── Sound ─────────────────────────────────────────────────────────────
+            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
+            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 4.5f };
+            reverb = new OutdoorExhaustReverb(sampleRate);
+
+            stepCount = 0;
+            Console.WriteLine("125cc Two-Stroke – expansion chamber tuned for ~8 500 RPM power peak");
+            Console.WriteLine($"  Exhaust cells: {exhaustCells}  |  header {headerCells}  diffuser {diffuserCells}" +
+                              $"  belly {bellyCells}  convergent {convergentCells}  stinger {stingerCells}");
+        }
+
+        public override float Process()
+        {
+            float engineRpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+
+            vehicle.ClutchInput = Clutch;
+
+            var (clutchTorque, effectiveInertia) = vehicle.Update(engineRpm, crankshaft.Inertia, (float)dt);
+            crankshaft.SetEffectiveInertia(effectiveInertia);
+            crankshaft.SetLoadTorque(clutchTorque);
+
+            crankshaft.Step((float)dt);
+            cylinder.PreStep((float)dt);
+
+            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
+            orificeAreas[throttleAreaIdx]  = throttledArea;
+            orificeAreas[intakeValveIdx]   = cylinder.IntakeValveArea;
+            orificeAreas[exhaustValveIdx]  = cylinder.ExhaustValveArea;
+            boundaries.SetOrificeAreas(orificeAreas);
+
+            solver.Step();
+            stepCount++;
+
+            float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
+            float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
+
+            float exhaustDry = exhaustSound.Process(exhaustFlow);
+            float intakeDry  = intakeSound.Process(intakeFlow);
+
+            if (stepCount % 2000 == 0)
+            {
+                float rpm     = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+                float powerKw = crankshaft.AveragePower * 1e-3f;
+                float torqueNm = crankshaft.AverageTorque;
+                Console.WriteLine($"Step {stepCount,7} | RPM={rpm,6:F0} | Power={powerKw,5:F2} kW" +
+                                  $" | Torque={torqueNm,5:F1} Nm | Gear={vehicle.CurrentGear}" +
+                                  $" | Speed={vehicle.SpeedKmh,4:F0} km/h");
+            }
+
+            return reverb.Process((intakeDry + exhaustDry) * 0.5f);
+        }
+
+        // ── Drawing ───────────────────────────────────────────────────────────────
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+
+            float intakeY  = winH / 2f - 40f;
+            float exhaustY = winH / 2f + 80f;
+            float openEndX = 40f;
+
+            // Intake stub
+            float x = openEndX;
+            float w = 120f;
+            DrawPipe(target, pipeSystem, 0, intakeY, x, x + w);
+
+            // Throttle body
+            float throttleX = x + w + 5f;
+            var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
+            {
+                FillColor = Color.Yellow,
+                Position  = new Vector2f(throttleX, intakeY - 15f)
+            };
+            target.Draw(throttleRect);
+
+            // Plenum
+            float plenW = 40f, plenH = 60f;
+            float plenX = throttleX + 10f;
+            DrawVolume(target, intakePlenum, plenX + plenW / 2f, intakeY - plenH / 2f, plenW, plenH);
+
+            // Runner
+            float runnerStartX = plenX + plenW + 5f;
+            DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerStartX + 100f);
+
+            // Cylinder
+            float cylCX  = runnerStartX + 150f;
+            float cylTopY = intakeY - 120f;
+            DrawCylinder(target, cylinder, cylCX, cylTopY, 80f, 240f);
+
+            // Exhaust pipe (expansion chamber)
+            float exhStartX = cylCX + 40f + 20f;
+            DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, winW - 60f, areaScale: 800f);
+
+            // HUD labels
+            float rpm     = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+            float powerKw = crankshaft.AveragePower * 1e-3f;
+            float torqueNm = crankshaft.AverageTorque;
+
+            DrawLabel(target, $"RPM: {rpm:F0}",           new Vector2f(20, 90),  Color.White, 24);
+            DrawLabel(target, $"Power: {powerKw:F2} kW",  new Vector2f(20, 115), Color.White, 24);
+            DrawLabel(target, $"Torque: {torqueNm:F1} Nm",new Vector2f(20, 140), Color.White, 20);
+
+            string gearText = vehicle.CurrentGear == 0 ? "N" : vehicle.CurrentGear.ToString();
+            DrawLabel(target, $"Gear: {gearText}",        new Vector2f(20, 162), Color.Cyan, 20);
+            DrawLabel(target, $"Speed: {vehicle.SpeedKmh:F0} km/h",
+                                                          new Vector2f(20, 184), Color.Cyan, 20);
+            DrawLabel(target, vehicle.Engagement > 0.99f ? "Clutch: Locked" : "Clutch: Slipping",
+                                                          new Vector2f(20, 204), Color.Cyan, 14);
+
+            // Dyno curve
+            UpdateDynoCurve(rpm, powerKw, torqueNm);
+            DrawDynoCurve(target, winW - 410f, winH - 260f, 400f, 250f, rpm, powerKw);
+        }
+    }
+}