added valves display

.vscode/launch.json

@@ -0,0 +1,26 @@
+{
+    "version": "0.2.0",
+    "configurations": [
+        {
+            // Use IntelliSense to find out which attributes exist for C# debugging
+            // Use hover for the description of the existing attributes
+            // For further information visit https://github.com/dotnet/vscode-csharp/blob/main/debugger-launchjson.md
+            "name": ".NET Core Launch (console)",
+            "type": "coreclr",
+            "request": "launch",
+            "preLaunchTask": "build",
+            // If you have changed target frameworks, make sure to update the program path.
+            "program": "${workspaceFolder}/bin/Debug/net10.0/FluidSim.dll",
+            "args": [],
+            "cwd": "${workspaceFolder}",
+            // For more information about the 'console' field, see https://aka.ms/VSCode-CS-LaunchJson-Console
+            "console": "internalConsole",
+            "stopAtEntry": false
+        },
+        {
+            "name": ".NET Core Attach",
+            "type": "coreclr",
+            "request": "attach"
+        }
+    ]
+}

.vscode/tasks.json

@@ -0,0 +1,41 @@
+{
+    "version": "2.0.0",
+    "tasks": [
+        {
+            "label": "build",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "build",
+                "${workspaceFolder}/FluidSim.csproj",
+                "/property:GenerateFullPaths=true",
+                "/consoleloggerparameters:NoSummary;ForceNoAlign"
+            ],
+            "problemMatcher": "$msCompile"
+        },
+        {
+            "label": "publish",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "publish",
+                "${workspaceFolder}/FluidSim.csproj",
+                "/property:GenerateFullPaths=true",
+                "/consoleloggerparameters:NoSummary;ForceNoAlign"
+            ],
+            "problemMatcher": "$msCompile"
+        },
+        {
+            "label": "watch",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "watch",
+                "run",
+                "--project",
+                "${workspaceFolder}/FluidSim.csproj"
+            ],
+            "problemMatcher": "$msCompile"
+        }
+    ]
+}

Components/Crankshaft.cs

@@ -0,0 +1,54 @@
+// Components/Crankshaft.cs
+using System;
+
+namespace FluidSim.Components
+{
+    public class Crankshaft
+    {
+        public double AngularVelocity { get; set; }   // rad/s
+        public double CrankAngle { get; set; }        // rad, 0 … 4π (four‑stroke cycle)
+        public double PreviousAngle { get; set; }     // ← now has public setter
+
+        public double Inertia { get; set; } = 0.2;
+        public double FrictionConstant { get; set; } = 2.0;      // N·m
+        public double FrictionViscous { get; set; } = 0.005;     // N·m per rad/s
+
+        private double externalTorque;
+
+        public Crankshaft(double initialRPM = 400.0)
+        {
+            AngularVelocity = initialRPM * 2.0 * Math.PI / 60.0;
+            CrankAngle = 0.0;
+            PreviousAngle = 0.0;
+        }
+
+        public void AddTorque(double torque) => externalTorque += torque;
+
+        public void Step(double dt)
+        {
+            // Catch NaN before it propagates
+            if (double.IsNaN(AngularVelocity) || double.IsInfinity(AngularVelocity))
+                AngularVelocity = 0.0;
+            if (double.IsNaN(externalTorque) || double.IsInfinity(externalTorque))
+                externalTorque = 0.0;
+
+            PreviousAngle = CrankAngle;
+
+            double friction = FrictionConstant * Math.Sign(AngularVelocity) + FrictionViscous * AngularVelocity;
+            double netTorque = externalTorque - friction;
+            double alpha = netTorque / Inertia;
+            AngularVelocity += alpha * dt;
+
+            if (AngularVelocity < 0) AngularVelocity = 0;
+
+            CrankAngle += AngularVelocity * dt;
+
+            if (CrankAngle >= 4.0 * Math.PI)
+                CrankAngle -= 4.0 * Math.PI;
+            else if (CrankAngle < 0)
+                CrankAngle += 4.0 * Math.PI;
+
+            externalTorque = 0.0;
+        }
+    }
+}

Components/EngineCylinder.cs

@@ -0,0 +1,249 @@
+using System;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+    public class EngineCylinder
+    {
+        public Volume0D Cylinder { get; private set; }
+        private Crankshaft crankshaft;
+
+        private double bore, stroke, conRodLength, compressionRatio;
+        private double pistonArea;
+
+        public double V_disp { get; private set; }
+        public double V_clear { get; private set; }
+        public bool ignition = false;
+
+        // ---- Exhaust valve ----
+        private double exhMaxOrificeArea;
+        private double exhValveOpenStart = 130.0 * Math.PI / 180.0;
+        private double exhValveOpenEnd   = 390.0 * Math.PI / 180.0;
+        private double exhValveRampWidth = 30.0  * Math.PI / 180.0;
+        public double ExhaustOrificeArea => ExhaustValveLift() * exhMaxOrificeArea;
+        public double ExhaustValveLiftCurrent => ExhaustValveLift();
+
+        // ---- Intake valve ----
+        private double intMaxOrificeArea;
+        private double intValveOpenStart = 340.0 * Math.PI / 180.0;
+        private double intValveOpenEnd   = 600.0 * Math.PI / 180.0;
+        private double intValveRampWidth = 30.0  * Math.PI / 180.0;
+        public double IntakeOrificeArea => IntakeValveLift() * intMaxOrificeArea;
+        public double IntakeValveLiftCurrent => IntakeValveLift();
+
+        // ---- Combustion ----
+        public double TargetPeakPressure { get; set; } = 50.0 * 101325.0;
+        private const double PeakTemperature = 2500.0;
+        private bool burnInProgress = false;
+        private double burnStartAngle;         // cycle angle (0–4π)
+        private double burnDuration = 40.0 * Math.PI / 180.0;
+        private double targetBurnEnergy;
+        private double preIgnitionMass, preIgnitionInternalEnergy;
+
+        private Random rand = new Random();
+        public double MisfireProbability { get; set; } = 0.02;
+        private bool misfireCurrent = false;
+
+        public int CombustionCount { get; private set; }
+        public int MisfireCount { get; private set; }
+
+        // Cycle‑aware angle (0 – 4π)
+        public double CycleAngle => crankshaft.CrankAngle % (4.0 * Math.PI);
+        private double prevCycleAngle;
+
+        // Piston position fraction (0 = TDC, 1 = BDC)
+        public double PistonPositionFraction
+        {
+            get
+            {
+                double currentVol = Cylinder.Volume;
+                if (currentVol <= V_clear) return 0.0;
+                if (currentVol >= V_clear + V_disp) return 1.0;
+                return (currentVol - V_clear) / V_disp;
+            }
+        }
+
+        public EngineCylinder(Crankshaft crankshaft,
+                              double bore, double stroke, double compressionRatio,
+                              double exhPipeArea, double intPipeArea, int sampleRate)
+        {
+            this.crankshaft = crankshaft;
+            this.bore = bore;
+            this.stroke = stroke;
+            conRodLength = 2.0 * stroke;
+            this.compressionRatio = compressionRatio;
+            exhMaxOrificeArea = exhPipeArea * 0.5;
+            intMaxOrificeArea = intPipeArea * 0.5;
+            pistonArea = Math.PI / 4.0 * bore * bore;
+
+            V_disp = pistonArea * stroke;
+            V_clear = V_disp / (compressionRatio - 1.0);
+
+            // Start at BDC with fresh ambient charge
+            double V_bdc = V_clear + V_disp;
+            double p_amb = 101325.0;
+            double T_amb = 300.0;
+            double rho0 = p_amb / (287.0 * T_amb);
+            double mass0 = rho0 * V_bdc;
+            double energy0 = p_amb * V_bdc / (1.4 - 1.0);
+
+            Cylinder = new Volume0D(V_bdc, p_amb, T_amb, sampleRate)
+            {
+                Gamma = 1.4,
+                GasConstant = 287.0
+            };
+            Cylinder.Volume = V_bdc;
+            Cylinder.Mass = mass0;
+            Cylinder.InternalEnergy = energy0;
+
+            prevCycleAngle = CycleAngle;
+
+            preIgnitionMass = Cylinder.Mass;
+            preIgnitionInternalEnergy = Cylinder.InternalEnergy;
+        }
+
+        // ---- Piston kinematics ----
+        private (double volume, double dvdt) PistonKinematics(double cycleAngle)
+        {
+            double theta = cycleAngle % (2.0 * Math.PI);
+            double R = stroke / 2.0;
+            double cosT = Math.Cos(theta);
+            double sinT = Math.Sin(theta);
+            double L = conRodLength;
+
+            double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT);
+            double V = V_clear + pistonArea * s;
+
+            double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
+            double dVdθ = pistonArea * (R * sinT + (R * R * sinT * cosT) / sqrtTerm);
+            double dvdt = dVdθ * crankshaft.AngularVelocity;
+            return (V, dvdt);
+        }
+
+        // ---- Valve lifts (cycle‑aware) ----
+        private double ExhaustValveLift()
+        {
+            double a = CycleAngle;
+            if (a < exhValveOpenStart || a > exhValveOpenEnd) return 0.0;
+            double duration = exhValveOpenEnd - exhValveOpenStart;
+            double ramp = exhValveRampWidth;
+            double t = (a - exhValveOpenStart) / duration;
+            double rampFrac = ramp / duration;
+            if (t < rampFrac) return t / rampFrac;
+            if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
+            return 1.0;
+        }
+
+        private double IntakeValveLift()
+        {
+            double a = CycleAngle;
+            if (a < intValveOpenStart || a > intValveOpenEnd) return 0.0;
+            double duration = intValveOpenEnd - intValveOpenStart;
+            double ramp = intValveRampWidth;
+            double t = (a - intValveOpenStart) / duration;
+            double rampFrac = ramp / duration;
+            if (t < rampFrac) return t / rampFrac;
+            if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
+            return 1.0;
+        }
+
+        // ---- Wiebe burn fraction ----
+        private double WiebeFraction(double angleFromIgnition)
+        {
+            if (angleFromIgnition >= burnDuration) return 1.0;
+            double a = 5.0, m = 2.0;
+            double x = angleFromIgnition / burnDuration;
+            return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1));
+        }
+
+        // ---- Torque from pressure ----
+        private double ComputeTorque()
+        {
+            double p = Cylinder.Pressure;
+            double ambient = 101325.0;
+            double force = (p - ambient) * pistonArea;
+            if (force <= 0) return 0.0;
+
+            double theta = crankshaft.CrankAngle % (2.0 * Math.PI);
+            double R = stroke / 2.0;
+            double L = conRodLength;
+            double sinT = Math.Sin(theta);
+            double cosT = Math.Cos(theta);
+
+            double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
+            double lever = R * sinT * (1.0 + (R * cosT) / sqrtTerm);
+            return force * lever;
+        }
+
+        // ---- TDC detection (power stroke, at angle 0 mod 4π) ----
+        private bool DetectTDCPowerStroke()
+        {
+            double current = CycleAngle;
+            double previous = prevCycleAngle;
+            prevCycleAngle = current;
+            return (previous > 3.8 * Math.PI && current < 0.2 * Math.PI);
+        }
+
+        public void Step(double dt)
+        {
+            bool crossingTDC = DetectTDCPowerStroke();
+
+            if (crossingTDC)
+            {
+                misfireCurrent = rand.NextDouble() < MisfireProbability;
+
+                // *** Always capture the state at TDC, whether we burn or not ***
+                preIgnitionMass = Cylinder.Mass;
+                preIgnitionInternalEnergy = Cylinder.InternalEnergy;
+
+                if (misfireCurrent)
+                {
+                    MisfireCount++;
+                }
+                else if (ignition)
+                {
+                    double V = Cylinder.Volume;
+                    targetBurnEnergy = TargetPeakPressure * V / (Cylinder.Gamma - 1.0);
+                    if (double.IsNaN(targetBurnEnergy))
+                        targetBurnEnergy = 101325.0 * V / (Cylinder.Gamma - 1.0);
+                    burnInProgress = true;
+                    burnStartAngle = CycleAngle;
+                    CombustionCount++;
+                }
+            }
+
+            if (burnInProgress)
+            {
+                double angleFromIgnition = CycleAngle - burnStartAngle;
+                if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI;
+
+                if (angleFromIgnition >= burnDuration)
+                {
+                    Cylinder.InternalEnergy = targetBurnEnergy;
+                    burnInProgress = false;
+                }
+                else
+                {
+                    double fraction = WiebeFraction(angleFromIgnition);
+                    Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction)
+                                              + targetBurnEnergy * fraction;
+                    Cylinder.Mass = preIgnitionMass;
+                }
+            }
+
+            var (vol, dvdt) = PistonKinematics(CycleAngle);
+            Cylinder.Volume = vol;
+            Cylinder.Dvdt = dvdt;
+
+            if (double.IsNaN(Cylinder.Pressure) || double.IsNaN(Cylinder.Temperature) || Cylinder.Mass < 1e-9)
+            {
+                double V = Math.Max(vol, V_clear);
+                Cylinder.Mass = 1.225 * V;
+                Cylinder.InternalEnergy = 101325.0 * V / (1.4 - 1.0);
+            }
+
+            double torque = ComputeTorque();
+            crankshaft.AddTorque(torque);
+        }
+    }
+}

Components/Nozzleflow.cs

@@ -0,0 +1,46 @@
+using System;
+
+namespace FluidSim.Components
+{
+    public static class NozzleFlow
+    {
+        /// <summary>
+        /// Computes the nozzle‑exit primitive state and mass flow rate from a
+        /// volume to a pipe, using isentropic relations. Follows ensim4's flow() logic.
+        /// </summary>
+        public static void Compute(double Pt_high, double Tt_high,
+                                   double P_low, double gamma, double R, double area,
+                                   out double rhoExit, out double uExit,
+                                   out double pExit, out double mdot)
+        {
+            double gm1 = gamma - 1.0;
+            double Pt_over_Ps = Pt_high / P_low;
+
+            // Mach number (subsonic, clamped to 1)
+            double M = Math.Sqrt(Math.Max(0.0,
+                (2.0 / gm1) * (Math.Pow(Pt_over_Ps, gm1 / gamma) - 1.0)));
+            if (M > 1.0) M = 1.0;
+
+            double T_star = Tt_high / (1.0 + 0.5 * gm1 * M * M);
+            double a_star = Math.Sqrt(gamma * R * T_star);
+            double u_star = M * a_star;
+            pExit = Pt_high * Math.Pow(1.0 + 0.5 * gm1 * M * M, -gamma / gm1);
+            rhoExit = pExit / (R * T_star);
+            uExit = u_star;   // positive away from high‑pressure side
+            mdot = rhoExit * uExit * area;
+        }
+
+        /// <summary>
+        /// Ambient cell for non‑reflecting open end (ensim4 calc_ambient_cell).
+        /// </summary>
+        public static void ComputeAmbientCell(double rhoInt, double uInt, double pInt,
+                                               double pAmbient, double gamma,
+                                               out double rhoAmb, out double uAmb,
+                                               out double pAmb)
+        {
+            pAmb = pAmbient;
+            uAmb = uInt;
+            rhoAmb = rhoInt * Math.Pow(pAmb / pInt, 1.0 / gamma);
+        }
+    }
+}

Components/Pipe1D.cs

@@ -1,13 +1,15 @@
 using System;
+using System.Numerics;
 using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
     public enum BoundaryType
     {
-        VolumeCoupling,
         OpenEnd,
-        ClosedEnd
+        ZeroPressureOpen,
+        ClosedEnd,
+        GhostCell
     }
 
     public class Pipe1D
@@ -17,404 +19,598 @@ namespace FluidSim.Components
         public double Area => _area;
         public double DampingMultiplier { get; set; } = 1.0;
 
-        private int _n;
-        private double _dx, _dt, _gamma, _area, _diameter;
-        private double[] _rho, _rhou, _E;
+        private int _n;                     // number of real cells
+        private float _dx, _dt;            // spatial step, global time step
+        private float _area, _diameter;    // cross‑section, diameter
+        private float _gamma;              // ratio of specific heats (1.4)
 
-        // Volume‑coupling ghost states for boundaries A and B
-        private double _rhoA, _pA;
-        private double _rhoB, _pB;
-        private bool _aBCSet, _bBCSet;
+        // Conserved variables – arrays sized [_n] (only real cells, ghosts handled externally)
+        private float[] _rho;
+        private float[] _rhou;
+        private float[] _E;
 
-        private BoundaryType _aBCType = BoundaryType.VolumeCoupling;
-        private BoundaryType _bBCType = BoundaryType.VolumeCoupling;
+        // Flux arrays for faces 0 .. _n (face i is between cell i-1 and i)
+        private float[] _fluxM;   // mass flux
+        private float[] _fluxP;   // momentum flux
+        private float[] _fluxE;   // energy flux
 
-        private double _aAmbientPressure = 101325.0;
-        private double _bAmbientPressure = 101325.0;
+        // Ghost cell states
+        private float _rhoGhostL, _uGhostL, _pGhostL;
+        private float _rhoGhostR, _uGhostR, _pGhostR;
+        private bool _ghostLSet, _ghostRSet;
 
-        private const double CflTarget = 0.8;
-        private const double ReferenceSoundSpeed = 340.0;
+        private BoundaryType _aBCType = BoundaryType.GhostCell;
+        private BoundaryType _bBCType = BoundaryType.GhostCell;
 
-        public int GetCellCount() => _n;
-        public double GetCellDensity(int i) => _rho[i];
-        public double GetCellPressure(int i) => Pressure(i);
-        public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
+        private float _aAmbientPressure = 101325f;
+        private float _bAmbientPressure = 101325f;
 
-        public BoundaryType ABCType => _aBCType;
-        public BoundaryType BBCType => _bBCType;
+        // CFL / sub-stepping
+        private const float CflTarget = 0.8f;
+        private const float ReferenceSoundSpeed = 340f;
+        private float _lastMassFlow = 0f;
+
+        // Pre‑computed for damping
+        private float _laminarCoeff;  // 8*mu / r^2, then multiplied by DampingMultiplier
+
+        // ---- Energy loss (Newton cooling) ----
+        private float _ambientEnergyReference;   // total energy density at ambient (Pamb / (γ-1))
+        public float EnergyRelaxationRate { get; set; } = 0.0f;   // 1/s
 
         public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
         {
-            double dtGlobal = 1.0 / sampleRate;
+            float dtGlobal = 1f / sampleRate;
             int nCells;
+            float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget;
 
             if (forcedCellCount > 1)
-            {
                 nCells = forcedCellCount;
-            }
             else
             {
-                double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
-                nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
-                while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
+                nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero));
+                while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1)
                     nCells++;
             }
 
             _n = nCells;
-            _dx = length / _n;
+            _dx = (float)(length / nCells);
             _dt = dtGlobal;
-            _area = area;
-            _gamma = 1.4;
+            _area = (float)area;
+            _diameter = (float)(2.0 * Math.Sqrt(area / Math.PI));
+            _gamma = 1.4f;
 
-            // Hydraulic diameter for a circular pipe
-            _diameter = 2.0 * Math.Sqrt(area / Math.PI);
+            _rho = new float[_n];
+            _rhou = new float[_n];
+            _E = new float[_n];
 
-            _rho = new double[_n];
-            _rhou = new double[_n];
-            _E = new double[_n];
+            // +1 because there are _n+1 faces
+            _fluxM = new float[_n + 1];
+            _fluxP = new float[_n + 1];
+            _fluxE = new float[_n + 1];
+
+            // Pre‑compute laminar damping coefficient (using air at 20°C)
+            float mu_air = 1.8e-5f;
+            float radius = _diameter * 0.5f;
+            _laminarCoeff = 8f * mu_air / (radius * radius); // will be multiplied by DampingMultiplier at each step
+
+            // Ambient reference energy (internal energy per unit volume at 101325 Pa)
+            _ambientEnergyReference = 101325f / (_gamma - 1f);   // ≈ 253312.5 J/m³
 
             PortA = new Port();
             PortB = new Port();
         }
 
+        // ==================== PUBLIC API ============================
         public void SetABoundaryType(BoundaryType type) => _aBCType = type;
         public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
-        public void SetAAmbientPressure(double p) => _aAmbientPressure = p;
-        public void SetBAmbientPressure(double p) => _bAmbientPressure = p;
+        public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
+        public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p;
+
+        public float GetFaceMassFlux(int faceIndex)
+        {
+            if (faceIndex < 0 || faceIndex > _n) return 0f;
+            return _fluxM[faceIndex];
+        }
+
+        public void SetGhostLeft(double rho, double u, double p)
+        {
+            _rhoGhostL = (float)rho;
+            _uGhostL = (float)u;
+            _pGhostL = (float)p;
+            _ghostLSet = true;
+        }
+        public void SetGhostRight(double rho, double u, double p)
+        {
+            _rhoGhostR = (float)rho;
+            _uGhostR = (float)u;
+            _pGhostR = (float)p;
+            _ghostRSet = true;
+        }
+        public void ClearGhostFlag()
+        {
+            _ghostLSet = false;
+            _ghostRSet = false;
+        }
 
         public void SetUniformState(double rho, double u, double p)
         {
-            double e = p / ((_gamma - 1) * rho);
-            double Etot = rho * e + 0.5 * rho * u * u;
+            float r = (float)rho;
+            float vel = (float)u;
+            float pr = (float)p;
+            float e = pr / ((_gamma - 1f) * r);
+            float Etot = r * e + 0.5f * r * vel * vel;
             for (int i = 0; i < _n; i++)
             {
-                _rho[i] = rho;
-                _rhou[i] = rho * u;
+                _rho[i] = r;
+                _rhou[i] = r * vel;
                 _E[i] = Etot;
             }
         }
 
-        public void SetCellState(int i, double rho, double u, double p)
+        public int GetCellCount() => _n;
+        public double GetCellDensity(int i) => _rho[i];
+        public double GetCellVelocity(int i)
         {
-            if (i < 0 || i >= _n) return;
-            _rho[i] = rho;
-            _rhou[i] = rho * u;
-            double e = p / ((_gamma - 1) * rho);
-            _E[i] = rho * e + 0.5 * rho * u * u;
+            float rho = Math.Max(_rho[i], 1e-12f);
+            return _rhou[i] / rho;
         }
-
-        public void SetAVolumeState(double rhoVol, double pVol)
+        public double GetCellPressure(int i)
         {
-            _rhoA = rhoVol;
-            _pA = pVol;
-            _aBCSet = true;
-        }
-
-        public void SetBVolumeState(double rhoVol, double pVol)
-        {
-            _rhoB = rhoVol;
-            _pB = pVol;
-            _bBCSet = true;
-        }
-
-        public void ClearBC() => _aBCSet = _bBCSet = false;
-
-        public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
-        {
-            double maxW = 0.0;
-            for (int i = 0; i < _n; i++)
-            {
-                double rho = _rho[i];
-                double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
-                double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
-                double local = u + c;
-                if (local > maxW) maxW = local;
-            }
-            maxW = Math.Max(maxW, 1e-8);
-            return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
-        }
-
-        public void SimulateSingleStep(double dtSub)
-        {
-            int n = _n;
-            double[] Fm = new double[n + 1];
-            double[] Fp = new double[n + 1];
-            double[] Fe = new double[n + 1];
-
-            // ---------- Boundary A (face 0, left) ----------
-            double rhoIntA = _rho[0];
-            double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
-            double pIntA = Pressure(0);
-
-            switch (_aBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_aBCSet)
-                    {
-                        HLLCFlux(_rhoA, 0.0, _pA,
-                                 rhoIntA, uIntA, pIntA,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    }
-                    else
-                    {
-                        Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false,
-                                  out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-            }
-
-            // ---------- Internal faces ----------
-            for (int i = 0; i < n - 1; i++)
-            {
-                double rhoL = _rho[i];
-                double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
-                double pL = Pressure(i);
-
-                double rhoR = _rho[i + 1];
-                double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
-                double pR = Pressure(i + 1);
-
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR,
-                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
-            }
-
-            // ---------- Boundary B (face n, right) ----------
-            double rhoIntB = _rho[n - 1];
-            double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12);
-            double pIntB = Pressure(n - 1);
-
-            switch (_bBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_bBCSet)
-                    {
-                        HLLCFlux(rhoIntB, uIntB, pIntB,
-                                 _rhoB, 0.0, _pB,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    }
-                    else
-                    {
-                        Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true,
-                                  out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-            }
-
-            // ---- Cell update with linear laminar damping ----
-            double radius = _diameter / 2.0;
-            double mu_air = 1.8e-5;
-            double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
-
-            for (int i = 0; i < n; i++)
-            {
-                double dM = (Fm[i + 1] - Fm[i]) / _dx;
-                double dP = (Fp[i + 1] - Fp[i]) / _dx;
-                double dE = (Fe[i + 1] - Fe[i]) / _dx;
-
-                _rho[i] -= dtSub * dM;
-                _rhou[i] -= dtSub * dP;
-                _E[i] -= dtSub * dE;
-
-                double rho = Math.Max(_rho[i], 1e-12);
-                double dampingRate = laminarCoeff / rho;
-                double dampingFactor = Math.Exp(-dampingRate * dtSub);
-                _rhou[i] *= dampingFactor;
-
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
-                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
-                double pMin = 100.0;
-                double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
-                if (_E[i] < eMin) _E[i] = eMin;
-            }
-
-            // ---------- Port quantities ----------
-            double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
-            double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
-
-            PortA.MassFlowRate = mdotA_sub;
-            PortB.MassFlowRate = mdotB_sub;
-            PortA.Pressure = pIntA;
-            PortB.Pressure = pIntB;
-            PortA.Density = _rho[0];
-            PortB.Density = _rho[n - 1];
-
-            // Corrected enthalpy for both directions
-            if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet)
-            {
-                PortA.SpecificEnthalpy = mdotA_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(0)
-                    : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12);
-            }
-            if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet)
-            {
-                PortB.SpecificEnthalpy = mdotB_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(_n - 1)
-                    : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12);
-            }
-        }
-
-        private double GetCellTotalSpecificEnthalpy(int i)
-        {
-            double rho = Math.Max(_rho[i], 1e-12);
-            double u = _rhou[i] / rho;
-            double p = Pressure(i);
-            double h = _gamma / (_gamma - 1.0) * p / rho;
-            return h + 0.5 * u * u;
-        }
-
-        private double Pressure(int i) =>
-            (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
-
-        // ========== Characteristic‑based Open End ==========
-        private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
-            if (uInt <= -cInt)   // supersonic inflow – use interior state as ghost
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt <= 0)   // subsonic inflow
-            {
-                // Reservoir condition: p = pAmb, T = 300K, u = 0
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-                return;
-            }
-            else   // subsonic outflow (uInt > 0)
-            {
-                // Ghost pressure forced to pAmb
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary)
-                double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
-
-                // Prevent spurious inflow by clipping to zero
-                if (uGhost < 0) uGhost = 0;
-
-                HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-            }
-        }
-
-        private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            if (uInt >= cInt)   // supersonic outflow (extrapolation)
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt >= 0)   // subsonic outflow
-            {
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary)
-                double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
-
-                // Clip to zero to prevent inflow
-                if (uGhost > 0) uGhost = 0;
-
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe);
-            }
-            else   // subsonic inflow
-            {
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe);
-            }
-        }
-
-        // ========== Closed end (mirror) ==========
-        private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary,
-                                   out double fm, out double fp, out double fe)
-        {
-            double rhoGhost = rhoInt;
-            double pGhost = pInt;
-            double uGhost = -uInt;   // mirror velocity
-
-            if (isRightBoundary)
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
-            else
-                HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
-        }
-
-        // ========== Standard HLLC flux ==========
-        private void HLLCFlux(double rL, double uL, double pL,
-                              double rR, double uR, double pR,
-                              out double fm, out double fp, out double fe)
-        {
-            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12));
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12));
-            double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
-            double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR;
-            double SL = Math.Min(uL - cL, uR - cR);
-            double SR = Math.Max(uL + cL, uR + cR);
-
-            double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
-                      / (rL * (SL - uL) - rR * (SR - uR));
-
-            double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
-            double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
-
-            if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
-            else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
-            else if (Ss >= 0)
-            {
-                double rsL = rL * (SL - uL) / (SL - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
-                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
-                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
-            }
-            else
-            {
-                double rsR = rR * (SR - uR) / (SR - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
-                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
-                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
-            }
+            float rho = Math.Max(_rho[i], 1e-12f);
+            return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
         }
 
         public double GetPressureAtFraction(double fraction)
         {
             int i = (int)(fraction * (_n - 1));
             i = Math.Clamp(i, 0, _n - 1);
-            return Pressure(i);
+            return GetCellPressure(i);
+        }
+
+        public void SetCellState(int i, double rho, double u, double p)
+        {
+            if (i < 0 || i >= _n) return;
+            float r = (float)rho;
+            float vel = (float)u;
+            float pr = (float)p;
+            _rho[i] = r;
+            _rhou[i] = r * vel;
+            float e = pr / ((_gamma - 1f) * r);
+            _E[i] = r * e + 0.5f * r * vel * vel;
+        }
+
+        public double GetOpenEndMassFlow()
+        {
+            if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen)
+                return 0.0;
+
+            int lastCell = _n - 1;
+            float rho = _rho[lastCell];
+            float u = _rhou[lastCell] / Math.Max(rho, 1e-12f);
+            float p = PressureScalar(lastCell);
+
+            float c = MathF.Sqrt(_gamma * p / rho);
+            float uFace = u;
+            float rhoFace = rho;
+            float pFace = p;
+
+            if (uFace > 0 && uFace < c) // subsonic outflow
+            {
+                float s = p / MathF.Pow(rho, _gamma);
+                float rhoAmb = MathF.Pow(_bAmbientPressure / s, 1f / _gamma);
+                float cAmb = MathF.Sqrt(_gamma * _bAmbientPressure / rhoAmb);
+                float J_plus = u + 2f * c / (_gamma - 1f);
+                float uFaceNew = J_plus - 2f * cAmb / (_gamma - 1f);
+                if (uFaceNew > 0) uFace = uFaceNew;
+                if (uFace < 0) uFace = 0;
+                rhoFace = rhoAmb;
+                pFace = _bAmbientPressure;
+            }
+
+            return rhoFace * uFace * _area;
+        }
+
+        public double GetAndStoreMassFlowForDerivative()
+        {
+            float current = (float)GetOpenEndMassFlow();
+            double derivative = (current - _lastMassFlow) / _dt;
+            _lastMassFlow = current;
+            return derivative;
+        }
+
+        public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8f)
+        {
+            float maxW = 0f;
+            for (int i = 0; i < _n; i++)
+            {
+                float rho = _rho[i];
+                float u = MathF.Abs(_rhou[i] / Math.Max(rho, 1e-12f));
+                float p = PressureScalar(i);
+                float c = MathF.Sqrt(_gamma * p / Math.Max(rho, 1e-12f));
+                float local = u + c;
+                if (local > maxW) maxW = local;
+            }
+            maxW = Math.Max(maxW, 1e-8f);
+            return Math.Max(1, (int)Math.Ceiling((float)dtGlobal * maxW / ((float)cflTarget * _dx)));
+        }
+
+        // ==================== MAIN SIMULATION ==================================
+        public void SimulateSingleStep(double dtSub)
+        {
+            float dt = (float)dtSub;
+            int n = _n;
+
+            // --- 1. Left boundary face (index 0) – scalar -----------------------
+            float rhoL = _rho[0];
+            float uL = _rhou[0] / Math.Max(rhoL, 1e-12f);
+            float pL = PressureScalar(0);
+            ComputeLeftBoundaryFlux(rhoL, uL, pL, out _fluxM[0], out _fluxP[0], out _fluxE[0]);
+
+            // --- 2. Internal faces (1 .. n-1) – SIMD ---------------------------
+            int vectorSize = Vector<float>.Count;
+            int lastSimdFace = n - vectorSize;   // highest face index that starts a full vector block
+            for (int f = 1; f <= lastSimdFace; f += vectorSize)
+            {
+                SimdInternalFluxBlock(f, vectorSize);
+            }
+            // Scalar remainder for faces f .. n-1
+            for (int f = Math.Max(1, lastSimdFace + 1); f < n; f++)
+            {
+                float rhoLi = _rho[f - 1];
+                float uLi = _rhou[f - 1] / Math.Max(rhoLi, 1e-12f);
+                float pLi = PressureScalar(f - 1);
+                float rhoRi = _rho[f];
+                float uRi = _rhou[f] / Math.Max(rhoRi, 1e-12f);
+                float pRi = PressureScalar(f);
+                HLLCFluxScalar(rhoLi, uLi, pLi, rhoRi, uRi, pRi,
+                               out _fluxM[f], out _fluxP[f], out _fluxE[f]);
+            }
+
+            // --- 3. Right boundary face (index n) – scalar --------------------
+            float rhoR = _rho[n - 1];
+            float uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12f);
+            float pR = PressureScalar(n - 1);
+            ComputeRightBoundaryFlux(rhoR, uR, pR, out _fluxM[n], out _fluxP[n], out _fluxE[n]);
+
+            // --- 4. Cell update + damping + energy loss – SIMD -----------------
+            SimdCellUpdate(dt);
+        }
+
+        // ==================== PRIVATE SCALAR HELPERS ===========================
+        private float PressureScalar(int i)
+        {
+            float rho = Math.Max(_rho[i], 1e-12f);
+            return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
+        }
+
+        private void ComputeLeftBoundaryFlux(float rhoInt, float uInt, float pInt,
+            out float fm, out float fp, out float fe)
+        {
+            if (_aBCType == BoundaryType.GhostCell && _ghostLSet)
+                HLLCFluxScalar(_rhoGhostL, _uGhostL, _pGhostL, rhoInt, uInt, pInt, out fm, out fp, out fe);
+            else if (_aBCType == BoundaryType.OpenEnd)
+                OpenEndFluxLeft(rhoInt, uInt, pInt, _aAmbientPressure, out fm, out fp, out fe);
+            else if (_aBCType == BoundaryType.ZeroPressureOpen)
+            {
+                float rhoFace = rhoInt;
+                float uFace = uInt;
+                float pFace = _aAmbientPressure;
+                HLLCFluxScalar(rhoFace, uFace, pFace, rhoInt, uInt, pInt, out fm, out fp, out fe);
+            }
+            else if (_aBCType == BoundaryType.ClosedEnd)
+                ClosedEndFlux(rhoInt, uInt, pInt, false, out fm, out fp, out fe);
+            else
+            { fm = 0; fp = pInt; fe = 0; }
+        }
+
+        private void ComputeRightBoundaryFlux(float rhoInt, float uInt, float pInt,
+            out float fm, out float fp, out float fe)
+        {
+            if (_bBCType == BoundaryType.GhostCell && _ghostRSet)
+                HLLCFluxScalar(rhoInt, uInt, pInt, _rhoGhostR, _uGhostR, _pGhostR, out fm, out fp, out fe);
+            else if (_bBCType == BoundaryType.OpenEnd)
+                OpenEndFluxRight(rhoInt, uInt, pInt, _bAmbientPressure, out fm, out fp, out fe);
+            else if (_bBCType == BoundaryType.ZeroPressureOpen)
+            {
+                float rhoFace = rhoInt;
+                float uFace = uInt;
+                float pFace = _bAmbientPressure;
+                HLLCFluxScalar(rhoInt, uInt, pInt, rhoFace, uFace, pFace, out fm, out fp, out fe);
+            }
+            else if (_bBCType == BoundaryType.ClosedEnd)
+                ClosedEndFlux(rhoInt, uInt, pInt, true, out fm, out fp, out fe);
+            else
+            { fm = 0; fp = pInt; fe = 0; }
+        }
+
+        // ==================== SCALAR HLLC & BOUNDARY FLUX ======================
+        private void HLLCFluxScalar(float rL, float uL, float pL, float rR, float uR, float pR,
+                                    out float fm, out float fp, out float fe)
+        {
+            float cL = MathF.Sqrt(_gamma * pL / Math.Max(rL, 1e-12f));
+            float cR = MathF.Sqrt(_gamma * pR / Math.Max(rR, 1e-12f));
+            float EL = pL / ((_gamma - 1f) * rL) + 0.5f * uL * uL;
+            float ER = pR / ((_gamma - 1f) * rR) + 0.5f * uR * uR;
+            float SL = Math.Min(uL - cL, uR - cR);
+            float SR = Math.Max(uL + cL, uR + cR);
+
+            float denom = rL * (SL - uL) - rR * (SR - uR);
+            float Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
+
+            float FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
+            float FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
+
+            if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
+            else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
+            else if (Ss >= 0)
+            {
+                float rsL = rL * (SL - uL) / (SL - Ss);
+                float ps = pL + rL * (SL - uL) * (Ss - uL);
+                float EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
+                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
+            }
+            else
+            {
+                float rsR = rR * (SR - uR) / (SR - Ss);
+                float ps = pR + rR * (SR - uR) * (Ss - uR);
+                float EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
+                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
+            }
+        }
+
+        private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb,
+                                     out float fm, out float fp, out float fe)
+        {
+            float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
+            if (uInt <= -cInt) // supersonic inflow
+            {
+                fm = rhoInt * uInt;
+                fp = rhoInt * uInt * uInt + pInt;
+                fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
+                return;
+            }
+            if (uInt <= 0) // subsonic inflow
+            {
+                float T0 = 300f, R = 287f;
+                float ghostRho = pAmb / (R * T0);
+                HLLCFluxScalar(ghostRho, 0f, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
+                return;
+            }
+            // subsonic outflow
+            float s = pInt / MathF.Pow(rhoInt, _gamma);
+            float ghostRho2 = MathF.Pow(pAmb / s, 1f / _gamma);
+            float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho2);
+            float J_minus = uInt - 2f * cInt / (_gamma - 1f);
+            float uGhost = J_minus + 2f * cGhost / (_gamma - 1f);
+            if (uGhost < 0) uGhost = 0;
+            HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
+        }
+
+        private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb,
+                                      out float fm, out float fp, out float fe)
+        {
+            float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
+            if (uInt >= cInt) // supersonic outflow
+            {
+                fm = rhoInt * uInt;
+                fp = rhoInt * uInt * uInt + pInt;
+                fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
+                return;
+            }
+            if (uInt >= 0) // subsonic outflow
+            {
+                float s = pInt / MathF.Pow(rhoInt, _gamma);
+                float ghostRho = MathF.Pow(pAmb / s, 1f / _gamma);
+                float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho);
+                float J_plus = uInt + 2f * cInt / (_gamma - 1f);
+                float uGhost = J_plus - 2f * cGhost / (_gamma - 1f);
+                if (uGhost > 0) uGhost = 0;
+                HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe);
+                return;
+            }
+            // subsonic inflow
+            float T0 = 300f, R = 287f;
+            float ghostRho2 = pAmb / (R * T0);
+            HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho2, 0f, pAmb, out fm, out fp, out fe);
+        }
+
+        private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight,
+                                   out float fm, out float fp, out float fe)
+        {
+            float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt;
+            if (isRight)
+                HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
+            else
+                HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
+        }
+
+        // ==================== SIMD INTERNAL FACE ROUTINE ========================
+        private void SimdInternalFluxBlock(int startFace, int count)
+        {
+            int leftIdx = startFace - 1;
+            int rightIdx = startFace;
+
+            Vector<float> rL = new Vector<float>(_rho, leftIdx);
+            Vector<float> ruL = new Vector<float>(_rhou, leftIdx);
+            Vector<float> EL = new Vector<float>(_E, leftIdx);
+
+            Vector<float> rR = new Vector<float>(_rho, rightIdx);
+            Vector<float> ruR = new Vector<float>(_rhou, rightIdx);
+            Vector<float> ER = new Vector<float>(_E, rightIdx);
+
+            Vector<float> uL = ruL / rL;
+            Vector<float> uR = ruR / rR;
+
+            Vector<float> half = new Vector<float>(0.5f);
+            Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
+            Vector<float> gammaVec = new Vector<float>(_gamma);
+
+            Vector<float> pL = gammaMinus1 * (EL - half * ruL * ruL / rL);
+            Vector<float> pR = gammaMinus1 * (ER - half * ruR * ruR / rR);
+
+            Vector<float> cL = Vector.SquareRoot(gammaVec * pL / rL);
+            Vector<float> cR = Vector.SquareRoot(gammaVec * pR / rR);
+
+            Vector<float> SL = Vector.Min(uL - cL, uR - cR);
+            Vector<float> SR = Vector.Max(uL + cL, uR + cR);
+
+            Vector<float> num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR);
+            Vector<float> den = rL * (SL - uL) - rR * (SR - uR);
+            Vector<float> Ss = num / den;
+
+            Vector<float> eL = EL / rL;
+            Vector<float> eR = ER / rR;
+
+            // Left flux
+            Vector<float> Fm_L = ruL;
+            Vector<float> Fp_L = ruL * uL + pL;
+            Vector<float> Fe_L = (EL + pL) * uL;
+
+            // Right flux
+            Vector<float> Fm_R = ruR;
+            Vector<float> Fp_R = ruR * uR + pR;
+            Vector<float> Fe_R = (ER + pR) * uR;
+
+            // Star‑left fluxes
+            Vector<float> diffL = SL - uL;
+            Vector<float> dL_den = SL - Ss;
+            Vector<float> rsL = rL * diffL / dL_den;
+            Vector<float> psSL = pL + rL * diffL * (Ss - uL);
+            Vector<float> EsL = eL + (Ss - uL) * (Ss + pL / (rL * diffL));
+            Vector<float> Fm_starL = rsL * Ss;
+            Vector<float> Fp_starL = rsL * Ss * Ss + psSL;
+            Vector<float> Fe_starL = (rsL * EsL + psSL) * Ss;
+
+            // Star‑right fluxes
+            Vector<float> diffR = SR - uR;
+            Vector<float> dR_den = SR - Ss;
+            Vector<float> rsR = rR * diffR / dR_den;
+            Vector<float> psSR = pR + rR * diffR * (Ss - uR);
+            Vector<float> EsR = eR + (Ss - uR) * (Ss + pR / (rR * diffR));
+            Vector<float> Fm_starR = rsR * Ss;
+            Vector<float> Fp_starR = rsR * Ss * Ss + psSR;
+            Vector<float> Fe_starR = (rsR * EsR + psSR) * Ss;
+
+            var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector<float>.Zero);
+            var maskSRle0 = Vector.LessThanOrEqual(SR, Vector<float>.Zero);
+            var maskMiddle = ~(maskSLge0 | maskSRle0);
+            var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector<float>.Zero);
+            var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector<float>.Zero);
+
+            Vector<float> fm = Vector.ConditionalSelect(maskSLge0, Fm_L,
+                              Vector.ConditionalSelect(maskSRle0, Fm_R,
+                              Vector.ConditionalSelect(maskStarL, Fm_starL,
+                              Vector.ConditionalSelect(maskStarR, Fm_starR, Vector<float>.Zero))));
+
+            Vector<float> fp = Vector.ConditionalSelect(maskSLge0, Fp_L,
+                              Vector.ConditionalSelect(maskSRle0, Fp_R,
+                              Vector.ConditionalSelect(maskStarL, Fp_starL,
+                              Vector.ConditionalSelect(maskStarR, Fp_starR, Vector<float>.Zero))));
+
+            Vector<float> fe = Vector.ConditionalSelect(maskSLge0, Fe_L,
+                              Vector.ConditionalSelect(maskSRle0, Fe_R,
+                              Vector.ConditionalSelect(maskStarL, Fe_starL,
+                              Vector.ConditionalSelect(maskStarR, Fe_starR, Vector<float>.Zero))));
+
+            fm.CopyTo(_fluxM, startFace);
+            fp.CopyTo(_fluxP, startFace);
+            fe.CopyTo(_fluxE, startFace);
+        }
+
+        // ==================== SIMD CELL UPDATE + DAMPING + ENERGY LOSS =========
+        private void SimdCellUpdate(float dt)
+        {
+            float dt_dx = dt / _dx;
+            Vector<float> vDtDx = new Vector<float>(dt_dx);
+            float coeff = _laminarCoeff * (float)DampingMultiplier;
+            Vector<float> vCoeff = new Vector<float>(coeff);
+            Vector<float> vDt = new Vector<float>(dt);
+
+            int vectorSize = Vector<float>.Count;
+            int n = _n;
+            int lastSimdCell = n - vectorSize;
+
+            // Pre‑defined constants used in clamping
+            Vector<float> half = new Vector<float>(0.5f);
+            Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
+            Vector<float> ambientEnergyVec = new Vector<float>(_ambientEnergyReference);
+            Vector<float> energyRelaxRateVec = new Vector<float>(EnergyRelaxationRate);
+
+            for (int i = 0; i <= lastSimdCell; i += vectorSize)
+            {
+                // Load conserved
+                Vector<float> r = new Vector<float>(_rho, i);
+                Vector<float> ru = new Vector<float>(_rhou, i);
+                Vector<float> E = new Vector<float>(_E, i);
+
+                // Load fluxes at faces i (left) and i+1 (right)
+                Vector<float> flxM_L = new Vector<float>(_fluxM, i);
+                Vector<float> flxM_R = new Vector<float>(_fluxM, i + 1);
+                Vector<float> flxP_L = new Vector<float>(_fluxP, i);
+                Vector<float> flxP_R = new Vector<float>(_fluxP, i + 1);
+                Vector<float> flxE_L = new Vector<float>(_fluxE, i);
+                Vector<float> flxE_R = new Vector<float>(_fluxE, i + 1);
+
+                // Update conserved: Q_new = Q - dt/dx * (flux_right - flux_left)
+                Vector<float> newR = r - vDtDx * (flxM_R - flxM_L);
+                Vector<float> newRu = ru - vDtDx * (flxP_R - flxP_L);
+                Vector<float> newE = E - vDtDx * (flxE_R - flxE_L);
+
+                // Damping
+                Vector<float> dampingFactor = Vector.Exp(-vCoeff / r * vDt);
+                newRu *= dampingFactor;
+
+                // Energy loss (Newton cooling toward ambient)
+                Vector<float> relaxFactor = Vector.Exp(-energyRelaxRateVec * vDt);
+                newE = ambientEnergyVec + (newE - ambientEnergyVec) * relaxFactor;
+
+                // Clamp density
+                newR = Vector.Max(newR, new Vector<float>(1e-12f));
+
+                // Enforce minimal pressure (100 Pa) -> minimal energy
+                Vector<float> kinE = half * newRu * newRu / newR;
+                Vector<float> eMin = new Vector<float>(100f) / gammaMinus1 + kinE;
+                newE = Vector.Max(newE, eMin);
+
+                newR.CopyTo(_rho, i);
+                newRu.CopyTo(_rhou, i);
+                newE.CopyTo(_E, i);
+            }
+
+            // Scalar remainder
+            float relaxRate = EnergyRelaxationRate;
+            for (int i = Math.Max(0, lastSimdCell + 1); i < n; i++)
+            {
+                float r = _rho[i];
+                float ru = _rhou[i];
+                float E = _E[i];
+
+                float dM = _fluxM[i + 1] - _fluxM[i];
+                float dP = _fluxP[i + 1] - _fluxP[i];
+                float dE_flux = _fluxE[i + 1] - _fluxE[i];
+
+                float newR = r - dt_dx * dM;
+                float newRu = ru - dt_dx * dP;
+                float newE = E - dt_dx * dE_flux;
+
+                // Damping
+                float dampingFactor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
+                newRu *= dampingFactor;
+
+                // Energy loss
+                float relaxFactor = MathF.Exp(-relaxRate * dt);
+                newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
+
+                // Clamps
+                if (newR < 1e-12f) newR = 1e-12f;
+                float kin = 0.5f * newRu * newRu / newR;
+                float emin = 100f / (_gamma - 1f) + kin;
+                if (newE < emin) newE = emin;
+
+                _rho[i] = newR;
+                _rhou[i] = newRu;
+                _E[i] = newE;
+            }
         }
     }
 }

Components/Volume0D.cs

@@ -1,28 +1,27 @@
-using System;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
+using System;
 
 namespace FluidSim.Components
 {
     public class Volume0D
     {
-        public Port Port { get; private set; }
-
-        public double Mass { get; private set; }
-        public double InternalEnergy { get; private set; }
+        public double Mass { get; set; }
+        public double InternalEnergy { get; set; }
 
         public double Gamma { get; set; } = 1.4;
         public double GasConstant { get; set; } = 287.0;
 
         public double Volume { get; set; }
-        public double dVdt { get; set; }
+        public double Dvdt { get; set; }
 
         private double _dt;
 
-        public double Density => Mass / Volume;
-        public double Pressure => (Gamma - 1.0) * InternalEnergy / Volume;
-        public double Temperature => Pressure / (Density * GasConstant);
-        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
+        public double Density => Mass / Math.Max(Volume, 1e-12);
+        public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12);
+        public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
+        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12);
+
+        public double MassFlowRateIn { get; set; }
+        public double SpecificEnthalpyIn { get; set; }
 
         public Volume0D(double initialVolume, double initialPressure,
                         double initialTemperature, int sampleRate,
@@ -31,54 +30,40 @@ namespace FluidSim.Components
             GasConstant = gasConstant;
             Gamma = gamma;
             Volume = initialVolume;
-            dVdt = 0.0;
+            Dvdt = 0.0;
             _dt = 1.0 / sampleRate;
 
             double rho0 = initialPressure / (GasConstant * initialTemperature);
             Mass = rho0 * Volume;
             InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
-
-            Port = new Port();
-            PushStateToPort();
         }
 
-        public void PushStateToPort()
-        {
-            Port.Pressure = Pressure;
-            Port.Density = Density;
-            Port.Temperature = Temperature;
-            Port.SpecificEnthalpy = SpecificEnthalpy;
-        }
-
-        // Original integrate (uses the constructor’s sample rate)
-        public void Integrate()
-        {
-            Integrate(_dt);
-        }
-
-        public void SetPressure(double newPressure)
-        {
-            InternalEnergy = newPressure * Volume / (Gamma - 1.0);
-            // Mass stays the same, so density is unchanged
-        }
-
-        // New overload: integrate with a custom time step (for sub‑steps)
         public void Integrate(double dtOverride)
         {
-            double mdot = Port.MassFlowRate;
-            double h_in = Port.SpecificEnthalpy;
-
-            double dm = mdot * dtOverride;
-            double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride;
+            double dm = MassFlowRateIn * dtOverride;
+            double dE = (MassFlowRateIn * SpecificEnthalpyIn) * dtOverride - Pressure * Dvdt * dtOverride;
 
             Mass += dm;
             InternalEnergy += dE;
 
-            // Hard physical bounds – prevent NaN and unphysical states
-            if (Mass < 1e-12) Mass = 1e-12;
-            if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
+            // ---- ABSOLUTE SAFEGUARD: keep at least 1 µg of gas at ambient pressure ----
+            double minMass = 1e-9;
+            double V = Math.Max(Volume, 1e-12);
+            if (Mass < minMass)
+            {
+                Mass = minMass;
+                InternalEnergy = 5000.0 * V / (Gamma - 1.0);   // 0.05 bar, not ambient
+            }
+            else if (InternalEnergy < 0.0)
+            {
+                InternalEnergy = 101325.0 * V / (Gamma - 1.0);
+            }
 
-            PushStateToPort();
+            // Final non‑negative clamp
+            if (Mass < 0.0) Mass = 0.0;
+            if (InternalEnergy < 0.0) InternalEnergy = 0.0;
         }
+
+        public void Integrate() => Integrate(_dt);
     }
 }

Core/Constants.cs

@@ -0,0 +1,11 @@
+namespace FluidSim.Core
+{
+    public static class Constants
+    {
+        public const double Gamma   = 1.4;
+        public const double R_gas   = 287.0;         // J/(kg·K)
+        public const double P_amb   = 101325.0;      // Pa
+        public const double T_amb   = 300.0;         // K
+        public static readonly double Rho_amb = P_amb / (R_gas * T_amb);  // ≈ 1.177 kg/m³
+    }
+}

Core/NozzleFlow.cs

@@ -0,0 +1,72 @@
+using System;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+    public static class NozzleFlow
+    {
+        public static void Compute(Volume0D vol, double area, double downstreamPressure,
+                                   out double massFlow, out double rhoFace, out double uFace, out double pFace,
+                                   double gamma = 1.4)
+        {
+            massFlow = 0.0;
+            rhoFace = 0.0;
+            uFace = 0.0;
+            pFace = 0.0;
+
+            if (vol == null || vol.Mass <= 0 || vol.Volume <= 0)
+                return;
+
+            double p0 = vol.Pressure;
+            double T0 = vol.Temperature;
+            double R = vol.GasConstant;
+            double rho0 = vol.Density;
+
+            if (double.IsNaN(p0) || double.IsNaN(T0) || double.IsNaN(rho0) ||
+                p0 <= 0 || T0 <= 0 || rho0 <= 0)
+                return;
+
+            double pr = downstreamPressure / p0;
+            double choked = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
+
+            // If pr > 1, flow is INTO the cylinder (reverse), so we swap the roles.
+            bool reverse = (pr > 1.0);
+            if (reverse)
+            {
+                // Treat the cylinder as the downstream, the pipe as the upstream.
+                double p_up = downstreamPressure;
+                double T_up = 300.0;        // pipe temperature (ambient)
+                double rho_up = downstreamPressure / (R * T_up);
+
+                double pr_rev = p0 / p_up;  // now cylinder / pipe
+                if (pr_rev < choked) pr_rev = choked;
+
+                double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr_rev, -(gamma - 1.0) / gamma) - 1.0));
+                if (double.IsNaN(M)) return;
+
+                // Flow from pipe INTO cylinder (positive mass flow into volume)
+                uFace = M * Math.Sqrt(gamma * R * T_up);
+                rhoFace = rho_up * Math.Pow(pr_rev, 1.0 / gamma);
+                pFace = p_up * pr_rev;
+                massFlow = rhoFace * uFace * area;
+                // massFlow is positive = into cylinder
+            }
+            else
+            {
+                // Normal flow out of cylinder
+                if (pr < choked) pr = choked;
+
+                double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
+                if (double.IsNaN(M)) return;
+
+                uFace = M * Math.Sqrt(gamma * R * T0);
+                rhoFace = rho0 * Math.Pow(pr, 1.0 / gamma);
+                pFace = p0 * pr;
+                massFlow = -rhoFace * uFace * area;   // negative = out of cylinder
+            }
+
+            if (double.IsNaN(massFlow) || double.IsInfinity(massFlow))
+                massFlow = 0.0;
+        }
+    }
+}

Core/OutdoorExhaustReverb.cs

@@ -0,0 +1,187 @@
+using System;
+
+namespace FluidSim.Core
+{
+    public class OutdoorExhaustReverb
+    {
+        // ========== Early reflection delays (stereo: left/right) ==========
+        private readonly DelayLine groundL, groundR;
+        private readonly DelayLine wall1L, wall1R;
+        private readonly DelayLine wall2L, wall2R;
+
+        // ========== Diffuse tail FDNs (left/right each with 8 channels) ==========
+        private const int FDN_CHANNELS = 8;
+        private readonly DelayLine[] fdnL, fdnR;
+        private readonly float[] stateL, stateR;
+        private readonly OrthonormalMixer mixerL, mixerR;
+        private readonly LowPassFilter[] filterL, filterR;
+
+        public float DryMix      { get; set; } = 1.0f;
+        public float EarlyMix    { get; set; } = 0.5f;
+        public float TailMix     { get; set; } = 0.9f;
+        public float Feedback    { get; set; } = 0.75f;   // safe range 0.7‑0.9
+        public float DampingFreq { get; set; } = 6000f;    // Hz
+
+        public OutdoorExhaustReverb(int sampleRate)
+        {
+            // Early reflections – left/right offset by ~1‑2 ms for stereo width
+            groundL = new DelayLine((int)(sampleRate * 0.008));   // 8 ms
+            groundR = new DelayLine((int)(sampleRate * 0.010));   // 10 ms
+            wall1L  = new DelayLine((int)(sampleRate * 0.045));
+            wall1R  = new DelayLine((int)(sampleRate * 0.047));
+            wall2L  = new DelayLine((int)(sampleRate * 0.080));
+            wall2R  = new DelayLine((int)(sampleRate * 0.082));
+
+            // FDN delay lengths – prime numbers, offset between L/R
+            int[] lengthsL = { 3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009 };
+            int[] lengthsR = { 3217, 4049, 5531, 7043, 8537, 10037, 11519, 13033 };
+            fdnL = new DelayLine[FDN_CHANNELS];
+            fdnR = new DelayLine[FDN_CHANNELS];
+            for (int i = 0; i < FDN_CHANNELS; i++)
+            {
+                int lenL = Math.Min(lengthsL[i], (int)(sampleRate * 0.25));
+                int lenR = Math.Min(lengthsR[i], (int)(sampleRate * 0.25));
+                fdnL[i] = new DelayLine(lenL);
+                fdnR[i] = new DelayLine(lenR);
+            }
+
+            stateL = new float[FDN_CHANNELS];
+            stateR = new float[FDN_CHANNELS];
+            mixerL = new OrthonormalMixer(FDN_CHANNELS);
+            mixerR = new OrthonormalMixer(FDN_CHANNELS);
+
+            filterL = new LowPassFilter[FDN_CHANNELS];
+            filterR = new LowPassFilter[FDN_CHANNELS];
+            for (int i = 0; i < FDN_CHANNELS; i++)
+            {
+                filterL[i] = new LowPassFilter(sampleRate, DampingFreq);
+                filterR[i] = new LowPassFilter(sampleRate, DampingFreq);
+            }
+        }
+
+        /// <summary>Stereo reverb – returns (left, right) sample pair.</summary>
+        public (float left, float right) ProcessStereo(float drySample)
+        {
+            // ---- Early reflections ----
+            float gL = groundL.ReadWrite(drySample * 0.8f);
+            float gR = groundR.ReadWrite(drySample * 0.8f);
+            float w1L = wall1L.ReadWrite(drySample * 0.5f);
+            float w1R = wall1R.ReadWrite(drySample * 0.5f);
+            float w2L = wall2L.ReadWrite(drySample * 0.4f);
+            float w2R = wall2R.ReadWrite(drySample * 0.4f);
+
+            float earlyL = (gL + w1L + w2L) * EarlyMix;
+            float earlyR = (gR + w1R + w2R) * EarlyMix;
+
+            // ---- Read diffuse tail ----
+            float[] delOutL = new float[FDN_CHANNELS];
+            float[] delOutR = new float[FDN_CHANNELS];
+            for (int i = 0; i < FDN_CHANNELS; i++)
+            {
+                delOutL[i] = fdnL[i].Read();
+                delOutR[i] = fdnR[i].Read();
+            }
+
+            // Mix via orthonormal matrix
+            float[] mixL = new float[FDN_CHANNELS];
+            float[] mixR = new float[FDN_CHANNELS];
+            mixerL.Process(delOutL, mixL);
+            mixerR.Process(delOutR, mixR);
+
+            // Feedback + air absorption
+            for (int i = 0; i < FDN_CHANNELS; i++)
+            {
+                stateL[i] = drySample * 0.15f + Feedback * mixL[i];
+                stateL[i] = filterL[i].Process(stateL[i]);
+                fdnL[i].Write(stateL[i]);
+
+                stateR[i] = drySample * 0.15f + Feedback * mixR[i];
+                stateR[i] = filterR[i].Process(stateR[i]);
+                fdnR[i].Write(stateR[i]);
+            }
+
+            float tailL = 0.0f, tailR = 0.0f;
+            for (int i = 0; i < FDN_CHANNELS; i++)
+            {
+                tailL += delOutL[i];
+                tailR += delOutR[i];
+            }
+            tailL *= TailMix;
+            tailR *= TailMix;
+
+            float left  = drySample * DryMix + earlyL + tailL;
+            float right = drySample * DryMix + earlyR + tailR;
+            return (left, right);
+        }
+
+        /// <summary>Mono fallback – sums left+right / 2.</summary>
+        public float Process(float drySample)
+        {
+            var (l, r) = ProcessStereo(drySample);
+            return (l + r) * 0.5f;
+        }
+
+        // ========== Helper classes ==========
+        private class DelayLine
+        {
+            private float[] buffer;
+            private int writePos;
+            public DelayLine(int length)
+            {
+                buffer = new float[Math.Max(length, 1)];
+            }
+            public float Read() => buffer[writePos];
+            public void Write(float value)
+            {
+                buffer[writePos] = value;
+                writePos = (writePos + 1) % buffer.Length;
+            }
+            public float ReadWrite(float value)
+            {
+                float outVal = buffer[writePos];
+                buffer[writePos] = value;
+                writePos = (writePos + 1) % buffer.Length;
+                return outVal;
+            }
+        }
+
+        private class LowPassFilter
+        {
+            private float b0, a1, y1;
+            private float sampleRate;
+            public LowPassFilter(int sampleRate, float cutoff)
+            {
+                this.sampleRate = sampleRate;
+                SetCutoff(cutoff);
+            }
+            public void SetCutoff(float cutoff)
+            {
+                float w = 2 * (float)Math.PI * cutoff / sampleRate;
+                float a0 = 1 + w;
+                b0 = w / a0;
+                a1 = (1 - w) / a0;
+            }
+            public float Process(float x)
+            {
+                float y = b0 * x - a1 * y1;
+                y1 = y;
+                return y;
+            }
+        }
+
+        private class OrthonormalMixer
+        {
+            private int size;
+            public OrthonormalMixer(int size) => this.size = size;
+
+            public void Process(float[] input, float[] output)
+            {
+                float sum = 0;
+                for (int i = 0; i < size; i++) sum += input[i];
+                float factor = 2.0f / size;
+                for (int i = 0; i < size; i++)
+                    output[i] = factor * sum - input[i];
+            }
+        }
+    }
+}

Core/PipeVolumeConnection.cs

@@ -0,0 +1,22 @@
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+    public class PipeVolumeConnection
+    {
+        public Volume0D Volume { get; }
+        public Pipe1D Pipe { get; }
+        public bool IsPipeLeftEnd { get; }
+        public double OrificeArea { get; set; }
+
+        public double LastMassFlowIntoVolume { get; set; }
+
+        public PipeVolumeConnection(Volume0D vol, Pipe1D pipe, bool isPipeLeftEnd, double orificeArea)
+        {
+            Volume = vol;
+            Pipe = pipe;
+            IsPipeLeftEnd = isPipeLeftEnd;
+            OrificeArea = orificeArea;
+        }
+    }
+}

Core/Solver.cs

@@ -1,7 +1,6 @@
 using System;
 using System.Collections.Generic;
 using FluidSim.Components;
-using FluidSim.Interfaces;
 
 namespace FluidSim.Core
 {
@@ -9,162 +8,113 @@ namespace FluidSim.Core
     {
         private readonly List<Volume0D> _volumes = new();
         private readonly List<Pipe1D> _pipes = new();
-        private readonly List<Connection> _connections = new();
+        private readonly List<PipeVolumeConnection> _connections = new();
 
         private double _dt;
+        private double _ambientPressure = 101325.0;
 
+        public void SetAmbientPressure(double p) => _ambientPressure = p;
         public void AddVolume(Volume0D v) => _volumes.Add(v);
         public void AddPipe(Pipe1D p) => _pipes.Add(p);
-        public void AddConnection(Connection c) => _connections.Add(c);
+        public void AddConnection(PipeVolumeConnection c) => _connections.Add(c);
         public void SetTimeStep(double dt) => _dt = dt;
 
-        /// <summary>
-        /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right).
-        /// </summary>
         public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0)
         {
             if (isA)
             {
                 pipe.SetABoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetAAmbientPressure(ambientPressure);
+                if (type == BoundaryType.OpenEnd) pipe.SetAAmbientPressure(ambientPressure);
             }
             else
             {
                 pipe.SetBBoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetBAmbientPressure(ambientPressure);
+                if (type == BoundaryType.OpenEnd) pipe.SetBAmbientPressure(ambientPressure);
             }
         }
 
         public float Step()
         {
-            // 1. Volumes publish state
-            foreach (var v in _volumes)
-                v.PushStateToPort();
-
-            // 2. Set volume BCs for volume‑coupled ends
+            // 1. For each connection, handle flow or closed wall
             foreach (var conn in _connections)
             {
-                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+                double area = conn.OrificeArea;
+                if (area < 1e-12)   // valve closed → treat as solid wall
                 {
-                    var pipe = GetPipe(conn.PortA);
-                    bool isA = pipe.PortA == conn.PortA;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortA, conn.PortB);
+                    conn.Volume.MassFlowRateIn = 0.0;
+                    conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy; // not used
+
+                    // Set ghost to a reflective wall (u = -u_pipe, same p, ρ)
+                    int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1;
+                    double rho = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6);
+                    double p   = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0);
+                    double u   = conn.Pipe.GetCellVelocity(cellIdx);
+                    if (conn.IsPipeLeftEnd)
+                        conn.Pipe.SetGhostLeft(rho, -u, p);
+                    else
+                        conn.Pipe.SetGhostRight(rho, -u, p);
+                    continue;
                 }
-                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+
+                // Valve open → use the nozzle model
+                double downstreamPressure = conn.IsPipeLeftEnd
+                    ? conn.Pipe.GetCellPressure(0)
+                    : conn.Pipe.GetCellPressure(conn.Pipe.GetCellCount() - 1);
+
+                NozzleFlow.Compute(conn.Volume, area, downstreamPressure,
+                                out double mdot, out double rhoFace, out double uFace, out double pFace,
+                                gamma: conn.Volume.Gamma);
+
+                // Clamp mdot to available mass
+                double maxMdot = conn.Volume.Mass / _dt;
+                conn.LastMassFlowIntoVolume = mdot;
+                if (mdot > maxMdot) mdot = maxMdot;
+                if (mdot < -maxMdot) mdot = -maxMdot;
+
+                conn.Volume.MassFlowRateIn = mdot;
+                // enthalpy: if inflow, use pipe enthalpy; if outflow, use cylinder enthalpy
+                if (mdot >= 0)
                 {
-                    var pipe = GetPipe(conn.PortB);
-                    bool isA = pipe.PortB == conn.PortB;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortB, conn.PortA);
+                    int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1;
+                    double pPipe = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0);
+                    double rhoPipe = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6);
+                    conn.Volume.SpecificEnthalpyIn = (conn.Volume.Gamma / (conn.Volume.Gamma - 1.0)) * pPipe / rhoPipe;
                 }
+                else
+                {
+                    conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy;
+                }
+
+                // Integrate the volume
+                conn.Volume.Integrate(_dt);
+
+                // Set ghost from nozzle face state (but don't allow zero density)
+                if (rhoFace < 1e-6) rhoFace = Constants.Rho_amb;
+                if (pFace < 100.0)  pFace = Constants.P_amb;
+                if (conn.IsPipeLeftEnd)
+                    conn.Pipe.SetGhostLeft(rhoFace, uFace, pFace);
+                else
+                    conn.Pipe.SetGhostRight(rhoFace, uFace, pFace);
             }
 
-            // 3. Sub‑steps
+            // 2. Sub‑step pipes
             int nSub = 1;
             foreach (var p in _pipes)
                 nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
             double dtSub = _dt / nSub;
 
             for (int sub = 0; sub < nSub; sub++)
-            {
                 foreach (var p in _pipes)
                     p.SimulateSingleStep(dtSub);
 
-                foreach (var conn in _connections)
-                {
-                    if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortA);
-                        bool isA = pipe.PortA == conn.PortA;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortA, conn.PortB, dtSub);
-                    }
-                    else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortB);
-                        bool isA = pipe.PortB == conn.PortB;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
-                    }
-                }
-
-                if (sub < nSub - 1)
-                {
-                    foreach (var v in _volumes)
-                        v.PushStateToPort();
-
-                    foreach (var conn in _connections)
-                    {
-                        if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortA);
-                            bool isA = pipe.PortA == conn.PortA;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortA, conn.PortB);
-                        }
-                        else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortB);
-                            bool isA = pipe.PortB == conn.PortB;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortB, conn.PortA);
-                        }
-                    }
-                }
-            }
-
-            // 5. Audio samples (none for now, but placeholder)
-            var audioSamples = new List<float>();
-            foreach (var conn in _connections)
-            {
-                if (conn is SoundConnection sc)
-                    audioSamples.Add(sc.GetAudioSample());
-            }
-
-            // 6. Clear BC flags
+            // 3. Clear ghost flags
             foreach (var p in _pipes)
-                p.ClearBC();
+                p.ClearGhostFlag();
 
-            return SoundProcessor.MixAndClip(audioSamples.ToArray());
-        }
-
-        private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
-        private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
-        private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
-        private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p);
-
-        private void SetVolumeBC(Port pipePort, Port volPort)
-        {
-            var pipe = GetPipe(pipePort);
-            if (pipe == null) return;
-            bool isA = pipe.PortA == pipePort;
-            if (isA)
-                pipe.SetAVolumeState(volPort.Density, volPort.Pressure);
-            else
-                pipe.SetBVolumeState(volPort.Density, volPort.Pressure);
-        }
-
-        private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
-        {
-            double mdot = pipePort.MassFlowRate;
-            volPort.MassFlowRate = -mdot;
-
-            if (mdot < 0)  // pipe → volume
-            {
-                volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
-            }
-            // else volume’s own enthalpy (from PushStateToPort) is used
-
-            GetVolume(volPort)?.Integrate(dtSub);
+            // 4. Return exhaust tailpipe mass flow
+            if (_pipes.Count > 0)
+                return (float)_pipes[0].GetOpenEndMassFlow();
+            return 0f;
         }
     }
 }

Core/SoundProcessor.cs

@@ -1,23 +1,48 @@
-namespace FluidSim.Core
-{
-    /// <summary>
-    /// Mixes multiple audio samples and applies a soft‑clipping tanh.
-    /// </summary>
-    public static class SoundProcessor
-    {
-        /// <summary>Overall gain applied after mixing (before tanh).</summary>
-        public static float MasterGain { get; set; } = 0.01f;
+using System;
+using FluidSim.Interfaces;
 
-        /// <summary>
-        /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1].
-        /// </summary>
-        public static float MixAndClip(params float[] samples)
+namespace FluidSim.Core
+{
+    public class SoundProcessor
+    {
+        private readonly double dt;
+        private readonly double scaleFactor;  // 1 / (4π r) and a user gain
+        private double prevMassFlowOut;
+
+        // Simple low‑pass for derivative smoothing (≈ 2‑3 ms)
+        private double smoothDMdt;
+        private readonly double alpha;
+
+        public float Gain { get; set; } = 1.0f;
+
+        public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0)
         {
-            float sum = 0f;
-            foreach (float s in samples)
-                sum += s;
-            sum *= MasterGain;
-            return sum;
+            dt = 1.0 / sampleRate;
+            scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters);
+
+            // Smoothing time constant ~ 2 ms, blocks single‑sample spikes
+            double tau = 0.002;
+            alpha = Math.Exp(-dt / tau);
+        }
+
+        public float Process(Port port)
+        {
+            // Outflow mass flow (positive = leaving pipe)
+            double flowOut = -port.MassFlowRate;
+
+            // Derivative
+            double rawDerivative = (flowOut - prevMassFlowOut) / dt;
+            prevMassFlowOut = flowOut;
+
+            // Smooth the derivative to kill isolated spikes
+            smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative;
+
+            // Far‑field monopole pressure
+            double pressure = smoothDMdt * scaleFactor * Gain;
+
+            // Soft clip to ±1 for audio output (safe limit)
+            float sample = (float)Math.Tanh(pressure);
+            return sample;
         }
     }
 }

FluidSim.csproj

@@ -5,7 +5,7 @@
     <TargetFramework>net10.0</TargetFramework>
     <ImplicitUsings>enable</ImplicitUsings>
     <Nullable>enable</Nullable>
-    <PublishAot>true</PublishAot>
+    <PublishAot>false</PublishAot>
     <InvariantGlobalization>true</InvariantGlobalization>
   </PropertyGroup>
 

Program.cs

@@ -12,37 +12,37 @@ public class Program
     private const double DrawFrequency = 60.0;
     private static Scenario scenario;
 
-    // Speed control
-    //private static double desiredSpeed = 1.0;
-    private static double desiredSpeed = 0.0001;
+    // Speed control (existing + new throttle)
+    private static double desiredSpeed = 0.01;
     private static double currentSpeed = desiredSpeed;
     private const double MinSpeed = 0.0001;
     private const double MaxSpeed = 1.0;
     private const double ScrollFactor = 1.1;
 
-    // Space‑toggle state
-    private static double lastDesiredSpeed = 0.1;   // remembers the last non‑1.0 scroll speed
-    private static bool isRealTime = true;          // true when desiredSpeed == 1.0
+    private static double lastDesiredSpeed = 0.1;
+    private static bool isRealTime = false;
+
+    // Throttle smoothing
+    private static double targetThrottle = 0.0;   // 1.0 when W is pressed, 0.0 otherwise
+    private static double currentThrottle = 0.0;
+    private const double ThrottleSmoothing = 20.0; // rate of change
 
     private static volatile bool running = true;
 
     public static void Main()
     {
         var mode = new VideoMode(new Vector2u(1280, 720));
-        var window = new RenderWindow(mode, "Pipe Resonator");
+        var window = new RenderWindow(mode, "FluidSim - Engine (W = throttle)");
         window.SetVerticalSyncEnabled(true);
         window.Closed += (_, _) => { running = false; window.Close(); };
         window.MouseWheelScrolled += OnMouseWheel;
         window.KeyPressed += OnKeyPressed;
 
         var soundEngine = new SoundEngine(bufferCapacity: 16384);
-        soundEngine.Volume = 70;
+        soundEngine.Volume = 100;
         soundEngine.Start();
 
-        //scenario = new PipeResonatorScenario();
-        //scenario = new HelmholtzResonatorScenario();
-        scenario = new SodShockTubeScenario();
-
+        scenario = new EngineScenario();
         scenario.Initialize(SampleRate);
 
         var stopwatch = Stopwatch.StartNew();
@@ -50,17 +50,13 @@ public class Program
         double drawInterval = 1.0 / DrawFrequency;
         double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
 
-        // Resampling buffer
-        List<float> simBuffer = new List<float>(4096);
+        var simBuffer = new List<float>(4096);
         double readIndex = 0.0;
-
         for (int i = 0; i < 4; i++)
             simBuffer.Add(scenario.Process());
 
         long totalSimSteps = simBuffer.Count;
         long totalOutputSamples = 0;
-
-        double lastRealTime = stopwatch.Elapsed.TotalSeconds;
         const int outputChunk = 256;
         float[] outputBuf = new float[outputChunk];
 
@@ -69,17 +65,22 @@ public class Program
             window.DispatchEvents();
 
             double currentRealTime = stopwatch.Elapsed.TotalSeconds;
-            double dtSpeed = currentRealTime - lastSpeedUpdateTime;
+            double dtClock = currentRealTime - lastSpeedUpdateTime;
             lastSpeedUpdateTime = currentRealTime;
 
-            // Smoothly transition currentSpeed → desiredSpeed
-            // When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp
-            double smoothingRate = 8.0;   // higher = faster catch‑up
-            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
+            // Smooth simulation speed
+            double speedSmoothing = 8.0;
+            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-speedSmoothing * dtClock));
 
-            // ---------- Generate audio ----------
+            // ---- THROTTLE INPUT ----
+            targetThrottle = Keyboard.IsKeyPressed(Keyboard.Key.W) ? 1.0 : 0.0;
+            currentThrottle += (targetThrottle - currentThrottle) * (1.0 - Math.Exp(-ThrottleSmoothing * dtClock));
+            // Push to engine scenario (if it's an EngineScenario)
+            if (scenario is EngineScenario engine)
+                engine.Throttle = currentThrottle;
+
+            // Generate audio
             double targetAudioClock = currentRealTime + 0.05;
-
             while (totalOutputSamples < targetAudioClock * SampleRate && running)
             {
                 int toGenerate = (int)Math.Min(
@@ -101,11 +102,9 @@ public class Program
                     int i0 = (int)readIndex;
                     int i1 = i0 + 1;
                     double frac = readIndex - i0;
-
                     float y0 = simBuffer[Math.Clamp(i0, 0, simBuffer.Count - 1)];
                     float y1 = simBuffer[Math.Clamp(i1, 0, simBuffer.Count - 1)];
                     outputBuf[i] = (float)(y0 + (y1 - y0) * frac);
-
                     readIndex += currentSpeed;
 
                     while (readIndex >= 1.0 && simBuffer.Count > 2)
@@ -117,23 +116,23 @@ public class Program
 
                 int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
                 totalOutputSamples += accepted;
-
                 if (accepted < toGenerate)
                     break;
             }
 
-            // ---------- Drawing & title ----------
+            // Drawing & title
             if (currentRealTime - lastDrawTime >= drawInterval)
             {
                 double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate);
                 double simTime = totalSimSteps / (double)SampleRate;
                 string toggleHint = isRealTime ? "[Space] slow mo" : "[Space] real time";
+                string throttleHint = Keyboard.IsKeyPressed(Keyboard.Key.W) ? "W held" : "W released";
                 window.SetTitle(
-                    $"{toggleHint}  Sim: {simTime:F2}s | " +
-                    $"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
-                    $"Actual: {actualSpeed:F2}x"
+                    $"{toggleHint}  {throttleHint}  " +
+                    $"Thr: {currentThrottle:F2}  " +
+                    $"Speed: {currentSpeed:F3}x → {desiredSpeed:F3}x  " +
+                    $"Act: {actualSpeed:F2}x"
                 );
-
                 window.Clear(Color.Black);
                 scenario.Draw(window);
                 window.Display();
@@ -145,22 +144,18 @@ public class Program
         window.Dispose();
     }
 
+    // (Mouse wheel, space toggle unchanged)
     private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
     {
         bool wasRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-
         if (e.Delta > 0)
             desiredSpeed *= ScrollFactor;
         else if (e.Delta < 0)
             desiredSpeed /= ScrollFactor;
 
         desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed);
-
-        // Update the remembered slow-mo speed (unless we are exactly at 1.0)
         if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6)
             lastDesiredSpeed = desiredSpeed;
-
-        // Update isRealTime flag
         isRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
     }
 
@@ -169,15 +164,9 @@ public class Program
         if (e.Code == Keyboard.Key.Space)
         {
             if (isRealTime)
-            {
-                // Switch to the remembered slow speed
                 desiredSpeed = lastDesiredSpeed;
-            }
             else
-            {
-                // Switch back to real time
                 desiredSpeed = 1.0;
-            }
             isRealTime = !isRealTime;
         }
     }

Scenarios/EngineScenario.cs

@@ -0,0 +1,325 @@
+using System;
+using FluidSim.Components;
+using FluidSim.Utils;
+using FluidSim.Interfaces;
+using SFML.Graphics;
+using SFML.System;
+
+namespace FluidSim.Core
+{
+    public class EngineScenario : Scenario
+    {
+        private Solver solver;
+        private Crankshaft crankshaft;
+        private EngineCylinder engineCyl;
+        private Pipe1D exhaustPipe;
+        private Pipe1D intakePipe;
+        private PipeVolumeConnection couplingExhaust;
+        private PipeVolumeConnection couplingIntake;
+        private SoundProcessor exhaustSoundProcessor;
+        private SoundProcessor intakeSoundProcessor;
+        private OutdoorExhaustReverb reverb;
+
+        private Port exhaustPort = new Port();
+        private Port intakePort = new Port();
+
+        private double dt;
+        private double exhPipeArea, intPipeArea;
+        private const double AmbientPressure = 101325.0;
+        private double time;
+        private int stepCount = 0;
+        private const int LogInterval = 1000;
+
+        public double Throttle { get; set; } = 0.15;
+        private const double FullLoadPeakPressure = 140.0 * Units.bar;
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            // Crankshaft
+            crankshaft = new Crankshaft(initialRPM: 2000.0)
+            {
+                Inertia = 0.05,
+                FrictionConstant = 0.2,
+                FrictionViscous = 0.025
+            };
+
+            // Exhaust pipe (longer, larger)
+            double exhLength = 0.5;
+            double exhRadius = 1.5 * Units.cm;
+            exhPipeArea = Math.PI * exhRadius * exhRadius;
+            exhaustPipe = new Pipe1D(exhLength, exhPipeArea, sampleRate, forcedCellCount: 30);
+            exhaustPipe.SetUniformState(1.225, 0.0, AmbientPressure);
+            exhaustPipe.DampingMultiplier = 0.0;
+            exhaustPipe.EnergyRelaxationRate = 100.0f;
+
+            // Intake pipe (shorter, narrower)
+            double intLength = 0.1;
+            double intRadius = 1 * Units.cm;
+            intPipeArea = Math.PI * intRadius * intRadius;
+            intakePipe = new Pipe1D(intLength, intPipeArea, sampleRate, forcedCellCount: 10);
+            intakePipe.SetUniformState(1.225, 0.0, AmbientPressure);
+
+            // Cylinder (starts at BDC, fresh charge)
+            engineCyl = new EngineCylinder(crankshaft,
+                bore: 56 * Units.mm, stroke: 57 * Units.mm, compressionRatio: 9.5,
+                exhPipeArea: exhPipeArea, intPipeArea: intPipeArea, sampleRate: sampleRate);
+            engineCyl.ignition = true;
+
+            // Set crank to BDC (180°) and sync
+            crankshaft.CrankAngle = Math.PI;
+            crankshaft.PreviousAngle = Math.PI;   // make sure this property is settable (public setter)
+
+            // Couplings
+            couplingExhaust = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0);
+            couplingIntake  = new PipeVolumeConnection(engineCyl.Cylinder, intakePipe, false, orificeArea: 0.0);
+
+            // Solver
+            solver = new Solver();
+            solver.SetTimeStep(dt);
+            solver.AddVolume(engineCyl.Cylinder);
+            solver.AddPipe(exhaustPipe);
+            solver.AddPipe(intakePipe);
+            solver.AddConnection(couplingExhaust);
+            solver.AddConnection(couplingIntake);
+            solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, AmbientPressure);
+            solver.SetPipeBoundary(intakePipe, true,  BoundaryType.GhostCell);        // cylinder side – left
+            solver.SetPipeBoundary(intakePipe, false, BoundaryType.OpenEnd, AmbientPressure);  // ambient side – right
+
+            // Sound
+            exhaustSoundProcessor = new SoundProcessor(sampleRate, exhRadius * 2);
+            exhaustSoundProcessor.Gain = 0.001f;
+
+            intakeSoundProcessor = new SoundProcessor(sampleRate, intRadius * 2);
+            intakeSoundProcessor.Gain = 0.001f;
+
+            // Reverb
+            reverb = new OutdoorExhaustReverb(sampleRate);
+            reverb.DryMix = 1.0f;
+            reverb.EarlyMix = 0.5f;
+            reverb.TailMix = 0.9f;
+            reverb.Feedback = 0.9f;
+            reverb.DampingFreq = 6000f;
+
+            Console.WriteLine("=== Engine with intake & cycle‑aware valves ===");
+        }
+
+        public override float Process()
+        {
+            double idleThrottle = 0.1;
+            if (crankshaft.AngularVelocity < 80) idleThrottle = 0.2;
+            double throttle = Math.Clamp(Throttle, idleThrottle, 1.0);
+            double targetPressure = throttle * FullLoadPeakPressure;
+            engineCyl.TargetPeakPressure = targetPressure;
+
+            engineCyl.Step(dt);
+            crankshaft.Step(dt);
+
+            couplingExhaust.OrificeArea = engineCyl.ExhaustOrificeArea;
+            couplingIntake.OrificeArea  = engineCyl.IntakeOrificeArea;
+
+            solver.Step();
+
+            UpdateExhaustPort();
+            UpdateIntakePort();
+            float dryExhaust = exhaustSoundProcessor.Process(exhaustPort);
+            float dryIntake = intakeSoundProcessor.Process(intakePort);
+            float dry = dryExhaust + dryIntake;
+
+            float wet = reverb.Process(dry);
+
+            if (++stepCount % LogInterval == 0) Log();
+
+            return wet;
+        }
+
+        private void Log()
+        {
+            double rpm = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI);
+            double cycleDeg = (engineCyl.CycleAngle * 180.0 / Math.PI) % 720.0;
+            string stroke = cycleDeg < 180.0 ? "Power" :
+                            cycleDeg < 360.0 ? "Exhaust" :
+                            cycleDeg < 540.0 ? "Intake" : "Compression";
+
+            // Cylinder
+            double pCyl = engineCyl.Cylinder.Pressure;
+            double TCyl = engineCyl.Cylinder.Temperature;
+            double VCyl = engineCyl.Cylinder.Volume;
+            double mCyl = engineCyl.Cylinder.Mass;
+            double exhArea = engineCyl.ExhaustOrificeArea * 1e6;   // mm²
+            double intArea = engineCyl.IntakeOrificeArea * 1e6;    // mm²
+
+            // Exhaust pipe
+            int exhLast = exhaustPipe.GetCellCount() - 1;
+            double pExhEnd = exhaustPipe.GetCellPressure(exhLast);
+            double mdotExhOut = exhaustPipe.GetOpenEndMassFlow();   // positive out
+
+            // Intake pipe
+            double mdotIntIn = couplingIntake.LastMassFlowIntoVolume;
+            double pIntAmbEnd = intakePort.Pressure; 
+
+            Console.WriteLine(
+                $"{stepCount,8} {stroke,-11} {cycleDeg,6:F1}°  " +
+                $"RPM:{rpm,5:F0}  " +
+                $"Cyl: p={pCyl/1e5,6:F3}bar T={TCyl,6:F0}K V={VCyl*1e6,6:F0}cm³ m={mCyl*1e3,6:F6}g  " +
+                $"Valves: Exh={exhArea,5:F0}mm² Int={intArea,5:F0}mm²  " +
+                $"Intake: p_end={pIntAmbEnd/1e5,6:F3}bar mdot_in={mdotIntIn,7:F4}kg/s  " +
+                $"Exhaust: p_end={pExhEnd/1e5,6:F3}bar mdot_out={mdotExhOut,7:F4}kg/s");
+        }
+
+        private void UpdateExhaustPort()
+        {
+            int last = exhaustPipe.GetCellCount() - 1;
+            double p   = exhaustPipe.GetCellPressure(last);
+            double rho = exhaustPipe.GetCellDensity(last);
+            double vel = exhaustPipe.GetCellVelocity(last);
+
+            // Safety clamps
+            rho = Math.Clamp(rho, 0.01, 50.0);
+            vel = Math.Clamp(vel, -500.0, 500.0);
+            p   = Math.Clamp(p, 1e4, 2e6);
+
+            double outflowMassFlow = rho * vel * exhPipeArea;
+
+            exhaustPort.Pressure       = p;
+            exhaustPort.Density        = rho;
+            exhaustPort.Temperature    = p / (rho * 287.05);
+            exhaustPort.MassFlowRate   = -outflowMassFlow;
+            exhaustPort.SpecificEnthalpy = 0.0;
+        }
+
+        private void UpdateIntakePort()
+        {
+            // Use the actual valve mass flow (positive = into cylinder)
+            double mdotIntoEngine = couplingIntake.LastMassFlowIntoVolume;
+
+            // Use cylinder pressure/density for the port state (or intake pipe last cell)
+            double pCyl = engineCyl.Cylinder.Pressure;
+            double rhoCyl = engineCyl.Cylinder.Density;
+
+            intakePort.Pressure       = Math.Max(pCyl, 100);
+            intakePort.Density        = Math.Max(rhoCyl, 1e-6);
+            intakePort.Temperature    = engineCyl.Cylinder.Temperature;
+            intakePort.MassFlowRate   = mdotIntoEngine;
+            intakePort.SpecificEnthalpy = 0.0;
+        }
+
+        // ==================== Drawing ====================
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+            float centerY = winH / 2f;
+
+            const float T_ambient = 293.15f;
+            const float T_hot = 1500f;
+            const float T_cold = 0f;
+            const float R = 287.05f;
+            float deltaHot = T_hot - T_ambient;
+            float deltaCold = T_ambient - T_cold;
+
+            float NormaliseTemperature(double T)
+            {
+                double t;
+                if (T >= T_ambient)
+                    t = (T - T_ambient) / deltaHot;
+                else
+                    t = (T - T_ambient) / deltaCold;
+                return (float)Math.Clamp(t, -1.0, 1.0);
+            }
+
+            // ---- Cylinder ----
+            float cylW = 80f, cylH = 150f;
+            var cylRect = new RectangleShape(new Vector2f(cylW, cylH));
+            cylRect.Position = new Vector2f(200f, centerY - cylH / 2f);
+            double tempCyl = engineCyl.Cylinder.Temperature;
+            float tnCyl = NormaliseTemperature(tempCyl);
+            byte rC = (byte)(tnCyl > 0 ? 255 * tnCyl : 0);
+            byte bC = (byte)(tnCyl < 0 ? -255 * tnCyl : 0);
+            byte gC = (byte)(255 * (1 - Math.Abs(tnCyl)));
+            cylRect.FillColor = new Color(rC, gC, bC);
+            target.Draw(cylRect);
+
+            // ---- Piston ----
+            float pistonWidth = cylW - 12f;
+            float pistonHeight = 16f;
+            float pistonFraction = (float)engineCyl.PistonPositionFraction;
+            float pistonTopY = cylRect.Position.Y + pistonFraction * (cylH - pistonHeight);
+            var pistonRect = new RectangleShape(new Vector2f(pistonWidth, pistonHeight))
+            {
+                Position = new Vector2f(cylRect.Position.X + 6f, pistonTopY),
+                FillColor = new Color(80, 80, 80)
+            };
+            target.Draw(pistonRect);
+
+            // ---------- NEW: Valve lift indicators ----------
+            float barWidth = 30f;
+            float barHeight = 10f;
+            float exhLift = (float)engineCyl.ExhaustValveLiftCurrent;
+            float intLift = (float)engineCyl.IntakeValveLiftCurrent;
+
+            // Exhaust valve indicator (right side of cylinder)
+            var exhBar = new RectangleShape(new Vector2f(barWidth, barHeight))
+            {
+                Position = new Vector2f(cylRect.Position.X + cylW - 10,
+                                         cylRect.Position.Y - 20 - exhLift * 20),
+                FillColor = new Color(200, 200, 200)
+            };
+            target.Draw(exhBar);
+
+            // Intake valve indicator (left side of cylinder)
+            var intBar = new RectangleShape(new Vector2f(barWidth, barHeight))
+            {
+                Position = new Vector2f(cylRect.Position.X - 20,
+                                         cylRect.Position.Y - 20 - intLift * 20),
+                FillColor = new Color(200, 200, 200)
+            };
+            target.Draw(intBar);
+
+            // ---- Exhaust pipe (rightwards) ----
+            DrawPipe(target, exhaustPipe, startX: 280f, endX: winW - 60f, centerY + 10 - cylRect.Size.Y / 2,
+                     T_ambient, T_hot, T_cold, R, NormaliseTemperature, true);
+
+            // ---- Intake pipe (leftwards) ----
+            DrawPipe(target, intakePipe, startX: 200f, endX: 20f, centerY + 10 - cylRect.Size.Y / 2,
+                     T_ambient, T_hot, T_cold, R, NormaliseTemperature, false);
+        }
+
+        private void DrawPipe(RenderWindow target, Pipe1D pipe,
+            float startX, float endX, float centerY,
+            float T_ambient, float T_hot, float T_cold, float R,
+            Func<double, float> normaliseTemp, bool leftToRight)
+        {
+            int n = pipe.GetCellCount();
+            float dir = leftToRight ? 1f : -1f;
+            float pipeLen = Math.Abs(endX - startX);
+            float dx = pipeLen / (n - 1) * dir;
+            float baseRadius = leftToRight ? 20f : 14f;   // exhaust thicker, intake thinner
+            var vertices = new Vertex[n * 2];
+            float ambPress = 101325f;
+
+            for (int i = 0; i < n; i++)
+            {
+                float x = startX + i * dx;
+                double p = pipe.GetCellPressure(i);
+                double rho = pipe.GetCellDensity(i);
+                double T = p / (rho * R);
+
+                float r = baseRadius * 0.2f * (float)(1.0 + (p - ambPress) / ambPress);
+                if (r < 2f) r = 2f;
+
+                float tn = normaliseTemp(T);
+                byte rC = (byte)(tn > 0 ? 255 * tn : 0);
+                byte bC = (byte)(tn < 0 ? -255 * tn : 0);
+                byte gC = (byte)(255 * (1 - Math.Abs(tn)));
+                var col = new Color(rC, gC, bC);
+
+                vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
+                vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
+            }
+
+            target.Draw(vertices, PrimitiveType.TriangleStrip);
+        }
+    }
+}

Scenarios/HelmholtzResonatorScenario.cs

@@ -1,6 +1,5 @@
 using System;
 using FluidSim.Components;
-using FluidSim.Interfaces;
 using FluidSim.Utils;
 using SFML.Graphics;
 using SFML.System;
@@ -12,7 +11,7 @@ namespace FluidSim.Core
         private Solver solver;
         private Volume0D cavity;
         private Pipe1D neck;
-        private Connection coupling;
+        private PipeVolumeConnection coupling;
         private int stepCount;
         private double time;
         private double dt;
@@ -38,12 +37,8 @@ namespace FluidSim.Core
             neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40);
             neck.SetUniformState(1.225, 0.0, ambientPressure);
 
-            coupling = new Connection(neck.PortA, cavity.Port)
-            {
-                Area = neckArea,
-                DischargeCoefficient = 0.62,
-                Gamma = 1.4
-            };
+            // Create the coupling between cavity and left end of the neck (PortA)
+            coupling = new PipeVolumeConnection(cavity, neck, isPipeLeftEnd: true, orificeArea: neckArea);
 
             solver = new Solver();
             solver.SetTimeStep(dt);
@@ -51,8 +46,8 @@ namespace FluidSim.Core
             solver.AddPipe(neck);
             solver.AddConnection(coupling);
 
-            // Port A (left) = volume coupling, Port B (right) = open end
-            solver.SetPipeBoundary(neck, isA: true, BoundaryType.VolumeCoupling);
+            // Left boundary (PortA) is volume‑coupled via ghost cell, right boundary (PortB) is open end
+            solver.SetPipeBoundary(neck, isA: true, BoundaryType.GhostCell);
             solver.SetPipeBoundary(neck, isA: false, BoundaryType.OpenEnd, ambientPressure);
         }
 
@@ -68,11 +63,11 @@ namespace FluidSim.Core
             if (stepCount % 20 == 0)
             {
                 double pCav = cavity.Pressure;
-                double mdotA = neck.PortA.MassFlowRate;   // positive = into pipe (leaving cavity)
+                // Mass flow rate is not directly available – we can compute from pressure difference or skip
                 Console.WriteLine(
                     $"t={time * 1e3:F2} ms step={stepCount} " +
                     $"P_cav={pCav:F1} Pa, P_open={pOpen:F1} Pa, " +
-                    $"mdot_A={mdotA * 1e3:F4} g/s, audio={audio:F4}");
+                    $"audio={audio:F4}");
             }
 
             return audio;
@@ -100,7 +95,7 @@ namespace FluidSim.Core
             float dx = neckLenPx / (n - 1);
             float baseRadius = 20f;
 
-            Vertex[] vertices = new Vertex[n * 2];
+            var vertices = new Vertex[n * 2];
             for (int i = 0; i < n; i++)
             {
                 float x = neckStartX + i * dx;