config tuning

.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"
+        }
+    ]
+}

Audio/AudioOutputStream.cs

@@ -3,22 +3,31 @@ using SFML.System;
 
 namespace FluidSim.Audio
 {
-    internal class RingBufferStream : SoundStream
+    public class AudioOutputStream : SoundStream
     {
-        private readonly RingBuffer ringBuffer;
+        private readonly SimulationRingBuffer _sourceBuffer;
+        private double _speed = 1.0;   // non‑volatile, accessed with Volatile.Read/Write
 
-        public RingBufferStream(RingBuffer buffer)
+        public AudioOutputStream(SimulationRingBuffer sourceBuffer)
         {
-            ringBuffer = buffer;
+            _sourceBuffer = sourceBuffer;
-            // 2 channels, 44.1 kHz, standard stereo mapping
+            // 2 channels, 44.1 kHz, stereo
             Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
         }
 
+        /// <summary>Playback speed (0.01 … 1.0 or higher for catch‑up).</summary>
+        public double Speed
+        {
+            get => Volatile.Read(ref _speed);
+            set => Volatile.Write(ref _speed, value);
+        }
+
         protected override bool OnGetData(out short[] samples)
         {
-            const int monoBlockSize = 512;          // number of mono samples we'll read
+            const int monoBlockSize = 512;
             float[] temp = new float[monoBlockSize];
-            int read = ringBuffer.Read(temp, monoBlockSize);
+
+            int read = _sourceBuffer.ReadInterpolated(temp, monoBlockSize, Speed);
             samples = new short[monoBlockSize * 2];
 
             if (read > 0)
@@ -31,6 +40,7 @@ namespace FluidSim.Audio
                     samples[i * 2 + 1] = final;  // right
                 }
             }
+            // Fill rest with silence
             for (int i = read * 2; i < samples.Length; i++)
                 samples[i] = 0;
 

Audio/RingBuffer.cs

@@ -1,41 +0,0 @@
-namespace FluidSim.Audio
-{
-    internal class RingBuffer
-    {
-        private readonly float[] buffer;
-        private volatile int readPos;
-        private volatile int writePos;
-
-        public RingBuffer(int capacity)
-        {
-            if ((capacity & (capacity - 1)) != 0)
-                throw new ArgumentException("Capacity must be a power of two.");
-            buffer = new float[capacity];
-        }
-
-        public int Count => (writePos - readPos) & (buffer.Length - 1);
-        public int Space => (readPos - writePos - 1) & (buffer.Length - 1);
-
-        public int Write(float[] data, int count)
-        {
-            int space = Space;
-            int toWrite = Math.Min(count, space);
-            int mask = buffer.Length - 1;
-            for (int i = 0; i < toWrite; i++)
-                buffer[(writePos + i) & mask] = data[i];
-            writePos = (writePos + toWrite) & mask;
-            return toWrite;
-        }
-
-        public int Read(float[] destination, int count)
-        {
-            int available = Count;
-            int toRead = Math.Min(count, available);
-            int mask = buffer.Length - 1;
-            for (int i = 0; i < toRead; i++)
-                destination[i] = buffer[(readPos + i) & mask];
-            readPos = (readPos + toRead) & mask;
-            return toRead;
-        }
-    }
-}

Audio/SimulationRingBuffer.cs

@@ -0,0 +1,98 @@
+namespace FluidSim.Audio
+{
+    public class SimulationRingBuffer
+    {
+        private readonly float[] _buffer;
+        private readonly int _capacity;
+        private int _writeHead;      // monotonic, producer only
+        private int _readHead;       // monotonic, consumer advances after consumption
+
+        // Consumer interpolation state
+        private double _readPosFrac;
+        private bool _consumerInit;
+
+        // Events for signalling
+        private readonly AutoResetEvent _spaceAvailable = new AutoResetEvent(false);
+        private readonly AutoResetEvent _dataAvailable = new AutoResetEvent(false);
+
+        public SimulationRingBuffer(int capacity)
+        {
+            if ((capacity & (capacity - 1)) != 0)
+                throw new ArgumentException("Capacity must be a power of two.");
+            _capacity = capacity;
+            _buffer = new float[capacity];
+        }
+
+        // ---------- Producer ----------
+        public int FreeSpace => _capacity - (_writeHead - Volatile.Read(ref _readHead));
+
+        /// <summary>Number of samples currently available for reading (integer count).</summary>
+        public int AvailableSamples => Volatile.Read(ref _writeHead) - Volatile.Read(ref _readHead);
+
+        public void Write(float sample)
+        {
+            while (FreeSpace == 0)
+                _spaceAvailable.WaitOne();
+
+            int w = _writeHead;
+            int mask = _capacity - 1;
+            _buffer[w & mask] = sample;
+            Volatile.Write(ref _writeHead, w + 1);
+            _dataAvailable.Set();
+        }
+
+        public int Write(float[] data, int count)
+        {
+            int free = FreeSpace;
+            int toWrite = Math.Min(count, free);
+            int w = _writeHead;
+            int mask = _capacity - 1;
+            for (int i = 0; i < toWrite; i++)
+                _buffer[(w + i) & mask] = data[i];
+            Volatile.Write(ref _writeHead, w + toWrite);
+            if (toWrite > 0)
+                _dataAvailable.Set();
+            return toWrite;
+        }
+
+        // ---------- Consumer ----------
+        public void ResetConsumer() => _consumerInit = false;
+
+        public int ReadInterpolated(float[] dest, int destCount, double speed)
+        {
+            if (!_consumerInit)
+            {
+                _readPosFrac = Volatile.Read(ref _readHead);
+                _consumerInit = true;
+            }
+
+            int mask = _capacity - 1;
+            int writeHead = Volatile.Read(ref _writeHead);
+            int produced = 0;
+
+            for (int i = 0; i < destCount; i++)
+            {
+                int idxFloor = (int)_readPosFrac;
+                int idxCeil = idxFloor + 1;
+                if (idxCeil >= writeHead)
+                    break;
+
+                float y0 = _buffer[idxFloor & mask];
+                float y1 = _buffer[idxCeil & mask];
+                double frac = _readPosFrac - idxFloor;
+                dest[i] = (float)(y0 + (y1 - y0) * frac);
+
+                _readPosFrac += speed;
+                produced++;
+            }
+
+            int newReadHead = (int)_readPosFrac;
+            if (newReadHead > Volatile.Read(ref _readHead))
+            {
+                Volatile.Write(ref _readHead, newReadHead);
+                _spaceAvailable.Set();
+            }
+            return produced;
+        }
+    }
+}

Audio/SoundEngine.cs

@@ -1,45 +1,45 @@
-namespace FluidSim.Audio;
+namespace FluidSim.Audio
-
-public class SoundEngine : IDisposable
 {
-    private readonly RingBuffer ringBuffer;
+    public class SoundEngine : IDisposable
-    private readonly RingBufferStream stream;
-    private bool isPlaying;
-
-    public SoundEngine(int bufferCapacity = 16384)
     {
-        ringBuffer = new RingBuffer(bufferCapacity);
+        private readonly AudioOutputStream _stream;
-        stream = new RingBufferStream(ringBuffer);
+        private bool _isPlaying;
+
+        public SoundEngine(SimulationRingBuffer sourceBuffer, int bufferCapacity = 16384)
+        {
+            _stream = new AudioOutputStream(sourceBuffer);
         }
 
         public void Start()
         {
-        if (isPlaying) return;
+            if (_isPlaying) return;
-        stream.Play();
+            _stream.Play();
-        isPlaying = true;
+            _isPlaying = true;
         }
 
         public void Stop()
         {
-        if (!isPlaying) return;
+            if (!_isPlaying) return;
-        stream.Stop();
+            _stream.Stop();
-        isPlaying = false;
+            _isPlaying = false;
-        float[] drain = new float[ringBuffer.Count];
-        ringBuffer.Read(drain, drain.Length);
         }
 
-    public int WriteSamples(float[] data, int count) =>
+        public double Speed
-        ringBuffer.Write(data, count);
+        {
+            get => _stream.Speed;
+            set => _stream.Speed = value;
+        }
 
         public float Volume
         {
-        get => stream.Volume;
+            get => _stream.Volume;
-        set => stream.Volume = value;
+            set => _stream.Volume = value;
         }
 
         public void Dispose()
         {
             Stop();
-        stream.Dispose();
+            _stream.Dispose();
+        }
     }
 }

Audio/SoundProcessor.cs

@@ -1,29 +0,0 @@
-using FluidSim.Components;
-using FluidSim.Interfaces;
-using System;
-
-namespace FluidSim.Audio
-{
-    public static class SoundProcessor
-    {
-        public static float MaxDeltaP { get; set; } = 100_000f;
-        public static float MaxArea { get; set; } = 1e-4f;
-        public static float MaxVelocity { get; set; } = 343f;
-        public static float ReferenceDensity { get; set; } = 1.225f;
-        public static float ReferenceSpeedOfSound { get; set; } = 343f;
-        public static float Gain { get; set; } = 1.0f;
-
-        public static double ComputeSample(Connection conn)
-        {
-            Port portA = conn.PortA;
-            Port portB = conn.PortB;
-
-            double pressureUp = portA.Pressure;
-            double pressureDown = portB.Pressure;
-
-            // No flow or no pressure drop → silence
-            double deltaP = pressureUp - pressureDown;
-            return deltaP / 1;
-        }
-    }
-}

Components/Atmosphere.cs

@@ -0,0 +1,39 @@
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    public class Atmosphere : IComponent
+    {
+        public float Pressure { get; set; } = 101325f;
+        public float Temperature { get; set; } = 300f;
+        public float GasConstant { get; set; } = 287f;
+        public float Gamma => 1.4f;
+
+        public float Density => Pressure / (GasConstant * Temperature);
+        public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / Density;
+
+        public Port Port { get; }
+
+        public Atmosphere()
+        {
+            Port = new Port { Owner = this };
+            UpdatePort();
+        }
+
+        public IReadOnlyList<Port> Ports => new[] { Port };
+
+        public void UpdateState(float dt)
+        {
+            UpdatePort();
+        }
+
+        private void UpdatePort()
+        {
+            Port.Pressure = Pressure;
+            Port.Density = Density;
+            Port.Temperature = Temperature;
+            Port.SpecificEnthalpy = SpecificEnthalpy;
+            Port.AirFraction = 1f;
+        }
+    }
+}

Components/Connection.cs

@@ -1,17 +0,0 @@
-using FluidSim.Interfaces;
-
-namespace FluidSim.Components
-{
-    /// <summary>Pure data link between two ports, with orifice parameters.</summary>
-    public class Connection
-    {
-        public Port PortA { get; }
-        public Port PortB { get; }
-
-        public double Area { get; set; } = 1e-5;          // effective orifice area (m²)
-        public double DischargeCoefficient { get; set; } = 0.62;
-        public double Gamma { get; set; } = 1.4;
-
-        public Connection(Port a, Port b) => (PortA, PortB) = (a, b);
-    }
-}

Components/Crankshaft.cs

@@ -0,0 +1,52 @@
+using System;
+
+namespace FluidSim.Components
+{
+    public class Crankshaft
+    {
+        public float AngularVelocity;   // rad/s
+        public float CrankAngle;        // rad, 0 … 4π
+        public float PreviousAngle;
+
+        public float Inertia = 0.2f;
+        public float FrictionConstant;   // N·m
+        public float FrictionViscous;    // N·m per rad/s
+
+        private float externalTorque;
+
+        public Crankshaft(float initialRPM = 400f)
+        {
+            AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
+            CrankAngle = 0f;
+            PreviousAngle = 0f;
+        }
+
+        public void AddTorque(float torque) => externalTorque += torque;
+
+        public void Step(float dt)
+        {
+            if (float.IsNaN(AngularVelocity) || float.IsInfinity(AngularVelocity))
+                AngularVelocity = 0f;
+            if (float.IsNaN(externalTorque) || float.IsInfinity(externalTorque))
+                externalTorque = 0f;
+
+            PreviousAngle = CrankAngle;
+
+            float friction = FrictionConstant * MathF.Sign(AngularVelocity)
+                           + FrictionViscous * AngularVelocity;
+            float netTorque = externalTorque - friction;
+            float alpha = netTorque / Inertia;
+            AngularVelocity += alpha * dt;
+
+            if (AngularVelocity < 0f) AngularVelocity = 0f;
+
+            CrankAngle += AngularVelocity * dt;
+            if (CrankAngle >= 4f * MathF.PI)
+                CrankAngle -= 4f * MathF.PI;
+            else if (CrankAngle < 0f)
+                CrankAngle += 4f * MathF.PI;
+
+            externalTorque = 0f;
+        }
+    }
+}

Components/Cylinder.cs

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

Components/Pipe1D.cs

@@ -1,250 +0,0 @@
-using System;
-using FluidSim.Interfaces;
-
-namespace FluidSim.Components
-{
-    public class Pipe1D
-    {
-        public Port PortA { get; }
-        public Port PortB { get; }
-        public double Area => _area;
-
-        private int _n;
-        private double _dx, _dt, _gamma = 1.4, _area;
-        private double[] _rho, _rhou, _E;
-        private double _hydraulicDiameter;
-
-        private double _rhoLeft, _pLeft, _rhoRight, _pRight;
-        private bool _leftBCSet, _rightBCSet;
-
-        public double FrictionFactor { get; set; }
-
-        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);
-
-        /// <summary>
-        /// Create a pipe with CFL‑stable automatic cell count.
-        /// </summary>
-        /// <param name="length">Pipe length [m].</param>
-        /// <param name="area">Cross‑sectional area [m²].</param>
-        /// <param name="sampleRate">Simulation step rate [Hz].</param>
-        /// <param name="c0">Speed of sound [m/s] (default 343).</param>
-        /// <param name="frictionFactor">Darcy friction factor (0 = inviscid).</param>
-        /// <param name="cflSafety">CFL safety factor ≤ 1 (0.8 recommended).</param>
-        public Pipe1D(double length, double area, int sampleRate,
-                      double c0 = 343.0, double frictionFactor = 0.02,
-                      double cflSafety = 0.8)
-        {
-            if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
-            _area = area;
-            _dt = 1.0 / sampleRate;
-            FrictionFactor = frictionFactor;
-
-            // Nyquist‑based cell count (wave resolution)
-            double nNyquist = Math.Ceiling(length * sampleRate / c0);
-            // CFL‑stable cell count: dx ≥ maxSpeed·dt / cflSafety, maxSpeed = 2·c0 (supersonic safe)
-            double maxSpeed = 2.0 * c0;
-            double dxMinStable = maxSpeed * _dt / cflSafety;
-            double nStable = Math.Floor(length / dxMinStable);
-
-            _n = Math.Max(2, (int)Math.Min(nNyquist, nStable));
-            _dx = length / _n;
-
-            _rho = new double[_n];
-            _rhou = new double[_n];
-            _E = new double[_n];
-
-            _hydraulicDiameter = Math.Max(2.0 * Math.Sqrt(_area / Math.PI), 1e-9);
-
-            PortA = new Port();
-            PortB = new Port();
-        }
-
-        public void SetUniformState(double rho, double u, double p)
-        {
-            double e = p / ((_gamma - 1) * rho);
-            double Etot = rho * e + 0.5 * rho * u * u;
-            for (int i = 0; i < _n; i++)
-            {
-                _rho[i] = rho;
-                _rhou[i] = rho * u;
-                _E[i] = Etot;
-            }
-        }
-
-        public double GetLeftPressure() => Pressure(0);
-        public double GetRightPressure() => Pressure(_n - 1);
-        public double GetLeftDensity() => _rho[0];
-        public double GetRightDensity() => _rho[_n - 1];
-
-        public void SetLeftVolumeState(double rhoVol, double pVol)
-        {
-            _rhoLeft = rhoVol;
-            _pLeft = pVol;
-            _leftBCSet = true;
-        }
-
-        public void SetRightVolumeState(double rhoVol, double pVol)
-        {
-            _rhoRight = rhoVol;
-            _pRight = pVol;
-            _rightBCSet = true;
-        }
-
-        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;
-        }
-
-        public void Simulate()
-        {
-            int n = _n;
-            double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];
-
-            // --- Left boundary (face 0) ---
-            if (_leftBCSet)
-            {
-                double rhoL = _rhoLeft, uL = 0.0, pL = _pLeft;
-                double rhoR = _rho[0], uR = _rhou[0] / Math.Max(rhoR, 1e-12), pR = Pressure(0);
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]);
-            }
-            else
-            {
-                Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0;
-            }
-
-            // --- Internal faces ---
-            for (int i = 0; i < n - 1; i++)
-            {
-                double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
-                double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
-                HLLCFlux(_rho[i], uL, Pressure(i), _rho[i + 1], uR, Pressure(i + 1),
-                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
-            }
-
-            // --- Right boundary (face n) ---
-            if (_rightBCSet)
-            {
-                double rhoL = _rho[n - 1], uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12), pL = Pressure(n - 1);
-                double rhoR = _rhoRight, uR = 0.0, pR = _pRight;
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]);
-            }
-            else
-            {
-                Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0;
-            }
-
-            // --- Cell update (inviscid fluxes) ---
-            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] -= _dt * dM;
-                _rhou[i] -= _dt * dP;
-                _E[i] -= _dt * dE;
-
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
-                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
-                if (_E[i] < kinetic) _E[i] = kinetic;
-
-                // Emergency reset if NaN
-                if (double.IsNaN(_rho[i]) || double.IsNaN(_rhou[i]) || double.IsNaN(_E[i]))
-                {
-                    _rho[i] = 1.225;                   // reset to atmospheric air at 300 K
-                    _rhou[i] = 0.0;
-                    _E[i] = 101325.0 / (_gamma - 1.0); // internal energy at 1 atm
-                }
-            }
-
-            // --- Friction (Darcy–Weisbach, energy‑conserving) ---
-            if (FrictionFactor > 0)
-            {
-                double D = _hydraulicDiameter;
-                double twoD = 2.0 * D;
-                for (int i = 0; i < n; i++)
-                {
-                    double rho = _rho[i];
-                    double u = _rhou[i] / rho;
-                    double absU = Math.Abs(u);
-                    double src = FrictionFactor * rho * absU * u / twoD;
-                    double kinOld = 0.5 * rho * u * u;
-                    _rhou[i] -= _dt * src;
-                    double uNew = _rhou[i] / rho;
-                    double kinNew = 0.5 * rho * uNew * uNew;
-                    _E[i] += (kinOld - kinNew);
-                }
-            }
-
-            // --- Publish to ports ---
-            PortA.Pressure = Pressure(0);
-            PortA.Density = _rho[0];
-            PortB.Pressure = Pressure(_n - 1);
-            PortB.Density = _rho[_n - 1];
-
-            PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0;
-            PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0;
-
-            PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0);
-            PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1);
-
-            _leftBCSet = _rightBCSet = false;
-        }
-
-        double Pressure(int i) =>
-            (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
-
-        void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
-                      out double fm, out double fp, out double fe)
-        {
-            const double eps = 1e-12;
-            pL = Math.Max(pL, eps);
-            pR = Math.Max(pR, eps);
-
-            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, eps));
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, eps));
-            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 denom = rL * (SL - uL) - rR * (SR - uR);
-            double Ss;
-            if (Math.Abs(denom) < eps)
-                Ss = 0.5 * (uL + uR);
-            else
-                Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
-
-            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 diffSL = SL - uL;
-                if (Math.Abs(diffSL) < eps) diffSL = eps;
-                double rsL = rL * diffSL / (SL - Ss);
-                double ps = pL + rL * diffSL * (Ss - uL);
-                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * diffSL));
-                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
-            }
-            else
-            {
-                double diffSR = SR - uR;
-                if (Math.Abs(diffSR) < eps) diffSR = eps;
-                double rsR = rR * diffSR / (SR - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
-                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * diffSR));
-                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
-            }
-        }
-    }
-}

Components/SoundConnection.cs

@@ -1,9 +0,0 @@
-using FluidSim.Interfaces;
-
-namespace FluidSim.Components
-{
-    public class SoundConnection : Connection
-    {
-        public SoundConnection(Port a, Port b) : base(a, b) { }
-    }
-}

Components/Volume0D.cs

@@ -1,69 +1,133 @@
-using FluidSim.Interfaces;
+using System;
+using System.Collections.Generic;
+using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
-    public class Volume0D
+    public class Volume0D : IComponent
     {
-        public Port Port { get; private set; }
+        public List<Port> Ports { get; } = new List<Port>();
 
-        public double Mass { get; private set; }
+        private float _airMass;
-        public double InternalEnergy { get; private set; }
+        private float _exhaustMass;
+        public float InternalEnergy;
+        public float Volume;
+        public float Dvdt;
+        public float Gamma { get; set; } = 1.4f;
+        public float GasConstant { get; set; } = 287f;
+        public float AmbientPressure { get; set; } = 101325f;
 
-        public double Gamma { get; set; } = 1.4;
+        // ---------- Thermal relaxation to environment ----------
-        public double GasConstant { get; set; } = 287.0;
+        /// <summary>Rate of heat transfer to the surroundings (1/s). 0 = adiabatic.</summary>
+        public float EnergyRelaxationRate { get; set; } = 10f;
+        /// <summary>Temperature to relax toward (K). Default is room temperature.</summary>
+        public float AmbientTemperature { get; set; } = 300f;
 
-        public double Volume { get; set; }
+        public float Mass => _airMass + _exhaustMass;
-        public double dVdt { get; set; }
+        public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
+        public float Density => Mass / MathF.Max(Volume, 1e-12f);
+        public float Pressure => (Gamma - 1f) * InternalEnergy / MathF.Max(Volume, 1e-12f);
+        public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
+        public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / MathF.Max(Density, 1e-12f);
 
-        private double _dt;
+        public Volume0D(float initialVolume, float initialPressure,
-
+                        float initialTemperature, float gasConstant = 287f, float gamma = 1.4f)
-        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 Volume0D(double initialVolume, double initialPressure,
-                        double initialTemperature, int sampleRate,
-                        double gasConstant = 287.0, double gamma = 1.4)
         {
             GasConstant = gasConstant;
             Gamma = gamma;
             Volume = initialVolume;
-            dVdt = 0.0;
+            Dvdt = 0f;
-            _dt = 1.0 / sampleRate;
 
-            double rho0 = initialPressure / (GasConstant * initialTemperature);
+            float rho0 = initialPressure / (GasConstant * initialTemperature);
-            Mass = rho0 * Volume;
+            _airMass = rho0 * Volume;
-            InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
+            _exhaustMass = 0f;
-
+            InternalEnergy = (initialPressure * Volume) / (Gamma - 1f);
-            Port = new Port();
-            PushStateToPort();
         }
 
-        public void PushStateToPort()
+        public Port CreatePort()
         {
-            Port.Pressure = Pressure;
+            var port = new Port { Owner = this };
-            Port.Density = Density;
+            port.Pressure = Pressure;
-            Port.Temperature = Temperature;
+            port.Density = Density;
-            Port.SpecificEnthalpy = SpecificEnthalpy;
+            port.Temperature = Temperature;
+            port.SpecificEnthalpy = SpecificEnthalpy;
+            port.AirFraction = AirFraction;
+            Ports.Add(port);
+            return port;
         }
 
-        public void Integrate()
+        public void SetPressure(float pressure, float? temperature = null)
         {
-            double mdot = Port.MassFlowRate;
+            float V = MathF.Max(Volume, 1e-12f);
-            double h_in = Port.SpecificEnthalpy;
+            float T = temperature ?? Temperature;
+            float rho = pressure / (GasConstant * T);
+            float totalMass = rho * V;
+            float af = AirFraction;
+            _airMass = totalMass * af;
+            _exhaustMass = totalMass * (1f - af);
+            InternalEnergy = pressure * V / (Gamma - 1f);
+        }
 
-            double dm = mdot * _dt;
+        public void UpdateState(float dt)
-            double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
+        {
+            float totalMdotAir = 0f, totalMdotExhaust = 0f, totalEdot = 0f;
+            foreach (var port in Ports)
+            {
+                float mdot = port.MassFlowRate;
+                float af = mdot >= 0f ? port.AirFraction : AirFraction;
+                totalMdotAir += mdot * af;
+                totalMdotExhaust += mdot * (1f - af);
+                totalEdot += mdot * port.SpecificEnthalpy;
+            }
 
-            Mass += dm;
+            float dAir = totalMdotAir * dt;
+            float dExhaust = totalMdotExhaust * dt;
+            float dE = totalEdot * dt - Pressure * Dvdt * dt;
+
+            _airMass += dAir;
+            _exhaustMass += dExhaust;
             InternalEnergy += dE;
 
-            // Hard physical bounds – prevent NaN and unphysical states
+            // ---- Thermal relaxation ----
-            if (Mass < 1e-12) Mass = 1e-12;
+            if (EnergyRelaxationRate > 0f)
-            if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
+            {
-
+                float currentMass = Mass;
-            PushStateToPort();
+                if (currentMass > 1e-12f)
+                {
+                    // Target internal energy: current mass at ambient temperature
+                    float targetE = currentMass * GasConstant * AmbientTemperature / (Gamma - 1f);
+                    float relaxFactor = MathF.Exp(-EnergyRelaxationRate * dt);
+                    InternalEnergy = targetE + (InternalEnergy - targetE) * relaxFactor;
                 }
             }
+
+            float V = MathF.Max(Volume, 1e-12f);
+            float totalMass = _airMass + _exhaustMass;
+            if (totalMass < 1e-9f)
+            {
+                _airMass = 1e-9f;
+                _exhaustMass = 0f;
+                InternalEnergy = AmbientPressure * V / (Gamma - 1f);
+            }
+            else if (InternalEnergy < 0f)
+            {
+                InternalEnergy = AmbientPressure * V / (Gamma - 1f);
+            }
+
+            if (_airMass < 0f) _airMass = 0f;
+            if (_exhaustMass < 0f) _exhaustMass = 0f;
+
+            float p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afr = AirFraction;
+            foreach (var port in Ports)
+            {
+                port.Pressure = p;
+                port.Density = rho;
+                port.Temperature = T;
+                port.SpecificEnthalpy = h;
+                port.AirFraction = afr;
+            }
+        }
+
+        IReadOnlyList<Port> IComponent.Ports => Ports;
+    }
 }

Core/BoundarySystem.cs

@@ -0,0 +1,404 @@
+using FluidSim.Components;
+using FluidSim.Interfaces;
+using System;
+
+namespace FluidSim.Core
+{
+    public class BoundarySystem
+    {
+        // ---------- Private constants ----------
+        private const float Gamma = 1.4f;
+        private const float Gm1 = Gamma - 1f;   // 0.4
+        private const float Rgas = 287f;         // J/(kg·K)
+        private const float GammaOverGm1 = Gamma / Gm1; // 3.5
+
+        public struct OrificeDesc
+        {
+            public Port VolumePort;
+            public int PipeIndex;
+            public bool IsLeftEnd;
+            public int AreaIndex;
+            public float DischargeCoeff;
+
+            // --- Inertance support ---
+            public bool UseInertance;
+            public float EffectiveLength;
+            public float CurrentMdot;       // kg/s, positive = volume → pipe
+
+            // --- Loss coefficient (linear resistance) ---
+            public float LossCoefficient;   // N·s/m⁵  or  kg/(m⁴·s)
+        }
+
+        public struct OpenEndDesc
+        {
+            public int PipeIndex;
+            public bool IsLeftEnd;
+            public float AmbientPressure;
+            public float Gamma;
+            public float PipeArea;
+            public float LastMassFlowRate;
+            public float LastFacePressure;
+        }
+
+        private OrificeDesc[] _orifices;
+        private OpenEndDesc[] _openEnds;
+        private float[] _orificeAreas;
+        private PipeSystem _pipeSystem;
+
+        public BoundarySystem(PipeSystem pipeSystem, int maxOrifices, int maxOpenEnds)
+        {
+            _pipeSystem = pipeSystem;
+            _orifices = new OrificeDesc[maxOrifices];
+            _openEnds = new OpenEndDesc[maxOpenEnds];
+            _orificeAreas = new float[maxOrifices];
+        }
+
+        public int OrificeCount { get; private set; }
+        public int OpenEndCount { get; private set; }
+
+        // ---------- Add orifice (no inertance) ----------
+        public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
+                               int areaIndex, float dischargeCoeff = 1f,
+                               float lossCoefficient = 0f)
+        {
+            _orifices[OrificeCount] = new OrificeDesc
+            {
+                VolumePort = volumePort,
+                PipeIndex = pipeIndex,
+                IsLeftEnd = isLeftEnd,
+                AreaIndex = areaIndex,
+                DischargeCoeff = dischargeCoeff,
+                UseInertance = false,
+                EffectiveLength = 0f,
+                CurrentMdot = 0f,
+                LossCoefficient = lossCoefficient
+            };
+            OrificeCount++;
+        }
+
+        // ---------- Add orifice with inertance ----------
+        public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
+                                            int areaIndex, float dischargeCoeff,
+                                            float effectiveLength, float lossCoefficient = 0f)
+        {
+            // Reuse the base AddOrifice and then override fields
+            AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient);
+            ref var d = ref _orifices[OrificeCount - 1];
+            d.UseInertance = true;
+            d.EffectiveLength = effectiveLength;
+            d.LossCoefficient = lossCoefficient;   // store the linear resistance
+        }
+
+        public void AddOpenEnd(int pipeIndex, bool isLeftEnd,
+                               float ambientPressure, float pipeArea, float gamma = 1.4f)
+        {
+            int idx = OpenEndCount;
+            _openEnds[idx] = new OpenEndDesc
+            {
+                PipeIndex = pipeIndex,
+                IsLeftEnd = isLeftEnd,
+                AmbientPressure = ambientPressure,
+                Gamma = gamma,
+                PipeArea = pipeArea
+            };
+            OpenEndCount++;
+        }
+
+        public void SetOrificeAreas(float[] areas)
+        {
+            for (int i = 0; i < OrificeCount; i++)
+                _orificeAreas[i] = areas[i];
+        }
+
+        public float GetOpenEndMassFlow(int openEndIndex)
+        {
+            if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 0f;
+            return _openEnds[openEndIndex].LastMassFlowRate;
+        }
+
+        public float GetOpenEndPressure(int openEndIndex)
+        {
+            if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 101325f;
+            return _openEnds[openEndIndex].LastFacePressure;
+        }
+
+        // ---------- Resolve all orifices ----------
+        public void ResolveOrifices(float dt)
+        {
+            for (int i = 0; i < OrificeCount; i++)
+            {
+                ref var d = ref _orifices[i];
+                float area = _orificeAreas[d.AreaIndex];
+
+                // Gather volume state
+                float volP = d.VolumePort?.Pressure ?? 101325f;
+                float volRho = d.VolumePort?.Density ?? 1.2f;
+                float volT = d.VolumePort?.Temperature ?? 300f;
+                float volH = d.VolumePort?.SpecificEnthalpy ?? 0f;
+                float volAF = d.VolumePort?.AirFraction ?? 1f;
+
+                // Gather pipe interior state
+                var (pipeRho, pipeU, pipeP) = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
+                float pipeT = pipeP / MathF.Max(pipeRho * Rgas, 1e-12f);
+                float pipeAF = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
+
+                // ---- Handle closed orifice (area ≈ 0) as a wall ----
+                if (area < 1e-12f || d.VolumePort == null)
+                {
+                    var (rInt, uInt, pInt) = d.IsLeftEnd
+                        ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
+                        : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
+                    float afInt = d.IsLeftEnd
+                        ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
+                        : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
+
+                    if (d.IsLeftEnd)
+                        _pipeSystem.SetGhostLeft(d.PipeIndex, rInt, -uInt, pInt, afInt);
+                    else
+                        _pipeSystem.SetGhostRight(d.PipeIndex, rInt, -uInt, pInt, afInt);
+
+                    if (d.VolumePort != null) d.VolumePort.MassFlowRate = 0f;
+                    continue;
+                }
+
+                // ---- Preliminary isentropic solution ----
+                float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
+                if (volP >= pipeP)
+                {
+                    IsentropicOrifice.Compute(volP, volRho, volT, pipeP, Gamma, Rgas, area, d.DischargeCoeff,
+                        out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
+                }
+                else
+                {
+                    IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, Gamma, Rgas, area, d.DischargeCoeff,
+                        out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
+                    mdotEst = -mdotEst;
+                }
+
+                // ---- Compute final mass flow with limiters ----
+                float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
+
+                if (d.UseInertance)
+                {
+                    float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
+                    float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
+                    float dp = volP - pipeP;
+                    float Rlin = d.LossCoefficient;
+
+                    float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
+                    float mdotNew = d.CurrentMdot + dmdot_dt * dt;
+
+                    // Limit outflow from volume (if volume owner is Volume0D)
+                    if (d.VolumePort.Owner is Volume0D vol0)
+                    {
+                        float maxOut = vol0.Mass / dt;
+                        if (mdotNew > maxOut) mdotNew = maxOut;
+                        if (mdotNew < -maxOut) mdotNew = -maxOut;
+                    }
+
+                    // Limit inflow from pipe – pipe cell cannot be emptied in one step
+                    {
+                        int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
+                                                : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
+                        float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
+                        float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);  // true cell area, not orifice
+                        float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
+                        float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
+                        float maxFromPipe = pipeCellMass / dt;
+                        if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
+                    }
+
+                    // Velocity clamp to Mach 0.9
+                    float rhoFacePrelim = mdotNew >= 0 ? volRho : pipeRho;
+                    float uFacePrelim = MathF.Abs(mdotNew) / MathF.Max(rhoFacePrelim * area, 1e-12f);
+                    float cUp = mdotNew >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
+                    float maxU = 0.9f * cUp;
+                    if (uFacePrelim > maxU)
+                    {
+                        uFacePrelim = maxU;
+                        mdotNew = rhoFacePrelim * uFacePrelim * area * (mdotNew >= 0 ? 1f : -1f);
+                    }
+
+                    if (float.IsNaN(mdotNew)) mdotNew = 0f;
+
+                    d.CurrentMdot = mdotNew;
+                    mdotFinal = mdotNew;
+                    rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
+                    pFace = pFaceEst;
+                    uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
+                }
+                else
+                {
+                    // Standard quasi‑steady orifice
+                    mdotFinal = mdotEst;
+                    rhoFace = rhoFaceEst;
+                    uFace = uFaceEst;
+                    pFace = pFaceEst;
+
+                    // Limit outflow from volume (if Volume0D)
+                    if (d.VolumePort.Owner is Volume0D vol0)
+                    {
+                        float maxOut = vol0.Mass / dt;
+                        if (mdotFinal > maxOut) mdotFinal = maxOut;
+                    }
+
+                    // ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
+                    if (mdotFinal < 0)
+                    {
+                        int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
+                                                : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
+                        float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
+                        float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
+                        float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
+                        float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
+                        float maxFromPipe = pipeCellMass / dt;
+                        if (mdotFinal < -maxFromPipe)
+                            mdotFinal = -maxFromPipe;
+                    }
+
+                    d.CurrentMdot = mdotFinal;
+
+                    // Limit outflow from cylinder into pipe (positive mdot = volume → pipe)
+                    if (mdotFinal > 0f && d.VolumePort?.Owner is Cylinder cyl)
+                    {
+                        float maxOut = cyl.Mass / dt;
+                        if (mdotFinal > maxOut)
+                            mdotFinal = maxOut;
+                    }
+                }
+
+                // ---- Air fraction for ghost ----
+                if (mdotFinal >= 0)
+                    airFracGhost = volAF;
+                else
+                {
+                    airFracGhost = pipeAF;
+                    if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
+                }
+
+                // ---- Sign convention for velocity ----
+                if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace;
+                else if (mdotFinal >= 0 && !d.IsLeftEnd) uFace = -uFace;
+                else if (mdotFinal < 0 && d.IsLeftEnd) uFace = -uFace;
+                else if (mdotFinal < 0 && !d.IsLeftEnd) uFace = +uFace;
+
+                // ---- Set ghost cells ----
+                if (d.IsLeftEnd)
+                    _pipeSystem.SetGhostLeft(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
+                else
+                    _pipeSystem.SetGhostRight(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
+
+                // ---- Update volume port ----
+                if (d.VolumePort != null)
+                {
+                    d.VolumePort.MassFlowRate = -mdotFinal;
+
+                    if (-mdotFinal >= 0)   // mass entering volume (out of pipe)
+                    {
+                        float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
+                        d.VolumePort.SpecificEnthalpy = pipeH;
+                    }
+                    else   // mass leaving volume (into pipe)
+                    {
+                        d.VolumePort.SpecificEnthalpy = volH;
+                    }
+                }
+            }
+        }
+
+        // ---------- Resolve open ends ----------
+        public void ResolveOpenEnds(float dt)
+        {
+            for (int i = 0; i < OpenEndCount; i++)
+            {
+                ref var d = ref _openEnds[i];
+
+                var (rhoInt, uInt, pInt) = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorStateRight(d.PipeIndex);
+                float afInt = d.IsLeftEnd
+                    ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
+                    : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
+
+                float gamma = d.Gamma;
+                float gm1 = gamma - 1f;
+                float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f));
+                float pAmb = d.AmbientPressure;
+
+                // Characteristic solution (isentropic expansion to ambient)
+                float Jplus = uInt + 2f * cInt / gm1;
+                float Jminus = uInt - 2f * cInt / gm1;
+                float s = pInt / MathF.Pow(rhoInt, gamma);
+                float rhoIso = MathF.Pow(pAmb / s, 1f / gamma);
+                float cIso = MathF.Sqrt(gamma * pAmb / MathF.Max(rhoIso, 1e-12f));
+                float uIso = d.IsLeftEnd
+                    ? (Jminus + 2f * cIso / gm1)
+                    : (Jplus - 2f * cIso / gm1);
+
+                // Supersonic check
+                bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
+                if (!supersonic)
+                {
+                    supersonic = d.IsLeftEnd ? (uIso <= -cIso) : (uIso >= cIso);
+                }
+
+                float rhoGhost, uGhost, pGhost, afGhost;
+                if (supersonic)
+                {
+                    rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt;
+                }
+                else
+                {
+                    rhoGhost = rhoIso; uGhost = uIso; pGhost = pAmb;
+                    bool inflow = d.IsLeftEnd ? (uIso >= 0f) : (uIso <= 0f);
+                    afGhost = inflow ? 1f : afInt;
+                }
+
+                // ------- Mass flow limiter -------
+                int adjCell = d.IsLeftEnd
+                    ? _pipeSystem.GetPipeStart(d.PipeIndex)
+                    : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
+                float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
+                float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
+                float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
+                float cellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
+
+                float area = d.PipeArea;
+                float mdotRaw = rhoGhost * uGhost * area;   // positive out of pipe
+                if (d.IsLeftEnd) mdotRaw = -mdotRaw;        // now positive = out of pipe
+
+                // Outflow limit
+                if (mdotRaw > 0 && mdotRaw * dt > cellMass)
+                {
+                    mdotRaw = cellMass / dt;
+                }
+                // Inflow limit (allow up to 10× cell mass to avoid starving the pipe)
+                else if (mdotRaw < 0 && -mdotRaw * dt > 10f * cellMass)
+                {
+                    mdotRaw = -10f * cellMass / dt;
+                }
+
+                // Recompute uGhost from the limited mdot, keeping rhoGhost, pGhost
+                float mdotMag = MathF.Abs(mdotRaw);
+                uGhost = mdotMag / MathF.Max(rhoGhost * area, 1e-12f);
+                if (d.IsLeftEnd)
+                    uGhost = (mdotRaw >= 0f) ? -uGhost : uGhost;
+                else
+                    uGhost = (mdotRaw >= 0f) ? uGhost : -uGhost;
+
+                // Apply ghost
+                if (d.IsLeftEnd)
+                    _pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
+                else
+                    _pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
+
+                d.LastMassFlowRate = mdotRaw;
+                d.LastFacePressure = pGhost;
+            }
+        }
+    }
+}

Core/Constants.cs

@@ -0,0 +1,11 @@
+namespace FluidSim.Core
+{
+    public static class Constants
+    {
+        public const float Gamma = 1.4f;
+        public const float R_gas = 287f;
+        public const float P_amb = 101325f;
+        public const float T_amb = 300f;
+        public static readonly float Rho_amb = P_amb / (R_gas * T_amb);
+    }
+}

Core/GhostBuffer.cs

@@ -0,0 +1,27 @@
+namespace FluidSim.Core
+{
+    public class GhostBuffer
+    {
+        public float[] Rho, U, P, Y;
+        public int PipeCount { get; }
+
+        public GhostBuffer(int pipeCount)
+        {
+            PipeCount = pipeCount;
+            int size = pipeCount * 2;
+            Rho = new float[size];
+            U = new float[size];
+            P = new float[size];
+            Y = new float[size];
+        }
+
+        public void Set(int pipeIndex, bool isLeftEnd, float rho, float u, float p, float y)
+        {
+            int idx = pipeIndex * 2 + (isLeftEnd ? 0 : 1);
+            Rho[idx] = rho;
+            U[idx] = u;
+            P[idx] = p;
+            Y[idx] = y;
+        }
+    }
+}

Core/IsentropicOrifice.cs

@@ -0,0 +1,31 @@
+using System;
+
+namespace FluidSim.Core
+{
+    public static class IsentropicOrifice
+    {
+        public static void Compute(
+            float pUp, float rhoUp, float TUp,
+            float pDown, float gamma, float R, float area, float Cd,
+            out float mdot, out float rhoFace, out float uFace, out float pFace)
+        {
+            mdot = 0f; rhoFace = rhoUp; uFace = 0f; pFace = pUp;
+            if (area <= 0f || pUp <= 0f || rhoUp <= 0f || TUp <= 0f) return;
+
+            float pr = MathF.Min(pDown / pUp, 1f);
+            if (pr < 1e-6f) pr = 1e-6f;
+            float prCrit = MathF.Pow(2f / (gamma + 1f), gamma / (gamma - 1f));
+            if (pr < prCrit) pr = prCrit;
+
+            float exponent = (gamma - 1f) / gamma;
+            float M = MathF.Sqrt((2f / (gamma - 1f)) * (MathF.Pow(pr, -exponent) - 1f));
+            if (float.IsNaN(M)) M = 0f;
+
+            float aUp = MathF.Sqrt(gamma * R * TUp);
+            uFace = M * aUp;
+            rhoFace = rhoUp * MathF.Pow(pr, 1f / gamma);
+            pFace = pUp * pr;
+            mdot = rhoFace * uFace * area * Cd;
+        }
+    }
+}

Core/OrificeBoundary.cs

@@ -1,99 +0,0 @@
-using FluidSim.Components;
-
-namespace FluidSim.Core
-{
-    public static class OrificeBoundary
-    {
-        public static double MassFlow(double pA, double rhoA, double pB, double rhoB,
-                                       Connection conn)
-        {
-            if (double.IsNaN(pA) || double.IsNaN(rhoA) || double.IsNaN(pB) || double.IsNaN(rhoB) ||
-                double.IsInfinity(pA) || double.IsInfinity(rhoA) || double.IsInfinity(pB) || double.IsInfinity(rhoB) ||
-                pA <= 0 || rhoA <= 0 || pB <= 0 || rhoB <= 0)
-                return 0.0;
-
-            double dp = pA - pB;
-            double sign = Math.Sign(dp);
-            double absDp = Math.Abs(dp);
-            double rhoUp = dp >= 0 ? rhoA : rhoB;
-            double pUp = dp >= 0 ? pA : pB;
-            double pDown = dp >= 0 ? pB : pA;
-            double delta = 1e-6 * pUp;
-
-            if (absDp < delta)
-            {
-                double k = conn.DischargeCoefficient * conn.Area * Math.Sqrt(2 * rhoUp / delta);
-                return k * dp;
-            }
-            else
-            {
-                double pr = pDown / pUp;
-                double choked = Math.Pow(2.0 / (conn.Gamma + 1.0), conn.Gamma / (conn.Gamma - 1.0));
-                if (pr < choked)
-                {
-                    double term = Math.Sqrt(conn.Gamma *
-                        Math.Pow(2.0 / (conn.Gamma + 1.0), (conn.Gamma + 1.0) / (conn.Gamma - 1.0)));
-                    double flow = conn.DischargeCoefficient * conn.Area *
-                                  Math.Sqrt(rhoUp * pUp) * term;
-                    return sign * flow;
-                }
-                else
-                {
-                    double ex = 1.0 - Math.Pow(pr, (conn.Gamma - 1.0) / conn.Gamma);
-                    double flow = conn.DischargeCoefficient * conn.Area *
-                        Math.Sqrt(2.0 * rhoUp * pUp * (conn.Gamma / (conn.Gamma - 1.0)) *
-                                  pr * pr * ex);
-                    return sign * flow;
-                }
-            }
-        }
-
-        public static void PipeVolumeFlux(double pPipe, double rhoPipe, double uPipe,
-                                          double pVol, double rhoVol, double uVol,
-                                          Connection conn, double pipeArea,
-                                          bool isLeftBoundary,
-                                          out double massFlux, out double momFlux, out double energyFlux)
-        {
-            // ----- Compute STAGNATION pressures -----
-            double pStagPipe = pPipe + 0.5 * rhoPipe * uPipe * uPipe;
-            double pStagVol = pVol + 0.5 * rhoVol * uVol * uVol;   // uVol is always 0 for your volumes
-
-            // Mass flow driven by stagnation pressure difference (positive = pipe→volume)
-            double mdot = MassFlow(pStagPipe, rhoPipe, pStagVol, rhoVol, conn);
-
-            // Limit mass flow to the amount that can leave/enter the pipe cell
-            double maxMdot = rhoPipe * pipeArea * 343.0;
-            if (Math.Abs(mdot) > maxMdot) mdot = Math.Sign(mdot) * maxMdot;
-
-            bool flowLeavesPipe = mdot > 0;   // pipe → volume
-
-            double uFace, pFace, rhoFace;
-            double massFluxPerArea;
-
-            if (isLeftBoundary)
-            {
-                massFluxPerArea = -mdot / pipeArea;
-                if (flowLeavesPipe)
-                { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
-                else
-                { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
-            }
-            else // right boundary
-            {
-                massFluxPerArea = mdot / pipeArea;
-                if (flowLeavesPipe)
-                { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; }
-                else
-                { uFace = uVol; pFace = pVol; rhoFace = rhoVol; }
-            }
-
-            // Total enthalpy of the injected fluid
-            double specificEnthalpy = (1.4 / (1.4 - 1.0)) * pFace / Math.Max(rhoFace, 1e-12);
-            double totalEnthalpy = specificEnthalpy + 0.5 * uFace * uFace;
-
-            massFlux = massFluxPerArea;
-            momFlux = massFluxPerArea * uFace + pFace;
-            energyFlux = massFluxPerArea * totalEnthalpy;
-        }
-    }
-}

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;    // direct sound unchanged
+        public float EarlyMix { get; set; } = 0.12f;   // very little early reflection (ground bounce)
+        public float TailMix { get; set; } = 0.18f;   // subtle diffuse tail
+        public float Feedback { get; set; } = 0.35f;   // lower feedback – outdoor doesn't ring
+        public float DampingFreq { get; set; } = 2500f;   // air absorption – high frequencies die quickly
+
+        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 MathF.Tanh((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/Pipesystem.cs

@@ -0,0 +1,693 @@
+using System;
+using System.Diagnostics;
+using System.Numerics;
+
+namespace FluidSim.Core
+{
+    public class PipeSystem
+    {
+        // ---------- Master arrays ----------
+        private float[] _rho, _rhou, _E, _Y;
+        private readonly float[] _area;
+        private readonly float[] _dx;
+        private readonly int[] _pipeStart;
+        private readonly int[] _pipeEnd;
+        private readonly int _totalCells;          // original cell count (visible)
+        private readonly int _allCells;            // total allocated (padded to Vector<float>.Count)
+        private readonly int _pipeCount;
+
+        // Derived state – _p is kept for visualization
+        private float[] _p;
+
+        // Flux arrays for faces INTERNAL to a single pipe (size = _allCells + 1)
+        // Only valid for faces that are NOT pipe boundaries.
+        private float[] _fluxM, _fluxP, _fluxE, _fluxY;
+
+        // Per‑pipe boundary flux buffers (size = _pipeCount)
+        private float[] _leftFluxM, _leftFluxP, _leftFluxE, _leftFluxY;
+        private float[] _rightFluxM, _rightFluxP, _rightFluxE, _rightFluxY;
+
+        // Damping and relaxation
+        private float[] _dampingFactors;
+        private float[] _relaxFactors;
+        private bool _applyDamping;
+        private bool _applyRelax;
+
+        // Ghost buffer (per‑pipe ghost states)
+        private readonly GhostBuffer _ghost;
+
+        // Precomputed flag: true if a face is a pipe boundary (start or end)
+        private readonly bool[] _isPipeBoundaryFace;
+
+        // ---------- Physical constants ----------
+        private const float Gamma = 1.4f;
+        private const float Gm1 = 0.4f;
+        private const float Gm1Inv = 1f / Gm1;          // 2.5
+        private const float GammaOverGm1 = Gamma / Gm1; // 3.5
+        private float _coeffBase;
+        private float _relaxRate;
+        private float _ambientPressure = 101325f;
+        private float _ambientEnergyRef;
+
+        public float DampingMultiplier
+        {
+            set
+            {
+                _coeffBase = 0.1f * value;
+                _applyDamping = _coeffBase != 0f;
+            }
+        }
+        public float EnergyRelaxationRate
+        {
+            set
+            {
+                _relaxRate = value;
+                _applyRelax = _relaxRate != 0f;
+            }
+        }
+        public float AmbientPressure
+        {
+            set
+            {
+                _ambientPressure = value;
+                _ambientEnergyRef = value * Gm1Inv;
+            }
+        }
+
+        // ---------- Profiling ----------
+        public bool EnableProfiling { get; set; }
+        private long _profFluxTicks;
+        private long _profUpdateTicks;
+        private long _profCallCount;
+
+        // ---------- Construction ----------
+        public PipeSystem(int totalCells, int[] pipeStart, int[] pipeEnd,
+                          float[] area, float[] dx,
+                          float initialRho, float initialU, float initialP)
+        {
+            _pipeStart = pipeStart;
+            _pipeEnd = pipeEnd;
+            _pipeCount = pipeStart.Length;
+            _totalCells = totalCells;
+            _area = area;
+            _dx = dx;
+
+            // Pad to SIMD width so all vectorized loops cover the whole data
+            int vecSize = Vector<float>.Count;
+            _allCells = totalCells % vecSize == 0 ? totalCells : totalCells + vecSize - (totalCells % vecSize);
+
+            _rho = new float[_allCells];
+            _rhou = new float[_allCells];
+            _E = new float[_allCells];
+            _Y = new float[_allCells];
+            _p = new float[_allCells];          // pressure for drawing
+            int faceCount = _allCells + 1;
+            _fluxM = new float[faceCount];
+            _fluxP = new float[faceCount];
+            _fluxE = new float[faceCount];
+            _fluxY = new float[faceCount];
+
+            // Per‑pipe boundary flux buffers
+            _leftFluxM = new float[_pipeCount];
+            _leftFluxP = new float[_pipeCount];
+            _leftFluxE = new float[_pipeCount];
+            _leftFluxY = new float[_pipeCount];
+            _rightFluxM = new float[_pipeCount];
+            _rightFluxP = new float[_pipeCount];
+            _rightFluxE = new float[_pipeCount];
+            _rightFluxY = new float[_pipeCount];
+
+            _dampingFactors = new float[_allCells];
+            _relaxFactors = new float[_allCells];
+            _applyDamping = _coeffBase != 0f;
+            _applyRelax = _relaxRate != 0f;
+
+            _ghost = new GhostBuffer(_pipeCount);
+            _ambientEnergyRef = initialP * Gm1Inv;
+
+            // Mark faces that coincide with a pipe boundary (start or end)
+            _isPipeBoundaryFace = new bool[faceCount];
+            for (int p = 0; p < _pipeCount; p++)
+            {
+                _isPipeBoundaryFace[_pipeStart[p]] = true;
+                _isPipeBoundaryFace[_pipeEnd[p]] = true;
+            }
+
+            // Initialize uniform state
+            float initE = initialP / (Gm1 * initialRho);
+            float rhoE = initialRho * initE + 0.5f * initialRho * initialU * initialU;
+            for (int i = 0; i < totalCells; i++)
+            {
+                _rho[i] = initialRho;
+                _rhou[i] = initialRho * initialU;
+                _E[i] = rhoE;
+                _Y[i] = 1f;
+            }
+        }
+
+        // ---------- Ghost setters (for BoundarySystem) ----------
+        public void SetGhostLeft(int pipeIndex, float rho, float u, float p, float y)
+            => _ghost.Set(pipeIndex, true, rho, u, p, y);
+        public void SetGhostRight(int pipeIndex, float rho, float u, float p, float y)
+            => _ghost.Set(pipeIndex, false, rho, u, p, y);
+
+        // ---------- Public read methods ----------
+        public int TotalCells => _totalCells;
+        public int PipeCount => _pipeCount;
+        public int GetPipeStart(int pipeIdx) => _pipeStart[pipeIdx];
+        public int GetPipeEnd(int pipeIdx) => _pipeEnd[pipeIdx];
+        public float GetCellPressure(int i) => _p[i];
+        public float GetCellDensity(int i) => _rho[i];
+        public float GetCellDx(int i) => _dx[i];
+        public float GetCellArea(int i) => _area[i];
+        public float GetCellVelocity(int i)
+        {
+            float rho = _rho[i];
+            return rho > 1e-12f ? _rhou[i] / rho : 0f;
+        }
+        public float GetCellAirFraction(int i) => _Y[i];
+
+        public (float rho, float u, float p) GetInteriorStateLeft(int pipeIdx)
+        {
+            int i = _pipeStart[pipeIdx];
+            float rho = _rho[i];
+            float rhou = _rhou[i];
+            float u = rhou / MathF.Max(rho, 1e-12f);
+            float p = Gm1 * (_E[i] - 0.5f * rhou * u);
+            return (rho, u, p);
+        }
+        public (float rho, float u, float p) GetInteriorStateRight(int pipeIdx)
+        {
+            int i = _pipeEnd[pipeIdx] - 1;
+            float rho = _rho[i];
+            float rhou = _rhou[i];
+            float u = rhou / MathF.Max(rho, 1e-12f);
+            float p = Gm1 * (_E[i] - 0.5f * rhou * u);
+            return (rho, u, p);
+        }
+        public float GetInteriorAirFractionLeft(int pipeIdx) => _Y[_pipeStart[pipeIdx]];
+        public float GetInteriorAirFractionRight(int pipeIdx) => _Y[_pipeEnd[pipeIdx] - 1];
+
+        public void SetCellState(int i, float rho, float u, float p, float y = 1f)
+        {
+            if (i < 0 || i >= _totalCells) return;
+            _rho[i] = rho;
+            _rhou[i] = rho * u;
+            _E[i] = p * Gm1Inv + 0.5f * rho * u * u;
+            _Y[i] = y;
+        }
+
+        // ---------- Main step ----------
+        public void SimulateStep(float dt)
+        {
+            long t0 = 0, t1 = 0;
+            if (EnableProfiling)
+            {
+                _profCallCount++;
+                t0 = Stopwatch.GetTimestamp();
+            }
+
+            ComputeFluxes(dt);
+
+            if (EnableProfiling)
+            {
+                t1 = Stopwatch.GetTimestamp();
+                _profFluxTicks += (t1 - t0);
+                t0 = t1;
+            }
+
+            UpdateCells(dt);
+
+            if (EnableProfiling)
+            {
+                t1 = Stopwatch.GetTimestamp();
+                _profUpdateTicks += (t1 - t0);
+            }
+        }
+
+        // ---------- Flux computation ----------
+        private void ComputeFluxes(float dt)
+        {
+            float fm, fp, fe;
+            int vecSize = Vector<float>.Count;
+
+            // ---- 1. Left ghost boundaries → per‑pipe buffers ----
+            for (int p = 0; p < _pipeCount; p++)
+            {
+                int idx = _pipeStart[p];
+                int ghostIdx = p * 2;
+                float rL = _ghost.Rho[ghostIdx];
+                float uL = _ghost.U[ghostIdx];
+                float pL = _ghost.P[ghostIdx];
+                float YL = _ghost.Y[ghostIdx];
+                float cL = MathF.Sqrt(Gamma * pL / MathF.Max(rL, 1e-12f));
+
+                float rR = _rho[idx], rhouR = _rhou[idx];
+                float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
+                float uR = rhouR * invRhoR;
+                float pR = Gm1 * (_E[idx] - 0.5f * rhouR * uR);
+                float cR = MathF.Sqrt(Gamma * pR * invRhoR);
+                float YR = _Y[idx];
+
+                LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
+                _leftFluxM[p] = fm; _leftFluxP[p] = fp; _leftFluxE[p] = fe;
+
+                float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
+                ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
+                _leftFluxY[p] = fy;
+            }
+
+            // ---- 2. Right ghost boundaries → per‑pipe buffers ----
+            for (int p = 0; p < _pipeCount; p++)
+            {
+                int idx = _pipeEnd[p] - 1;
+                int ghostIdx = p * 2 + 1;
+                float rR = _ghost.Rho[ghostIdx];
+                float uR = _ghost.U[ghostIdx];
+                float pR = _ghost.P[ghostIdx];
+                float YR = _ghost.Y[ghostIdx];
+                float cR = MathF.Sqrt(Gamma * pR / MathF.Max(rR, 1e-12f));
+
+                float rL = _rho[idx], rhouL = _rhou[idx];
+                float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
+                float uL = rhouL * invRhoL;
+                float pL = Gm1 * (_E[idx] - 0.5f * rhouL * uL);
+                float cL = MathF.Sqrt(Gamma * pL * invRhoL);
+                float YL = _Y[idx];
+
+                LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
+                _rightFluxM[p] = fm; _rightFluxP[p] = fp; _rightFluxE[p] = fe;
+
+                float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
+                ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
+                _rightFluxY[p] = fy;
+            }
+
+            // ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
+            for (int face = 1; face < _totalCells; face++)
+            {
+                // Skip faces that belong to a pipe boundary (they are already handled)
+                if (_isPipeBoundaryFace[face])
+                    continue;
+
+                // Try to vectorize a block of contiguous non‑boundary faces
+                if (face + vecSize - 1 < _totalCells)
+                {
+                    bool canVectorize = true;
+                    for (int f = face; f < face + vecSize; f++)
+                    {
+                        if (_isPipeBoundaryFace[f])
+                        {
+                            canVectorize = false;
+                            break;
+                        }
+                    }
+
+                    if (canVectorize)
+                    {
+                        // --- Vectorised block ---
+                        var rhoL = new Vector<float>(_rho, face - 1);
+                        var rhouL = new Vector<float>(_rhou, face - 1);
+                        var EL = new Vector<float>(_E, face - 1);
+                        var YL = new Vector<float>(_Y, face - 1);
+                        var rhoR = new Vector<float>(_rho, face);
+                        var rhouR = new Vector<float>(_rhou, face);
+                        var ER = new Vector<float>(_E, face);
+                        var YR = new Vector<float>(_Y, face);
+
+                        var invRhoL = Vector<float>.One / Vector.Max(rhoL, new Vector<float>(1e-12f));
+                        var invRhoR = Vector<float>.One / Vector.Max(rhoR, new Vector<float>(1e-12f));
+                        var uL = rhouL * invRhoL;
+                        var uR = rhouR * invRhoR;
+                        var kinL = 0.5f * rhouL * uL;
+                        var kinR = 0.5f * rhouR * uR;
+                        var pL = Gm1 * (EL - kinL);
+                        var pR = Gm1 * (ER - kinR);
+                        var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
+                        var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
+
+                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
+                        var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
+
+                        var FmL = rhoL * uL;
+                        var FpL = rhoL * uL * uL + pL;
+                        var FeL = (rhoL * ELs + pL) * uL;
+
+                        var FmR = rhoR * uR;
+                        var FpR = rhoR * uR * uR + pR;
+                        var FeR = (rhoR * ERs + pR) * uR;
+
+                        var absUL = Vector.Abs(uL);
+                        var absUR = Vector.Abs(uR);
+                        var alpha = Vector.Max(absUL + cL, absUR + cR);
+
+                        var fmVec = 0.5f * (FmL + FmR) - 0.5f * alpha * (rhoR - rhoL);
+                        var fpVec = 0.5f * (FpL + FpR) - 0.5f * alpha * (rhouR - rhouL);
+                        var feVec = 0.5f * (FeL + FeR) - 0.5f * alpha * (rhoR * ERs - rhoL * ELs);
+
+                        var fyL = FmL * YL;
+                        var fyR = FmR * YR;
+                        var fyVec = 0.5f * (fyL + fyR) - 0.5f * alpha * (rhoR * YR - rhoL * YL);
+
+                        fmVec.CopyTo(_fluxM, face);
+                        fpVec.CopyTo(_fluxP, face);
+                        feVec.CopyTo(_fluxE, face);
+                        fyVec.CopyTo(_fluxY, face);
+
+                        face += vecSize - 1; // loop increment will add 1
+                        continue;
+                    }
+                }
+
+                // --- Scalar fallback for a single interior face ---
+                {
+                    int iL = face - 1, iR = face;
+                    float rL = _rho[iL], rhouL = _rhou[iL];
+                    float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
+                    float uL = rhouL * invRhoL;
+                    float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
+                    float cL = MathF.Sqrt(Gamma * pL * invRhoL);
+                    float YL = _Y[iL];
+
+                    float rR = _rho[iR], rhouR = _rhou[iR];
+                    float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
+                    float uR = rhouR * invRhoR;
+                    float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
+                    float cR = MathF.Sqrt(Gamma * pR * invRhoR);
+                    float YR = _Y[iR];
+
+                    LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
+                    _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
+
+                    float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
+                    ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
+                    _fluxY[face] = fy;
+                }
+            }
+        }
+
+        // ---------- Cell update (per pipe, using correct boundary fluxes) ----------
+        private void UpdateCells(float dt)
+        {
+            int vecSize = Vector<float>.Count;
+            float dtRelax = -_relaxRate * dt;
+
+            // Precompute damping and relaxation factors globally
+            if (_applyDamping)
+            {
+                for (int i = 0; i < _totalCells; i++)
+                {
+                    float rho = _rho[i];
+                    _dampingFactors[i] = rho > 1e-12f
+                        ? MathF.Exp(-_coeffBase * dt / rho)
+                        : 1f;
+                }
+            }
+            if (_applyRelax)
+            {
+                float relaxVal = MathF.Exp(dtRelax);
+                for (int i = 0; i < _totalCells; i++)
+                    _relaxFactors[i] = relaxVal;
+            }
+
+            // Update each pipe separately
+            for (int p = 0; p < _pipeCount; p++)
+            {
+                int start = _pipeStart[p];
+                int end = _pipeEnd[p];    // exclusive
+                int len = end - start;
+                if (len == 0) continue;
+
+                // ------- Left boundary cell (i = start) ------
+                {
+                    int i = start;
+                    float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
+
+                    // left face: always the pipe's left boundary flux
+                    float fluxM_L = _leftFluxM[p];
+                    float fluxP_L = _leftFluxP[p];
+                    float fluxE_L = _leftFluxE[p];
+                    float fluxY_L = _leftFluxY[p];
+
+                    // right face: depends on pipe length
+                    float fluxM_R, fluxP_R, fluxE_R, fluxY_R;
+                    if (len == 1)
+                    {
+                        // Only one cell: right face is the pipe's right boundary flux
+                        fluxM_R = _rightFluxM[p];
+                        fluxP_R = _rightFluxP[p];
+                        fluxE_R = _rightFluxE[p];
+                        fluxY_R = _rightFluxY[p];
+                    }
+                    else
+                    {
+                        // interior face (global flux at index i+1)
+                        fluxM_R = _fluxM[i + 1];
+                        fluxP_R = _fluxP[i + 1];
+                        fluxE_R = _fluxE[i + 1];
+                        fluxY_R = _fluxY[i + 1];
+                    }
+
+                    float dtdx = dt / _dx[i];
+                    float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                    float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                    float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                    float rhoYOld = rhoOld * YOld;
+                    float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                    if (_applyDamping) rhouNew *= _dampingFactors[i];
+                    if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
+
+                    rhoNew = MathF.Max(rhoNew, 1e-12f);
+                    float kin = 0.5f * rhouNew * rhouNew / rhoNew;
+                    float eMin = 100f * Gm1Inv + kin;
+                    ENew = MathF.Max(ENew, eMin);
+
+                    _rho[i] = rhoNew;
+                    _rhou[i] = rhouNew;
+                    _E[i] = ENew;
+                    _Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
+                }
+
+                // ------- Interior cells (i = start+1 to end-2) ------
+                if (len > 2)
+                {
+                    int iCell = start + 1;
+                    int iEnd = end - 1; // exclusive upper bound
+
+                    // Vectorised path for interior cells (if available)
+                    for (; iCell <= iEnd - vecSize; iCell += vecSize)
+                    {
+                        var rhoOld = new Vector<float>(_rho, iCell);
+                        var rhouOld = new Vector<float>(_rhou, iCell);
+                        var EOld = new Vector<float>(_E, iCell);
+                        var YOld = new Vector<float>(_Y, iCell);
+
+                        var fluxM_L = new Vector<float>(_fluxM, iCell);
+                        var fluxP_L = new Vector<float>(_fluxP, iCell);
+                        var fluxE_L = new Vector<float>(_fluxE, iCell);
+                        var fluxY_L = new Vector<float>(_fluxY, iCell);
+
+                        var fluxM_R = new Vector<float>(_fluxM, iCell + 1);
+                        var fluxP_R = new Vector<float>(_fluxP, iCell + 1);
+                        var fluxE_R = new Vector<float>(_fluxE, iCell + 1);
+                        var fluxY_R = new Vector<float>(_fluxY, iCell + 1);
+
+                        var dtdx = new Vector<float>(dt) / new Vector<float>(_dx, iCell);
+
+                        var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                        var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                        var ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                        var rhoYOld = rhoOld * YOld;
+                        var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                        if (_applyDamping)
+                            rhouNew *= new Vector<float>(_dampingFactors, iCell);
+                        if (_applyRelax)
+                        {
+                            var ambRef = new Vector<float>(_ambientEnergyRef);
+                            var relax = new Vector<float>(_relaxFactors, iCell);
+                            ENew = ambRef + (ENew - ambRef) * relax;
+                        }
+
+                        rhoNew = Vector.Max(rhoNew, new Vector<float>(1e-12f));
+                        var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
+                        var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
+                        ENew = Vector.Max(ENew, eMin);
+
+                        rhoNew.CopyTo(_rho, iCell);
+                        rhouNew.CopyTo(_rhou, iCell);
+                        ENew.CopyTo(_E, iCell);
+                        var yNew = rhoYNew / rhoNew;
+                        yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
+                        yNew.CopyTo(_Y, iCell);
+                    }
+
+                    // Scalar remainder for interior cells
+                    for (; iCell < iEnd; iCell++)
+                    {
+                        float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
+                        float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
+                        float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
+                        float dtdx = dt / _dx[iCell];
+
+                        float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                        float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                        float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                        float rhoYOld = rhoOld * YOld;
+                        float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                        if (_applyDamping) rhouNew *= _dampingFactors[iCell];
+                        if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell];
+
+                        rhoNew = MathF.Max(rhoNew, 1e-12f);
+                        float kin = 0.5f * rhouNew * rhouNew / rhoNew;
+                        float eMin = 100f * Gm1Inv + kin;
+                        ENew = MathF.Max(ENew, eMin);
+
+                        _rho[iCell] = rhoNew;
+                        _rhou[iCell] = rhouNew;
+                        _E[iCell] = ENew;
+                        _Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
+                    }
+                }
+
+                // ------- Right boundary cell (i = end-1, if len > 1) ------
+                if (len > 1)
+                {
+                    int i = end - 1;
+                    float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
+
+                    // left face
+                    float fluxM_L, fluxP_L, fluxE_L, fluxY_L;
+                    if (len == 2)
+                    {
+                        // Only two cells: left face is the pipe's left boundary flux
+                        fluxM_L = _leftFluxM[p];
+                        fluxP_L = _leftFluxP[p];
+                        fluxE_L = _leftFluxE[p];
+                        fluxY_L = _leftFluxY[p];
+                    }
+                    else
+                    {
+                        // interior face (global flux at i)
+                        fluxM_L = _fluxM[i];
+                        fluxP_L = _fluxP[i];
+                        fluxE_L = _fluxE[i];
+                        fluxY_L = _fluxY[i];
+                    }
+
+                    // right face: always the pipe's right boundary flux
+                    float fluxM_R = _rightFluxM[p];
+                    float fluxP_R = _rightFluxP[p];
+                    float fluxE_R = _rightFluxE[p];
+                    float fluxY_R = _rightFluxY[p];
+
+                    float dtdx = dt / _dx[i];
+                    float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
+                    float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
+                    float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
+                    float rhoYOld = rhoOld * YOld;
+                    float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
+
+                    if (_applyDamping) rhouNew *= _dampingFactors[i];
+                    if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
+
+                    rhoNew = MathF.Max(rhoNew, 1e-12f);
+                    float kin = 0.5f * rhouNew * rhouNew / rhoNew;
+                    float eMin = 100f * Gm1Inv + kin;
+                    ENew = MathF.Max(ENew, eMin);
+
+                    _rho[i] = rhoNew;
+                    _rhou[i] = rhouNew;
+                    _E[i] = ENew;
+                    _Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
+                }
+            }
+
+            // Recompute pressure for all cells (for visualization)
+            for (int i = 0; i < _totalCells; i++)
+            {
+                float rho = _rho[i];
+                float rhou = _rhou[i];
+                float u = rhou / MathF.Max(rho, 1e-12f);
+                _p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
+            }
+        }
+
+        // ---------- Scalar flux helpers ----------
+        private static void LaxFlux(float rL, float uL, float pL, float cL,
+                                    float rR, float uR, float pR, float cR,
+                                    out float fm, out float fp, out float fe)
+        {
+            float EL = pL * Gm1Inv / rL + 0.5f * uL * uL;
+            float ER = pR * Gm1Inv / rR + 0.5f * uR * uR;
+            float FmL = rL * uL;
+            float FpL = rL * uL * uL + pL;
+            float FeL = (rL * EL + pL) * uL;
+            float FmR = rR * uR;
+            float FpR = rR * uR * uR + pR;
+            float FeR = (rR * ER + pR) * uR;
+            float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
+            fm = 0.5f * (FmL + FmR) - 0.5f * alpha * (rR - rL);
+            fp = 0.5f * (FpL + FpR) - 0.5f * alpha * (rR * uR - rL * uL);
+            fe = 0.5f * (FeL + FeR) - 0.5f * alpha * (rR * ER - rL * EL);
+        }
+
+        private static void ScalarFlux(float rL, float uL, float YL,
+                                       float rR, float uR, float YR,
+                                       float alpha, out float fy)
+        {
+            float FyL = rL * uL * YL;
+            float FyR = rR * uR * YR;
+            fy = 0.5f * (FyL + FyR) - 0.5f * alpha * (rR * YR - rL * YL);
+        }
+
+        public int GetRequiredSubSteps(float dtGlobal, float cflTarget = 0.8f)
+        {
+            float maxW = 0f;
+            for (int i = 0; i < _totalCells; i++)
+            {
+                float rho = MathF.Max(_rho[i], 1e-12f);
+                float u = MathF.Abs(_rhou[i] / rho);
+                float p = Gm1 * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
+                float c = MathF.Sqrt(Gamma * p / rho);
+                float w = u + c;
+                if (w > maxW) maxW = w;
+            }
+            maxW = MathF.Max(maxW, 1e-8f);
+            float minDx = _dx.Min();  // need using System.Linq;
+            return Math.Max(1, (int)MathF.Ceiling(dtGlobal * maxW / (cflTarget * minDx)));
+        }
+
+        // ---------- Profiling report ----------
+        public string GetProfileReport()
+        {
+            if (!EnableProfiling || _profCallCount == 0)
+                return "Pipe profiling disabled or no data.";
+
+            double freq = Stopwatch.Frequency;
+            long totalTicks = _profFluxTicks + _profUpdateTicks;
+            if (totalTicks == 0) return "No pipe profile data collected.";
+
+            double totalMs = totalTicks * 1000.0 / freq;
+            double avgCallUs = totalMs * 1000.0 / _profCallCount;
+
+            double fluxMs = _profFluxTicks * 1000.0 / freq;
+            double updateMs = _profUpdateTicks * 1000.0 / freq;
+
+            double fluxAvgUs = fluxMs * 1000.0 / _profCallCount;
+            double updateAvgUs = updateMs * 1000.0 / _profCallCount;
+
+            string report = $"  Pipe kernel (over {_profCallCount} calls, total {totalMs:F2} ms, avg {avgCallUs:F2} µs/call):\n";
+            report += $"    Fluxes (incl. primitives): {fluxMs:F2} ms  ({_profFluxTicks * 100.0 / totalTicks:F1}%), avg {fluxAvgUs:F2} µs/call\n";
+            report += $"    Update cells:              {updateMs:F2} ms  ({_profUpdateTicks * 100.0 / totalTicks:F1}%), avg {updateAvgUs:F2} µs/call\n";
+
+            _profFluxTicks = 0;
+            _profUpdateTicks = 0;
+            _profCallCount = 0;
+
+            return report;
+        }
+    }
+}

Core/Solver.cs

@@ -1,4 +1,7 @@
-using FluidSim.Audio;
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
 using FluidSim.Components;
 using FluidSim.Interfaces;
 
@@ -6,120 +9,94 @@ namespace FluidSim.Core
 {
     public class Solver
     {
-        private readonly List<Volume0D> _volumes = new();
+        private readonly List<IComponent> _components = new();
-        private readonly List<Pipe1D> _pipes = new();
+        private PipeSystem _pipeSystem;
-        private readonly List<Connection> _connections = new();
+        private BoundarySystem _boundarySystem;
+        private double _dt;
 
-        public float LastSample { get; private set; }
+        public int SubStepCount { get; set; } = 4;
+        public bool EnableProfiling { get; set; } = false;
 
-        public void AddVolume(Volume0D v) => _volumes.Add(v);
+        private long _stepCount;
-        public void AddPipe(Pipe1D p) => _pipes.Add(p);
+        private long _ticksOrifice, _ticksOpenEnd, _ticksPipe, _ticksUpdate;
-        public void AddConnection(Connection c) => _connections.Add(c);
+
+        public void SetTimeStep(double dt) => _dt = dt;
+        public void AddComponent(IComponent component) => _components.Add(component);
+
+        public void SetPipeSystem(PipeSystem pipeSystem)
+        {
+            _pipeSystem = pipeSystem;
+        }
+        public void SetBoundarySystem(BoundarySystem boundarySystem)
+        {
+            _boundarySystem = boundarySystem;
+        }
 
         public void Step()
         {
-            // 1. Publish volume states to their own ports
+            if (_pipeSystem == null || _boundarySystem == null) return;
-            foreach (var v in _volumes)
-                v.PushStateToPort();
 
-            // 2. Handle direct volume‑to‑volume connections
+            int nSub = _pipeSystem.GetRequiredSubSteps((float)_dt, 0.8f);
-            foreach (var conn in _connections)
+            nSub = Math.Max(nSub, SubStepCount);   // never go below fixed minimum
+            float dtSub = (float)(_dt / nSub);
+
+            for (int sub = 0; sub < nSub; sub++)
             {
-                if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB))
+                long t0;
-                {
-                    Volume0D volA = _volumes.Find(v => v.Port == conn.PortA);
-                    Volume0D volB = _volumes.Find(v => v.Port == conn.PortB);
 
-                    if (volA == null || volB == null) continue;
+                t0 = Stopwatch.GetTimestamp();
+                _boundarySystem.ResolveOrifices(dtSub);
+                _ticksOrifice += Stopwatch.GetTimestamp() - t0;
 
-                    double pA = volA.Pressure, rhoA = volA.Density;
+                t0 = Stopwatch.GetTimestamp();
-                    double pB = volB.Pressure, rhoB = volB.Density;
+                _boundarySystem.ResolveOpenEnds(dtSub);
+                _ticksOpenEnd += Stopwatch.GetTimestamp() - t0;
 
-                    double mdot = OrificeBoundary.MassFlow(pA, rhoA, pB, rhoB, conn);
+                t0 = Stopwatch.GetTimestamp();
-
+                _pipeSystem.SimulateStep(dtSub);
-                    if (mdot > 0) // A → B
+                _ticksPipe += Stopwatch.GetTimestamp() - t0;
-                    {
-                        volA.Port.MassFlowRate = -mdot;
-                        volB.Port.MassFlowRate = mdot;
-                        volB.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // fluid from A
-                        volA.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // outflow carries its own enthalpy
-                    }
-                    else // B → A (mdot negative)
-                    {
-                        volA.Port.MassFlowRate = -mdot;  // positive
-                        volB.Port.MassFlowRate = mdot;  // negative
-                        volA.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // fluid from B
-                        volB.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // outflow carries its own
-                    }
-                }
             }
 
-            // 3. Pipe‑volume boundary conditions – unchanged
+            long tUS = Stopwatch.GetTimestamp();
-            foreach (var conn in _connections)
+            foreach (var comp in _components)
+                comp.UpdateState((float)_dt);
+            _ticksUpdate += Stopwatch.GetTimestamp() - tUS;
+
+            _stepCount++;
+            if (_stepCount % 5000 == 0 && EnableProfiling)
             {
-                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+                double freq = Stopwatch.Frequency;
-                    SetVolumeBC(conn.PortA, conn.PortB);
+                double total = _ticksOrifice + _ticksOpenEnd + _ticksPipe + _ticksUpdate;
-                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
+                double avgStepUs = (total / freq) * 1e6 / 5000.0;
-                    SetVolumeBC(conn.PortB, conn.PortA);
+
+                int orificeCalls = 5000 * nSub;
+                int updateCalls = 5000;
+
+                double orificeMs = _ticksOrifice * 1000.0 / freq;
+                double openEndMs = _ticksOpenEnd * 1000.0 / freq;
+                double pipeMs = _ticksPipe * 1000.0 / freq;
+                double updateMs = _ticksUpdate * 1000.0 / freq;
+
+                double orificeAvgUs = orificeMs * 1000.0 / orificeCalls;
+                double openEndAvgUs = openEndMs * 1000.0 / orificeCalls;
+                double pipeAvgUs = pipeMs * 1000.0 / orificeCalls;
+                double updateAvgUs = updateMs * 1000.0 / updateCalls;
+
+                Console.WriteLine($"--- Solver ({5000} steps, nSub={nSub}) ---");
+                Console.WriteLine($"  Average step: {avgStepUs:F2} µs");
+                Console.WriteLine($"  Orifice:  {orificeMs:F2} ms ({(double)_ticksOrifice / total * 100:F1}%), avg {orificeAvgUs:F2} µs/call");
+                Console.WriteLine($"  OpenEnd:  {openEndMs:F2} ms ({(double)_ticksOpenEnd / total * 100:F1}%), avg {openEndAvgUs:F2} µs/call");
+                Console.WriteLine($"  Pipe:     {pipeMs:F2} ms ({(double)_ticksPipe / total * 100:F1}%), avg {pipeAvgUs:F2} µs/call");
+                Console.WriteLine($"  Update:   {updateMs:F2} ms ({(double)_ticksUpdate / total * 100:F1}%), avg {updateAvgUs:F2} µs/call");
+
+                // Pipe internal breakdown (with per-phase averages)
+                if (_pipeSystem.EnableProfiling)
+                {
+                    Console.WriteLine(_pipeSystem.GetProfileReport());
                 }
 
-            // 4. Run pipe simulations
+                _ticksOrifice = _ticksOpenEnd = _ticksPipe = _ticksUpdate = 0;
-            foreach (var p in _pipes)
-                p.Simulate();
-
-            // 5. Transfer pipe‑to‑volume flows – unchanged
-            foreach (var conn in _connections)
-            {
-                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    TransferPipeToVolume(conn.PortA, conn.PortB);
-                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                    TransferPipeToVolume(conn.PortB, conn.PortA);
             }
-
-            // 6. Integrate volumes
-            foreach (var v in _volumes)
-                v.Integrate();
-
-            // 7. COMPUTE AUDIO SAMPLE from all sound connections
-            double sample = 0f;
-            foreach (var conn in _connections)
-            {
-                if (conn is SoundConnection)
-                {
-                    // Both ports have the same absolute mass flow; either works.
-                    sample += SoundProcessor.ComputeSample(conn);
-                }
-            }
-            LastSample = (float)Math.Tanh(sample);
-        }
-
-        bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
-        bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
-        Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
-
-        void SetVolumeBC(Port pipePort, Port volPort)
-        {
-            Pipe1D pipe = GetPipe(pipePort);
-            if (pipe == null) return;
-            bool isLeft = pipe.PortA == pipePort;
-
-            if (isLeft)
-                pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
-            else
-                pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
-        }
-
-        void TransferPipeToVolume(Port pipePort, Port volPort)
-        {
-            double mdot = pipePort.MassFlowRate;
-            volPort.MassFlowRate = -mdot;
-
-            if (mdot < 0)  // pipe → volume
-            {
-                // pipePort.SpecificEnthalpy is already total (h + ½u²)
-                volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
-            }
-            // else: volume → pipe, volume’s own static enthalpy is used (already set)
         }
     }
 }

Core/SoundProcessor.cs

@@ -0,0 +1,34 @@
+using System;
+
+namespace FluidSim.Core
+{
+    public class SoundProcessor
+    {
+        private readonly float dt;
+        private readonly float scaleFactor;   // 1 / (4π r)
+        private float flowLP, prevMassFlowOut, smoothDMdt;
+        private readonly float lpAlpha, alpha;
+
+        public float Gain = 1f;
+
+        public SoundProcessor(int sampleRate, float listenerDistance = 1f)
+        {
+            dt = 1f / sampleRate;
+            scaleFactor = 1f / (4f * MathF.PI * listenerDistance);
+            float tau = 0.02f;
+            alpha = MathF.Exp(-dt / tau);
+            float tauLP = 0.005f;
+            lpAlpha = MathF.Exp(-dt / tauLP);
+        }
+
+        public float Process(float massFlowOut)
+        {
+            flowLP = lpAlpha * flowLP + (1f - lpAlpha) * massFlowOut;
+            float rawDerivative = (flowLP - prevMassFlowOut) / dt;
+            prevMassFlowOut = flowLP;
+            smoothDMdt = alpha * smoothDMdt + (1f - alpha) * rawDerivative;
+            float pressure = smoothDMdt * scaleFactor * Gain;
+            return MathF.Tanh(pressure);
+        }
+    }
+}

Core/ThreadLoadTracker.cs

@@ -0,0 +1,34 @@
+using System;
+using System.Threading;
+
+namespace FluidSim
+{
+    /// <summary>
+    /// Tracks the duty cycle of a worker thread using an exponential moving average.
+    /// Thread‑safe: one writer (the sim thread), any reader (UI thread).
+    /// </summary>
+    public class ThreadLoadTracker
+    {
+        private double _loadPercent;   // 0 .. 100, accessed with Volatile.Read/Write
+        private const double Alpha = 0.1;   // smoothing factor (higher = faster response)
+
+        /// <summary>
+        /// Update the load percentage with a new observation.
+        /// </summary>
+        /// <param name="busyMs">Time spent on real work in the last cycle.</param>
+        /// <param name="totalMs">Total time of the last cycle (work + idle). If zero, ignored.</param>
+        public void Record(double busyMs, double totalMs)
+        {
+            if (totalMs <= 0) return;
+            double instantLoad = busyMs / totalMs * 100.0;
+
+            // Exponential moving average
+            double old = Volatile.Read(ref _loadPercent);
+            double newLoad = old + Alpha * (instantLoad - old);
+            Volatile.Write(ref _loadPercent, newLoad);
+        }
+
+        /// <summary>Current smoothed load percentage (0‑100).</summary>
+        public double LoadPercent => Volatile.Read(ref _loadPercent);
+    }
+}

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>
 
@@ -13,4 +13,13 @@
       <PackageReference Include="SFML.Net" Version="3.0.0" />
   </ItemGroup>
 
+  <ItemGroup>
+      <None Update="fonts\FiraCodeNerdFont-Medium.ttf">
+        <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+      </None>
+      <None Update="fonts\LiberationMono-Regular.ttf">
+          <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
+      </None>
+  </ItemGroup>
+
 </Project>

Interfaces/IComponent.cs

@@ -0,0 +1,10 @@
+using System.Collections.Generic;
+
+namespace FluidSim.Interfaces
+{
+    public interface IComponent
+    {
+        IReadOnlyList<Port> Ports { get; }
+        void UpdateState(float dt);
+    }
+}

Interfaces/Port.cs

@@ -2,19 +2,23 @@
 {
     public class Port
     {
-        public double Pressure;             // Pa
+        public float MassFlowRate;          // kg/s, positive INTO owning component
-        public double MassFlowRate;         // kg/s, positive INTO the component
+        public float SpecificEnthalpy;      // J/kg
-        public double SpecificEnthalpy;     // J/kg, enthalpy of fluid entering this port
+        public float Pressure;              // Pa
-        public double Density;              // kg/m³
+        public float Density;               // kg/m³
-        public double Temperature;          // K
+        public float Temperature;           // K
+        public float AirFraction;           // mass fraction (0 = exhaust, 1 = air)
+
+        public object? Owner { get; set; }
 
         public Port()
         {
-            Pressure = 101325.0;
+            MassFlowRate = 0f;
-            MassFlowRate = 0.0;
+            SpecificEnthalpy = 0f;
-            SpecificEnthalpy = 0.0;
+            Pressure = 101325f;
-            Density = 1.225;
+            Density = 1.225f;
-            Temperature = 300.0;
+            Temperature = 300f;
+            AirFraction = 1f;
         }
     }
 }

Program.cs

@@ -1,96 +1,216 @@
-using SFML.Graphics;
+using FluidSim.Audio;
-using SFML.Window;
+using FluidSim.Core;
+using FluidSim.Tests;
+using SFML.Graphics;
 using SFML.System;
+using SFML.Window;
+using System;
 using System.Diagnostics;
-using FluidSim.Scenarios;
+using System.Threading;
-using FluidSim.Audio;
+using System.Threading.Tasks;
 
 namespace FluidSim;
 
 public class Program
 {
     private const int SampleRate = 44100;
-    private static volatile bool running = true;
+    private const double DrawFrequency = 60.0;
-    // Global step counter – incremented every simulation step
+
-    private static long stepCount = 0;
+    // Playback speed
+    private static double _desiredSpeed = 0.001;
+    private static double _currentDisplaySpeed = _desiredSpeed;
+    private const double MinSpeed = 0.0001;
+    private const double MaxSpeed = 1.0;
+    private const double ScrollFactor = 1.1;
+    private static double _lastNormalSpeed = 0.1;
+    private static bool _isRealTime = false;
+
+    private static volatile bool _timeWarpActive;
+
+    // Thread load tracking
+    private static ThreadLoadTracker _loadTracker = new ThreadLoadTracker();
+
+    // Audio & simulation
+    private static SimulationRingBuffer _simRingBuffer = null!;
+    private static SoundEngine _soundEngine = null!;
+    private static Scenario _scenario = null!;   // cast to access ThrottleArea
+    private static Font? _overlayFont;
+    private static Text? _overlayText;
+
+    // Throttle control
+    private static float _throttleTarget = 1.0f;          // 0‑1, set by arrow keys
+    private static float _throttleCurrent = 0.0f;         // actual current fraction (lerped)
+    private const float ThrottleLerpRate = 10.0f;          // times per second (speed of movement)
+    private static bool _wKeyHeld = false;
+    private static float _lastThrottleUpdateTime;
+
+    private const int TargetMaxFill = (int)(SampleRate * 0.2);
 
     public static void Main()
     {
-        var mode = new VideoMode(new Vector2u(1280, 720));
+        var window = CreateWindow();
-        var window = new RenderWindow(mode, "Fluid Simulation");
+        LoadFont();
-        window.SetVerticalSyncEnabled(true);
+        _scenario = new SingleCylScenario();
-        window.Closed += (_, _) => { running = false; window.Close(); };
+        _scenario.Initialize(SampleRate);
+        _lastThrottleUpdateTime = 0.0f;
 
-        var soundEngine = new SoundEngine(bufferCapacity: 2048);
+        _simRingBuffer = new SimulationRingBuffer(8192);
-        soundEngine.Volume = 70;
+        _soundEngine = new SoundEngine(_simRingBuffer) { Volume = 100 };
-        soundEngine.Start();
+        _soundEngine.Start();
+
+        var cts = new CancellationTokenSource();
+        Task.Run(() => SimulationLoop(cts.Token), cts.Token);
 
         var stopwatch = Stopwatch.StartNew();
-
+        double lastDrawTime = 0.0;
-        // --- Warmup: fill audio buffer with silence ---
-        int warmupSamples = SampleRate / 2;                // 0.5 s
-        float[] warmup = new float[warmupSamples];
-        for (int i = 0; i < warmupSamples; i++)
-            warmup[i] = 0;
-        soundEngine.WriteSamples(warmup, warmupSamples);
-
-        // Reset timer after warmup – this is the “real‑time zero”
-        stopwatch.Restart();
-        stepCount = 0;                                     // simulation steps start now
-
-        // --- Initialise the simulation scenario ---
-        Simulation.Initialize(SampleRate);
-
-        const int chunkSize = 2048;
-        float[] buffer = new float[chunkSize];
-
-        double lastLogTime = 0.0;                         // for periodic speed printout
 
         while (window.IsOpen)
         {
             window.DispatchEvents();
 
-            // --- Compute how many audio samples are needed since last frame ---
+            double now = stopwatch.Elapsed.TotalSeconds;
-            double currentTime = stopwatch.Elapsed.TotalSeconds;
-            double elapsed = currentTime;                  // since stopwatch was reset
-            int samplesNeeded = (int)(elapsed * SampleRate) - (int)(stepCount);
-            // (stepCount is total generated samples, so we just need the remainder)
 
-            // --- Generate the required number of simulation steps ---
+            // ---- Playback speed smoothing ----
-            while (samplesNeeded > 0 && running)
+            double targetSpeed = _timeWarpActive ? 1.0 : _desiredSpeed;
+            _currentDisplaySpeed += (targetSpeed - _currentDisplaySpeed) *
+                (1.0 - Math.Exp(-8.0 * (now - lastDrawTime)));
+            _soundEngine.Speed = _currentDisplaySpeed;
+
+            // ---- Throttle update ----
+            float dtThrottle = (float)now - _lastThrottleUpdateTime;
+            _lastThrottleUpdateTime = (float)now;
+
+            float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
+
+            // Snap to zero instantly when target is zero (key released)
+            if (throttleDesiredFraction == 0.0)
             {
-                int toGenerate = Math.Min(samplesNeeded, chunkSize);
+                _throttleCurrent = 0.0f;
-                for (int i = 0; i < toGenerate; i++)
-                {
-                    buffer[i] = Simulation.Process();
-                    stepCount++;
             }
-                soundEngine.WriteSamples(buffer, toGenerate);
+            else
-                samplesNeeded -= toGenerate;
+            {
+                float smoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
+                _throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
             }
 
-            // --- Display speed ---
+            _scenario.Throttle = _throttleCurrent;
-            double simTime = stepCount / (double)SampleRate;
-            double wallTime = stopwatch.Elapsed.TotalSeconds;
-            double speed = (wallTime > 0) ? simTime / wallTime : 0.0;
 
-            // Update window title with instant speed
+            // ---- Drawing ----
-            window.SetTitle($"FluidSim | Speed: {speed:F3}× | Steps: {stepCount}");
+            if (now - lastDrawTime >= 1.0 / DrawFrequency)
-
-            // Console log once per second
-            if (wallTime - lastLogTime >= 1.0)
             {
-                Console.WriteLine($"Speed: {speed:F3}×  ({stepCount} steps, {wallTime:F2}s wall)");
+                if (_overlayText != null)
-                lastLogTime = wallTime;
+                {
+                    string toggleHint = _isRealTime ? "[Space] slow mo" : "[Space] real time";
+                    _overlayText.DisplayedString =
+                        $"{toggleHint}  Speed: {_currentDisplaySpeed:F3}x  RT: {(_currentDisplaySpeed * 100.0):F1}%  Sim load: {_loadTracker.LoadPercent:F0}%\n" +
+                        $"Throttle: {_throttleCurrent * 100:F0}%  Target: {_throttleTarget * 100:F0}%  [W] {(_wKeyHeld ? "BLIP" : "---")}";
                 }
 
-            // --- Rendering (placeholder) ---
                 window.Clear(Color.Black);
+                _scenario.Draw(window);
+                if (_overlayText != null) window.Draw(_overlayText);
                 window.Display();
+                lastDrawTime = now;
+            }
         }
 
-        // --- Cleanup ---
+        cts.Cancel();
-        soundEngine.Dispose();
+        _soundEngine.Dispose();
         window.Dispose();
     }
+
+    private static void SimulationLoop(CancellationToken token)
+    {
+        while (!token.IsCancellationRequested)
+        {
+            long cycleStart = Stopwatch.GetTimestamp();
+
+            long workStart = Stopwatch.GetTimestamp();
+            float sample = _scenario.Process();
+            _simRingBuffer.Write(sample);
+            long workEnd = Stopwatch.GetTimestamp();
+
+            while (_simRingBuffer.AvailableSamples > TargetMaxFill &&
+                   !token.IsCancellationRequested)
+            {
+                Thread.Sleep(1);
+            }
+
+            long cycleEnd = Stopwatch.GetTimestamp();
+
+            double busyMs = (workEnd - workStart) / (double)Stopwatch.Frequency * 1000.0;
+            double totalMs = (cycleEnd - cycleStart) / (double)Stopwatch.Frequency * 1000.0;
+            _loadTracker.Record(busyMs, totalMs);
+        }
+    }
+
+    // ---------- Window & input ----------
+    private static RenderWindow CreateWindow()
+    {
+        var mode = new VideoMode(new Vector2u(1280, 720));
+        var window = new RenderWindow(mode, "FluidSim");
+        window.SetVerticalSyncEnabled(false);
+        window.SetFramerateLimit(60);
+        window.Closed += (_, _) => window.Close();
+        window.MouseWheelScrolled += OnMouseWheel;
+        window.KeyPressed += OnKeyPressed;
+        window.KeyReleased += OnKeyReleased;
+        return window;
+    }
+
+    private static void LoadFont()
+    {
+        try { _overlayFont = new Font("fonts/FiraCodeNerdFont-Medium.ttf"); }
+        catch { _overlayFont = null; }
+
+        if (_overlayFont != null)
+            _overlayText = new Text(_overlayFont)
+            {
+                FillColor = Color.White,
+                Position = new Vector2f(10, 10)
+            };
+    }
+
+    private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
+    {
+        if (_timeWarpActive) return;
+        if (e.Delta > 0) _desiredSpeed *= ScrollFactor;
+        else if (e.Delta < 0) _desiredSpeed /= ScrollFactor;
+        _desiredSpeed = Math.Clamp(_desiredSpeed, MinSpeed, MaxSpeed);
+        _lastNormalSpeed = _desiredSpeed;
+        _isRealTime = Math.Abs(_desiredSpeed - 1.0) < 1e-6;
+    }
+
+    private static void OnKeyPressed(object? sender, KeyEventArgs e)
+    {
+        switch (e.Code)
+        {
+            case Keyboard.Key.Space:
+                _timeWarpActive = !_timeWarpActive;
+                if (!_timeWarpActive)
+                {
+                    _desiredSpeed = _lastNormalSpeed;
+                    _isRealTime = false;
+                }
+                break;
+
+            case Keyboard.Key.W:
+                _wKeyHeld = true;
+                break;
+
+            case Keyboard.Key.Up:
+                _throttleTarget = MathF.Min(1.0f, _throttleTarget + 0.05f);
+                break;
+
+            case Keyboard.Key.Down:
+                _throttleTarget = MathF.Max(0.0f, _throttleTarget - 0.05f);
+                break;
+        }
+    }
+
+    private static void OnKeyReleased(object? sender, KeyEventArgs e)
+    {
+        if (e.Code == Keyboard.Key.W)
+            _wKeyHeld = false;
+    }
 }

Scenarios/HelmholtzScenario.cs

@@ -0,0 +1,118 @@
+using FluidSim.Components;
+using FluidSim.Core;
+using FluidSim.Interfaces;
+using SFML.Graphics;
+using SFML.System;
+using System;
+
+namespace FluidSim.Tests
+{
+    public class HelmholtzScenario : Scenario
+    {
+        private Volume0D cavity;
+        private Port cavityPort;
+
+        private PipeSystem pipeSystem;
+        private int[] pipeStart = { 0 };
+        private int[] pipeEnd;
+
+        private BoundarySystem boundaries;
+        private int cavityOrificeIdx = 0;
+        private int openEndIdx = 0;
+
+        private Solver solver;
+        private double dt;
+        private int stepCount;
+
+        private SoundProcessor soundProcessor;
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            // --- Realistic Helmholtz resonator dimensions ---
+            float cavityVolume = 1e-3f;          // 1 liter
+            float neckLength = 0.05f;            // 5 cm
+            float neckDiameter = 0.02f;          // 2 cm diameter
+            float neckArea = MathF.PI * 0.25f * neckDiameter * neckDiameter;
+            int neckCells = 20;
+
+            // --- Volume (cavity) ---
+            float initialPressure = 1.1f * 101325f;   // slight overpressure
+            float initialTemperature = 300f;
+            cavity = new Volume0D(cavityVolume, initialPressure, initialTemperature);
+            cavityPort = cavity.CreatePort();
+
+            // --- Pipe (neck) ---
+            float[] areas = new float[neckCells];
+            float[] dxs = new float[neckCells];
+            float dx = neckLength / neckCells;
+            for (int i = 0; i < neckCells; i++)
+            {
+                areas[i] = neckArea;
+                dxs[i] = dx;
+            }
+            pipeEnd = new[] { neckCells };
+
+            float rho0 = 101325f / (287f * 300f);
+            pipeSystem = new PipeSystem(neckCells, pipeStart, pipeEnd, areas, dxs, rho0, 0f, 101325f);
+
+            // --- Boundary system ---
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
+
+            // Standard orifice with built‑in minor loss (K = 0.5) – no inertance needed
+            boundaries.AddOrificeWithInertance(
+                cavityPort, pipeIndex: 0, isLeftEnd: true,
+                areaIndex: cavityOrificeIdx,
+                dischargeCoeff: 0.9f,
+                effectiveLength: neckLength   // physical neck length
+            );
+
+            // Open end at right side of pipe
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, neckArea);
+
+            float[] orificeAreas = new float[1] { neckArea };
+            boundaries.SetOrificeAreas(orificeAreas);
+
+            // --- Solver ---
+            solver = new Solver { SubStepCount = 8, EnableProfiling = false };
+            solver.SetTimeStep(dt);
+            solver.SetPipeSystem(pipeSystem);
+            solver.SetBoundarySystem(boundaries);
+            solver.AddComponent(cavity);
+
+            // --- Sound ---
+            soundProcessor = new SoundProcessor(sampleRate, 1f) { Gain = 2f };
+
+            Console.WriteLine("Helmholtz resonator ready.");
+            stepCount = 0;
+        }
+
+        public override float Process()
+        {
+            solver.Step();
+            stepCount++;
+
+            float flow = boundaries.GetOpenEndMassFlow(openEndIdx);
+            float sample = soundProcessor.Process(flow);
+
+            return sample;
+        }
+
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+
+            float cavityCenterX = 100f;
+            float cavityWidth = 80f, cavityHeight = 100f;
+            float cavityTopY = winH / 2f - cavityHeight / 2f;
+            DrawVolume(target, cavity, cavityCenterX, cavityTopY - 40f, cavityWidth, cavityHeight);
+
+            float pipeStartX = cavityCenterX + cavityWidth / 2f + 10f;
+            float pipeEndX = winW - 50f;
+            float pipeCenterY = winH / 2f;
+            DrawPipe(target, pipeSystem, 0, pipeCenterY, pipeStartX, pipeEndX);
+        }
+    }
+}

Scenarios/Scenario.cs

@@ -0,0 +1,161 @@
+using SFML.Graphics;
+using SFML.System;
+using FluidSim.Core;
+using FluidSim.Components;
+
+namespace FluidSim.Tests
+{
+    public abstract class Scenario
+    {
+        protected const float AmbientPressure = 101325f;
+        protected const float AmbientTemperature = 300f;
+        public float Throttle { get; set; }
+
+        public abstract void Initialize(int sampleRate);
+        public abstract float Process();
+        public abstract void Draw(RenderWindow target);
+
+        protected Color PressureColor(float pressurePa)
+        {
+            float bar = pressurePa / 1e5f;
+            byte r, g, b;
+            if (bar < 1f)
+            {
+                float f = Math.Clamp(bar, 0f, 1f);
+                r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
+            }
+            else
+            {
+                float f = Math.Min((bar - 1f) / 9f, 1f);
+                r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
+            }
+            return new Color(r, g, b);
+        }
+
+        protected Color TemperatureColor(float t)
+        {
+            t = Math.Clamp(t, 0f, 2000f);
+            byte r, g, b;
+            if (t < AmbientTemperature)
+            {
+                float f = t / AmbientTemperature;
+                r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
+            }
+            else
+            {
+                float f = (t - AmbientTemperature) / (2000f - AmbientTemperature);
+                r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
+            }
+            return new Color(r, g, b);
+        }
+
+        protected void DrawVolume(RenderWindow target, Volume0D volume,
+                                  float centerX, float topY, float width, float height)
+        {
+            var rect = new RectangleShape(new Vector2f(width, height))
+            {
+                FillColor = PressureColor(volume.Pressure),
+                Position = new Vector2f(centerX - width / 2f, topY)
+            };
+            target.Draw(rect);
+            var border = new RectangleShape(new Vector2f(width, height))
+            {
+                FillColor = Color.Transparent,
+                OutlineColor = Color.White,
+                OutlineThickness = 1f,
+                Position = rect.Position
+            };
+            target.Draw(border);
+        }
+
+        protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
+                                    float centerX, float topY, float width, float maxHeight)
+        {
+            float fraction = cylinder.PistonFraction;
+            float currentHeight = maxHeight * fraction;
+            var wall = new RectangleShape(new Vector2f(width, maxHeight))
+            {
+                FillColor = new Color(60, 60, 60),
+                Position = new Vector2f(centerX - width / 2f, topY)
+            };
+            target.Draw(wall);
+            var gas = new RectangleShape(new Vector2f(width, currentHeight))
+            {
+                FillColor = PressureColor(cylinder.Pressure),
+                Position = new Vector2f(centerX - width / 2f, topY)
+            };
+            target.Draw(gas);
+            var piston = new RectangleShape(new Vector2f(width, 4f))
+            {
+                FillColor = Color.White,
+                Position = new Vector2f(centerX - width / 2f, topY + currentHeight)
+            };
+            target.Draw(piston);
+            float valveW = 6f, valveH = 10f, valveY = topY + 4f;
+            var iv = new RectangleShape(new Vector2f(valveW, valveH))
+            {
+                FillColor = cylinder.IntakeValveArea > 0f ? Color.Green : Color.Red,
+                Position = new Vector2f(centerX - width / 2f - valveW - 2f, valveY)
+            };
+            target.Draw(iv);
+            var ev = new RectangleShape(new Vector2f(valveW, valveH))
+            {
+                FillColor = cylinder.ExhaustValveArea > 0f ? Color.Green : Color.Red,
+                Position = new Vector2f(centerX + width / 2f + 2f, valveY)
+            };
+            target.Draw(ev);
+        }
+
+        protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex,
+                                float pipeCenterY, float pipeStartX, float pipeEndX)
+        {
+            int start = pipeSystem.GetPipeStart(pipeIndex);
+            int end = pipeSystem.GetPipeEnd(pipeIndex);
+            int n = end - start;
+            if (n < 2) return;
+
+            float pipeLen = pipeEndX - pipeStartX;
+            float dx = pipeLen / (n - 1);
+            float baseRadius = 25f;
+
+            var centers = new float[n];
+            var radii = new float[n];
+            var temps = new float[n];
+            for (int i = 0; i < n; i++)
+            {
+                int cell = start + i;
+                float p = pipeSystem.GetCellPressure(cell);
+                float rho = pipeSystem.GetCellDensity(cell);
+                temps[i] = p / MathF.Max(rho * 287f, 1e-12f);
+                float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
+                radii[i] = baseRadius * (1f + dev * 2f);
+                if (radii[i] < 2f) radii[i] = 2f;
+                centers[i] = pipeStartX + i * dx;
+            }
+
+            int segments = 8;
+            var va = new VertexArray(PrimitiveType.TriangleStrip);
+            for (int i = 0; i < n; i++)
+            {
+                float x = centers[i], r = radii[i];
+                Color col = TemperatureColor(temps[i]);
+                va.Append(new Vertex(new Vector2f(x, pipeCenterY - r), col));
+                va.Append(new Vertex(new Vector2f(x, pipeCenterY + r), col));
+                if (i < n - 1)
+                {
+                    for (int s = 1; s <= segments; s++)
+                    {
+                        float t = s / (float)segments;
+                        float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
+                        float ri = radii[i] + (radii[i + 1] - radii[i]) * t;
+                        float Ti = temps[i] + (temps[i + 1] - temps[i]) * t;
+                        Color colS = TemperatureColor(Ti);
+                        va.Append(new Vertex(new Vector2f(xi, pipeCenterY - ri), colS));
+                        va.Append(new Vertex(new Vector2f(xi, pipeCenterY + ri), colS));
+                    }
+                }
+            }
+            target.Draw(va);
+        }
+    }
+}

Scenarios/Simulation.cs

@@ -1,93 +0,0 @@
-using FluidSim.Audio;
-using FluidSim.Components;
-using FluidSim.Core;
-using FluidSim.Utilities;
-
-namespace FluidSim.Scenarios
-{
-    public static class Simulation
-    {
-        private static Solver solver;
-        private static Volume0D cavity, ambient;
-        private static Pipe1D neck;
-        private static Connection connNeckCavity, connNeckAmbient;
-        private static int stepCount;
-        private static double time;
-        private static double dt;
-
-        public static void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            // --- Cavity (the “bottle”) ---
-            double V_cav = 1.0 * Units.L;                 // 1 litre
-            cavity = new Volume0D(V_cav, 10 * Units.atm, Units.Celsius(20), sampleRate);
-
-            // --- Ambient (a huge “room”) ---
-            double V_amb = 1000.0;                        // 1000 m³ ≈ constant pressure
-            ambient = new Volume0D(V_amb, 1.0 * Units.atm, Units.Celsius(20), sampleRate);
-
-            // --- Neck (short pipe) ---
-            double length_neck = 80.0 * Units.mm;
-            double diam_neck = 10 * Units.mm;
-            double area_neck = Units.AreaFromRadius(diam_neck, Units.mm);                             // few cells enough for a short neck
-            neck = new Pipe1D(length_neck, area_neck, sampleRate);
-            neck.SetUniformState(ambient.Density, 0.0, ambient.Pressure);
-            neck.FrictionFactor = 0.02;                    // slight damping
-
-            // --- Connections ---
-            // Cavity-to-neck (full area)
-            connNeckCavity = new Connection(cavity.Port, neck.PortA)
-            {
-                Area = area_neck,
-                DischargeCoefficient = 1.0,
-                Gamma = 1.4
-            };
-            // Neck-to-ambient (SoundConnection to capture the radiated tone)
-            connNeckAmbient = new SoundConnection(neck.PortB, ambient.Port)
-            {
-                Area = area_neck,
-                DischargeCoefficient = 1.0,
-                Gamma = 1.4
-            };
-
-            // --- Solver ---
-            solver = new Solver();
-            solver.AddVolume(cavity);
-            solver.AddVolume(ambient);
-            solver.AddPipe(neck);
-            solver.AddConnection(connNeckCavity);
-            solver.AddConnection(connNeckAmbient);
-
-            // --- Sound tuning ---
-            SoundProcessor.MaxDeltaP = 0.1f * (float)Units.atm; // small Δp expected
-            SoundProcessor.MaxArea = (float)area_neck;
-            SoundProcessor.MaxVelocity = 343f;
-            SoundProcessor.ReferenceDensity = 1.2f;
-            SoundProcessor.ReferenceSpeedOfSound = 343f;
-            SoundProcessor.Gain = 10.0f;  // amplify because Δp is small
-        }
-
-        public static float Process()
-        {
-            solver.Step();
-            time += dt;
-            stepCount++;
-            Log();
-            return solver.LastSample;
-        }
-
-        public static void Log()
-        {
-            if (stepCount <= 200 || stepCount % (int)(0.5 / dt) == 0)
-            {
-                Console.WriteLine(
-                    $"t = {time * 1e3:F3} ms  " +
-                    $"Sample = {solver.LastSample:F6},  " +
-                    $"P_cav = {cavity.Pressure / 1e5:F6} bar,  " +
-                    $"flow_cav = {cavity.Port.MassFlowRate / 1e3:F4} g/s,  " +
-                    $"flow_neck = {neck.PortB.MassFlowRate * 1e3:F4} g/s");
-            }
-        }
-    }
-}

Scenarios/SingleCylScenario.cs

@@ -0,0 +1,246 @@
+using FluidSim.Components;
+using FluidSim.Core;
+using FluidSim.Interfaces;
+using FluidSim.Utils;
+using SFML.Graphics;
+using SFML.System;
+using System;
+
+namespace FluidSim.Tests
+{
+    public class SingleCylScenario : Scenario
+    {
+        private Crankshaft crankshaft;
+        private Cylinder cylinder;
+
+        private PipeSystem pipeSystem;
+        private BoundarySystem boundaries;
+        private Solver solver;
+
+        private Volume0D intakePlenum;
+        private Port plenumInlet, plenumOutlet;
+        private Volume0D exhaustCollector;
+        private Port colIn, colOut;
+
+        private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx;
+        private float[] orificeAreas;
+        private int intakeOpenIdx, exhaustOpenIdx;
+
+        private SoundProcessor exhaustSound, intakeSound;
+        private OutdoorExhaustReverb reverb;
+
+        private double dt;
+        private int stepCount;
+
+        // Use a private field for the maximum throttle area, avoiding any base‑class conflicts
+        private float _maxThrottleArea;
+
+        // pipe area for open end calculations
+        private float pipeArea;
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            // Maximum throttle area – independent of base class
+            _maxThrottleArea = (float)Units.AreaFromDiameter(3 * Units.cm);   // 1 cm²
+
+            // ---- Crankshaft ----
+            crankshaft = new Crankshaft(2000);
+            crankshaft.Inertia = 0.01f;
+            crankshaft.FrictionConstant = 2f;
+            crankshaft.FrictionViscous = 0.0f;
+
+            // ---- Cylinder ----
+            float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 11f;
+            float ivo = 350f, ivc = 580f, evo = 120f, evc = 370f;
+            cylinder = new Cylinder(bore, stroke, conRod, compRatio,
+                                    ivo, ivc, evo, evc, crankshaft)
+            {
+                IntakeValveDiameter = 0.03f,
+                IntakeValveLift = 0.005f,
+                ExhaustValveDiameter = 0.028f,
+                ExhaustValveLift = 0.005f
+            };
+
+            // ---- Pipe system ----
+            int[] pipeStart = { 0, 10, 20 };
+            int[] pipeEnd   = { 10, 20, 70 };
+            int totalCells = pipeEnd[^1];   // automatically 70, stays in sync
+            float[] area = new float[totalCells];
+            float[] dx = new float[totalCells];
+            float pipeDiameter = 0.02f;          // 2 cm
+            pipeArea = MathF.PI * 0.25f * pipeDiameter * pipeDiameter;
+            float areaVal = pipeArea;
+            float intakeLenBefore = 0.2f, intakeLenRunner = 0.2f, exhaustLen = 0.4f;
+            for (int i = 0; i < totalCells; i++)
+            {
+                area[i] = areaVal;
+                if (i < 10) dx[i] = intakeLenBefore / 10f;
+                else if (i < 20) dx[i] = intakeLenRunner / 10f;
+                else dx[i] = exhaustLen / 50f;
+            }
+
+            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
+                                        1.225f, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 1.0f;
+            pipeSystem.EnergyRelaxationRate = 0.5f;
+            pipeSystem.AmbientPressure = 101325f;
+
+            // ---- Volumes ----
+            intakePlenum = new Volume0D(100e-6f, 101325f, 300f);  // 100 mL
+            plenumInlet = intakePlenum.CreatePort();
+            plenumOutlet = intakePlenum.CreatePort();
+            exhaustCollector = new Volume0D(10e-6f, 101325f, 800f);  // 10 mL (unused but present)
+            colIn = exhaustCollector.CreatePort();
+            colOut = exhaustCollector.CreatePort();
+
+            // ---- Boundary system ----
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
+
+            throttleAreaIdx = 0;
+            plenumRunnerAreaIdx = 1;
+            intakeValveIdx = 2;
+            exhaustValveIdx = 3;
+
+            // Intake open end (pipe0 left)
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, pipeArea);
+            intakeOpenIdx = 0;
+
+            // Throttle orifice (plenum inlet to pipe0 right)
+            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.2f);
+
+            // Plenum to runner (plenum outlet to pipe1 left)
+            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1f);
+
+            // Intake valve (cylinder intake to pipe1 right)
+            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1f);
+
+            // Exhaust valve (cylinder exhaust to pipe2 left)
+            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1f);
+
+            // Exhaust open end (pipe2 right)
+            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, pipeArea);
+            exhaustOpenIdx = 1;
+
+            orificeAreas = new float[4];
+            orificeAreas[plenumRunnerAreaIdx] = areaVal; // fixed plenum->runner area
+
+            // ---- Solver ----
+            solver = new Solver { SubStepCount = 4, EnableProfiling = false };
+            solver.SetTimeStep(dt);
+            solver.SetPipeSystem(pipeSystem);
+            solver.SetBoundarySystem(boundaries);
+            solver.AddComponent(cylinder);
+            solver.AddComponent(intakePlenum);
+            solver.AddComponent(exhaustCollector);
+
+            // ---- Sound ----
+            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
+            intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
+            reverb = new OutdoorExhaustReverb(sampleRate);
+
+            stepCount = 0;
+            Console.WriteLine("TestScenario ready.");
+        }
+
+        public override float Process()
+        {
+            crankshaft.Step((float)dt);
+            cylinder.PreStep((float)dt);
+
+            // Update variable orifice areas – use the private _maxThrottleArea
+            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1f);
+            orificeAreas[throttleAreaIdx] = throttledArea;
+            orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
+            orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
+            boundaries.SetOrificeAreas(orificeAreas);
+
+            solver.Step();
+            stepCount++;
+
+            // Retrieve open‑end mass flows for sound synthesis
+            float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
+            float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
+
+            float exhaustDry = exhaustSound.Process(exhaustFlow);
+            float intakeDry = intakeSound.Process(intakeFlow);
+
+            if (stepCount % 1000 == 0)
+            {
+                float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
+                float crankDeg = crankshaft.CrankAngle;   // degrees (0–720)
+                Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1} deg, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar");
+                Console.WriteLine($"  intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
+
+                // Pipe 0 (intake before throttle)
+                var (r0L, u0L, p0L) = pipeSystem.GetInteriorStateLeft(0);
+                var (r0R, u0R, p0R) = pipeSystem.GetInteriorStateRight(0);
+                Console.WriteLine($"  Pipe0 L: rho={r0L:F4} u={u0L:F3} p={p0L/1e5:F3}bar | R: rho={r0R:F4} u={u0R:F3} p={p0R/1e5:F3}bar");
+
+                // Pipe 1 (runner)
+                var (r1L, u1L, p1L) = pipeSystem.GetInteriorStateLeft(1);
+                var (r1R, u1R, p1R) = pipeSystem.GetInteriorStateRight(1);
+                Console.WriteLine($"  Pipe1 L: rho={r1L:F4} u={u1L:F3} p={p1L/1e5:F3}bar | R: rho={r1R:F4} u={u1R:F3} p={p1R/1e5:F3}bar");
+
+                // Pipe 2 (exhaust)
+                var (r2L, u2L, p2L) = pipeSystem.GetInteriorStateLeft(2);
+                var (r2R, u2R, p2R) = pipeSystem.GetInteriorStateRight(2);
+                Console.WriteLine($"  Pipe2 L: rho={r2L:F4} u={u2L:F3} p={p2L/1e5:F3}bar | R: rho={r2R:F4} u={u2R:F3} p={p2R/1e5:F3}bar");
+
+                // Plenum and cylinder mass
+                Console.WriteLine($"  Plenum P={intakePlenum.Pressure/1e5:F3}bar, mass={intakePlenum.Mass:E4} kg");
+                Console.WriteLine($"  Cyl mass={cylinder.Mass:E4} kg");
+            }
+
+            return reverb.Process(intakeDry + exhaustDry);
+        }
+
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+
+            float intakeY = winH / 2f - 40f;
+            float exhaustY = winH / 2f + 80f;
+            float openEndX = 40f;
+
+            // Intake pipe before throttle (pipe 0)
+            float pipe1StartX = openEndX;
+            float pipe1EndX = pipe1StartX + 120f;
+            DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
+
+            // Throttle symbol
+            float throttleX = pipe1EndX + 5f;
+            var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
+            {
+                FillColor = Color.Yellow,
+                Position = new Vector2f(throttleX, intakeY - 15f)
+            };
+            target.Draw(throttleRect);
+
+            // Plenum
+            float plenW = 60f, plenH = 80f;
+            float plenLeftX = throttleX + 10f;
+            float plenCenterX = plenLeftX + plenW / 2f;
+            float plenTopY = intakeY - plenH / 2f;
+            DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
+
+            // Runner pipe (pipe 1)
+            float runnerStartX = plenLeftX + plenW + 5f;
+            float runnerEndX = runnerStartX + 100f;
+            DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
+
+            // Cylinder
+            float cylCX = runnerEndX + 50f;
+            float cylTopY = intakeY - 120f;
+            float cylW = 80f, cylMaxH = 240f;
+            DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
+
+            // Exhaust pipe (pipe 2)
+            float exhStartX = cylCX + cylW / 2f + 20f;
+            float exhEndX = winW - 60f;
+            DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
+        }
+    }
+}

Scenarios/TestScenario.cs

@@ -0,0 +1,102 @@
+using System;
+using SFML.Graphics;
+using SFML.System;
+using FluidSim.Core;
+
+namespace FluidSim.Tests
+{
+    public class TestScenario : Scenario
+    {
+        private PipeSystem pipeSystem;
+        private BoundarySystem boundaries;
+        private Solver solver;
+
+        private int[] pipeStart = { 0 };
+        private int[] pipeEnd;
+
+        private double dt;
+        private int stepCount;
+
+        // Sound output: use pressure at open end
+        private SoundProcessor openEndSound;
+        private int openEndIdx = 0;   // index of the open end in BoundarySystem (we added only one)
+
+        public override void Initialize(int sampleRate)
+        {
+            dt = 1.0 / sampleRate;
+
+            const int cellCount = 200;
+            float length = 2f;
+            float dia = 0.02f;
+            float area = MathF.PI * 0.25f * dia * dia;
+
+            float[] areas = new float[cellCount];
+            float[] dxs = new float[cellCount];
+            float dx = length / cellCount;
+            for (int i = 0; i < cellCount; i++)
+            {
+                areas[i] = area;
+                dxs[i] = dx;
+            }
+
+            pipeEnd = new[] { cellCount };
+
+            float rho0 = 101325f / (287f * 300f);
+            pipeSystem = new PipeSystem(cellCount, pipeStart, pipeEnd, areas, dxs,
+                                        rho0, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 0f;
+            pipeSystem.EnergyRelaxationRate = 0f;
+            pipeSystem.AmbientPressure = 101325f;
+
+            // Pressure bubble near right end
+            float pBubble = 10f * 101325f;
+            float TBubble = 2000f;
+            float rhoBubble = pBubble / (287f * TBubble);
+            for (int i = 0; i <= 10; i++)
+                pipeSystem.SetCellState(i, rhoBubble, 0f, pBubble);
+
+            // Boundaries: left closed, right open
+            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
+            boundaries.AddOrifice(null, pipeIndex: 0, isLeftEnd: true, areaIndex: 0, 1f);
+            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, area);
+            float[] orificeAreas = new float[1] { 0f };
+            boundaries.SetOrificeAreas(orificeAreas);
+
+            solver = new Solver { SubStepCount = 3};
+            solver.SetTimeStep(dt);
+            solver.SetPipeSystem(pipeSystem);
+            solver.SetBoundarySystem(boundaries);
+
+            solver.EnableProfiling = true;
+            pipeSystem.EnableProfiling = true;
+
+            // Simple sound processor: convert mass flow rate to audio
+            openEndSound = new SoundProcessor(sampleRate, 1f) { Gain = 2f };
+
+            Console.WriteLine("Pulse test ready.");
+            stepCount = 0;
+        }
+
+        public override float Process()
+        {
+            solver.Step();
+            stepCount++;
+
+            float flow = boundaries.GetOpenEndMassFlow(openEndIdx);
+            float sample = openEndSound.Process(flow);
+
+            return sample;
+        }
+
+        public override void Draw(RenderWindow target)
+        {
+            float winW = target.GetView().Size.X;
+            float winH = target.GetView().Size.Y;
+
+            float startX = 50f;
+            float endX = winW - 50f;
+            float y = winH / 2f;
+            DrawPipe(target, pipeSystem, 0, y, startX, endX);
+        }
+    }
+}

Utils/Units.cs

@@ -1,4 +1,6 @@
-namespace FluidSim.Utilities
+using System;
+
+namespace FluidSim.Utils
 {
     public static class Units
     {
@@ -19,10 +21,10 @@
 
         public static double Celsius(double tC) => tC + 273.15;
 
-        public static double AreaFromRadius(double radius, double unit = mm) =>
+        public static double AreaFromRadius(double radius) =>
-            Math.PI * (radius * unit) * (radius * unit);
+            Math.PI * (radius) * (radius);
 
-        public static double AreaFromDiameter(double diameter, double unit = mm) =>
+        public static double AreaFromDiameter(double diameter) =>
-            Math.PI * 0.25 * (diameter * unit) * (diameter * unit);
+            Math.PI * 0.25 * (diameter) * (diameter);
     }
 }