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

@@ -0,0 +1,53 @@
+using SFML.Audio;
+using SFML.System;
+
+namespace FluidSim.Audio
+{
+    public class AudioOutputStream : SoundStream
+    {
+        private readonly SimulationRingBuffer _sourceBuffer;
+        private double _speed = 1.0;   // non‑volatile, accessed with Volatile.Read/Write
+
+        public AudioOutputStream(SimulationRingBuffer sourceBuffer)
+        {
+            _sourceBuffer = sourceBuffer;
+            // 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;
+            float[] temp = new float[monoBlockSize];
+
+            int read = _sourceBuffer.ReadInterpolated(temp, monoBlockSize, Speed);
+            samples = new short[monoBlockSize * 2];
+
+            if (read > 0)
+            {
+                for (int i = 0; i < read; i++)
+                {
+                    float clamped = Math.Clamp(temp[i], -1f, 1f);
+                    short final = (short)(clamped * short.MaxValue);
+                    samples[i * 2] = final;      // left
+                    samples[i * 2 + 1] = final;  // right
+                }
+            }
+            // Fill rest with silence
+            for (int i = read * 2; i < samples.Length; i++)
+                samples[i] = 0;
+
+            return true;
+        }
+
+        protected override void OnSeek(Time timeOffset) =>
+            throw new NotSupportedException();
+    }
+}

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

@@ -0,0 +1,45 @@
+namespace FluidSim.Audio
+{
+    public class SoundEngine : IDisposable
+    {
+        private readonly AudioOutputStream _stream;
+        private bool _isPlaying;
+
+        public SoundEngine(SimulationRingBuffer sourceBuffer, int bufferCapacity = 16384)
+        {
+            _stream = new AudioOutputStream(sourceBuffer);
+        }
+
+        public void Start()
+        {
+            if (_isPlaying) return;
+            _stream.Play();
+            _isPlaying = true;
+        }
+
+        public void Stop()
+        {
+            if (!_isPlaying) return;
+            _stream.Stop();
+            _isPlaying = false;
+        }
+
+        public double Speed
+        {
+            get => _stream.Speed;
+            set => _stream.Speed = value;
+        }
+
+        public float Volume
+        {
+            get => _stream.Volume;
+            set => _stream.Volume = value;
+        }
+
+        public void Dispose()
+        {
+            Stop();
+            _stream.Dispose();
+        }
+    }
+}

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/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,420 +0,0 @@
-using System;
-using FluidSim.Interfaces;
-
-namespace FluidSim.Components
-{
-    public enum BoundaryType
-    {
-        VolumeCoupling,
-        OpenEnd,
-        ClosedEnd
-    }
-
-    public class Pipe1D
-    {
-        public Port PortA { get; }
-        public Port PortB { get; }
-        public double Area => _area;
-        public double DampingMultiplier { get; set; } = 1.0;
-
-        private int _n;
-        private double _dx, _dt, _gamma, _area, _diameter;
-        private double[] _rho, _rhou, _E;
-
-        // Volume‑coupling ghost states for boundaries A and B
-        private double _rhoA, _pA;
-        private double _rhoB, _pB;
-        private bool _aBCSet, _bBCSet;
-
-        private BoundaryType _aBCType = BoundaryType.VolumeCoupling;
-        private BoundaryType _bBCType = BoundaryType.VolumeCoupling;
-
-        private double _aAmbientPressure = 101325.0;
-        private double _bAmbientPressure = 101325.0;
-
-        private const double CflTarget = 0.8;
-        private const double ReferenceSoundSpeed = 340.0;
-
-        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);
-
-        public BoundaryType ABCType => _aBCType;
-        public BoundaryType BBCType => _bBCType;
-
-        public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
-        {
-            double dtGlobal = 1.0 / sampleRate;
-            int nCells;
-
-            if (forcedCellCount > 1)
-            {
-                nCells = forcedCellCount;
-            }
-            else
-            {
-                double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
-                nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
-                while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
-                    nCells++;
-            }
-
-            _n = nCells;
-            _dx = length / _n;
-            _dt = dtGlobal;
-            _area = area;
-            _gamma = 1.4;
-
-            // Hydraulic diameter for a circular pipe
-            _diameter = 2.0 * Math.Sqrt(area / Math.PI);
-
-            _rho = new double[_n];
-            _rhou = new double[_n];
-            _E = new double[_n];
-
-            PortA = new Port();
-            PortB = new Port();
-        }
-
-        public void SetABoundaryType(BoundaryType type) => _aBCType = type;
-        public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
-        public void SetAAmbientPressure(double p) => _aAmbientPressure = p;
-        public void SetBAmbientPressure(double p) => _bAmbientPressure = p;
-
-        public void 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 void SetCellState(int i, double rho, double u, double p)
-        {
-            if (i < 0 || i >= _n) return;
-            _rho[i] = rho;
-            _rhou[i] = rho * u;
-            double e = p / ((_gamma - 1) * rho);
-            _E[i] = rho * e + 0.5 * rho * u * u;
-        }
-
-        public void SetAVolumeState(double rhoVol, double pVol)
-        {
-            _rhoA = rhoVol;
-            _pA = pVol;
-            _aBCSet = true;
-        }
-
-        public void SetBVolumeState(double rhoVol, double pVol)
-        {
-            _rhoB = rhoVol;
-            _pB = pVol;
-            _bBCSet = true;
-        }
-
-        public void ClearBC() => _aBCSet = _bBCSet = false;
-
-        public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
-        {
-            double maxW = 0.0;
-            for (int i = 0; i < _n; i++)
-            {
-                double rho = _rho[i];
-                double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
-                double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
-                double local = u + c;
-                if (local > maxW) maxW = local;
-            }
-            maxW = Math.Max(maxW, 1e-8);
-            return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
-        }
-
-        public void SimulateSingleStep(double dtSub)
-        {
-            int n = _n;
-            double[] Fm = new double[n + 1];
-            double[] Fp = new double[n + 1];
-            double[] Fe = new double[n + 1];
-
-            // ---------- Boundary A (face 0, left) ----------
-            double rhoIntA = _rho[0];
-            double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
-            double pIntA = Pressure(0);
-
-            switch (_aBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_aBCSet)
-                    {
-                        HLLCFlux(_rhoA, 0.0, _pA,
-                                 rhoIntA, uIntA, pIntA,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    }
-                    else
-                    {
-                        Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false,
-                                  out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-            }
-
-            // ---------- Internal faces ----------
-            for (int i = 0; i < n - 1; i++)
-            {
-                double rhoL = _rho[i];
-                double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
-                double pL = Pressure(i);
-
-                double rhoR = _rho[i + 1];
-                double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
-                double pR = Pressure(i + 1);
-
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR,
-                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
-            }
-
-            // ---------- Boundary B (face n, right) ----------
-            double rhoIntB = _rho[n - 1];
-            double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12);
-            double pIntB = Pressure(n - 1);
-
-            switch (_bBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_bBCSet)
-                    {
-                        HLLCFlux(rhoIntB, uIntB, pIntB,
-                                 _rhoB, 0.0, _pB,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    }
-                    else
-                    {
-                        Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true,
-                                  out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-            }
-
-            // ---- Cell update with linear laminar damping ----
-            double radius = _diameter / 2.0;
-            double mu_air = 1.8e-5;
-            double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
-
-            for (int i = 0; i < n; i++)
-            {
-                double dM = (Fm[i + 1] - Fm[i]) / _dx;
-                double dP = (Fp[i + 1] - Fp[i]) / _dx;
-                double dE = (Fe[i + 1] - Fe[i]) / _dx;
-
-                _rho[i] -= dtSub * dM;
-                _rhou[i] -= dtSub * dP;
-                _E[i] -= dtSub * dE;
-
-                double rho = Math.Max(_rho[i], 1e-12);
-                double dampingRate = laminarCoeff / rho;
-                double dampingFactor = Math.Exp(-dampingRate * dtSub);
-                _rhou[i] *= dampingFactor;
-
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
-                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
-                double pMin = 100.0;
-                double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
-                if (_E[i] < eMin) _E[i] = eMin;
-            }
-
-            // ---------- Port quantities ----------
-            double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
-            double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
-
-            PortA.MassFlowRate = mdotA_sub;
-            PortB.MassFlowRate = mdotB_sub;
-            PortA.Pressure = pIntA;
-            PortB.Pressure = pIntB;
-            PortA.Density = _rho[0];
-            PortB.Density = _rho[n - 1];
-
-            // Corrected enthalpy for both directions
-            if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet)
-            {
-                PortA.SpecificEnthalpy = mdotA_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(0)
-                    : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12);
-            }
-            if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet)
-            {
-                PortB.SpecificEnthalpy = mdotB_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(_n - 1)
-                    : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12);
-            }
-        }
-
-        private double GetCellTotalSpecificEnthalpy(int i)
-        {
-            double rho = Math.Max(_rho[i], 1e-12);
-            double u = _rhou[i] / rho;
-            double p = Pressure(i);
-            double h = _gamma / (_gamma - 1.0) * p / rho;
-            return h + 0.5 * u * u;
-        }
-
-        private double Pressure(int i) =>
-            (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
-
-        // ========== Characteristic‑based Open End ==========
-        private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
-            if (uInt <= -cInt)   // supersonic inflow – use interior state as ghost
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt <= 0)   // subsonic inflow
-            {
-                // Reservoir condition: p = pAmb, T = 300K, u = 0
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-                return;
-            }
-            else   // subsonic outflow (uInt > 0)
-            {
-                // Ghost pressure forced to pAmb
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary)
-                double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
-
-                // Prevent spurious inflow by clipping to zero
-                if (uGhost < 0) uGhost = 0;
-
-                HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-            }
-        }
-
-        private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            if (uInt >= cInt)   // supersonic outflow (extrapolation)
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt >= 0)   // subsonic outflow
-            {
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary)
-                double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
-
-                // Clip to zero to prevent inflow
-                if (uGhost > 0) uGhost = 0;
-
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe);
-            }
-            else   // subsonic inflow
-            {
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe);
-            }
-        }
-
-        // ========== Closed end (mirror) ==========
-        private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary,
-                                   out double fm, out double fp, out double fe)
-        {
-            double rhoGhost = rhoInt;
-            double pGhost = pInt;
-            double uGhost = -uInt;   // mirror velocity
-
-            if (isRightBoundary)
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
-            else
-                HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
-        }
-
-        // ========== Standard HLLC flux ==========
-        private void HLLCFlux(double rL, double uL, double pL,
-                              double rR, double uR, double pR,
-                              out double fm, out double fp, out double fe)
-        {
-            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12));
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12));
-            double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
-            double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR;
-            double SL = Math.Min(uL - cL, uR - cR);
-            double SR = Math.Max(uL + cL, uR + cR);
-
-            double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
-                      / (rL * (SL - uL) - rR * (SR - uR));
-
-            double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
-            double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
-
-            if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
-            else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
-            else if (Ss >= 0)
-            {
-                double rsL = rL * (SL - uL) / (SL - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
-                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
-                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
-            }
-            else
-            {
-                double rsR = rR * (SR - uR) / (SR - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
-                double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
-                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
-            }
-        }
-
-        public double GetPressureAtFraction(double fraction)
-        {
-            int i = (int)(fraction * (_n - 1));
-            i = Math.Clamp(i, 0, _n - 1);
-            return Pressure(i);
-        }
-    }
-}

Components/Volume0D.cs

@@ -1,84 +1,133 @@
 using System;
+using System.Collections.Generic;
 using FluidSim.Interfaces;
-using FluidSim.Utils;
 
 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; }
-        public double InternalEnergy { get; private set; }
+        private float _airMass;
+        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;
-        public double GasConstant { get; set; } = 287.0;
+        // ---------- Thermal relaxation to environment ----------
+        /// <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 double dVdt { get; set; }
+        public float Mass => _airMass + _exhaustMass;
+        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 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)
+        public Volume0D(float initialVolume, float initialPressure,
+                        float initialTemperature, float gasConstant = 287f, float gamma = 1.4f)
         {
             GasConstant = gasConstant;
             Gamma = gamma;
             Volume = initialVolume;
-            dVdt = 0.0;
-            _dt = 1.0 / sampleRate;
+            Dvdt = 0f;
 
-            double rho0 = initialPressure / (GasConstant * initialTemperature);
-            Mass = rho0 * Volume;
-            InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
-
-            Port = new Port();
-            PushStateToPort();
+            float rho0 = initialPressure / (GasConstant * initialTemperature);
+            _airMass = rho0 * Volume;
+            _exhaustMass = 0f;
+            InternalEnergy = (initialPressure * Volume) / (Gamma - 1f);
         }
 
-        public void PushStateToPort()
+        public Port CreatePort()
         {
-            Port.Pressure = Pressure;
-            Port.Density = Density;
-            Port.Temperature = Temperature;
-            Port.SpecificEnthalpy = SpecificEnthalpy;
+            var port = new Port { Owner = this };
+            port.Pressure = Pressure;
+            port.Density = Density;
+            port.Temperature = Temperature;
+            port.SpecificEnthalpy = SpecificEnthalpy;
+            port.AirFraction = AirFraction;
+            Ports.Add(port);
+            return port;
         }
 
-        // Original integrate (uses the constructor’s sample rate)
-        public void Integrate()
+        public void SetPressure(float pressure, float? temperature = null)
         {
-            Integrate(_dt);
+            float V = MathF.Max(Volume, 1e-12f);
+            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);
         }
 
-        public void SetPressure(double newPressure)
+        public void UpdateState(float dt)
         {
-            InternalEnergy = newPressure * Volume / (Gamma - 1.0);
-            // Mass stays the same, so density is unchanged
+            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;
             }
 
-        // New overload: integrate with a custom time step (for sub‑steps)
-        public void Integrate(double dtOverride)
-        {
-            double mdot = Port.MassFlowRate;
-            double h_in = Port.SpecificEnthalpy;
+            float dAir = totalMdotAir * dt;
+            float dExhaust = totalMdotExhaust * dt;
+            float dE = totalEdot * dt - Pressure * Dvdt * dt;
 
-            double dm = mdot * dtOverride;
-            double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride;
-
-            Mass += dm;
+            _airMass += dAir;
+            _exhaustMass += dExhaust;
             InternalEnergy += dE;
 
-            // Hard physical bounds – prevent NaN and unphysical states
-            if (Mass < 1e-12) Mass = 1e-12;
-            if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
-
-            PushStateToPort();
+            // ---- Thermal relaxation ----
+            if (EnergyRelaxationRate > 0f)
+            {
+                float currentMass = Mass;
+                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,100 +0,0 @@
-using System;
-using FluidSim.Interfaces;
-
-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,5 +1,7 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
 using FluidSim.Components;
 using FluidSim.Interfaces;
 
@@ -7,164 +9,94 @@ namespace FluidSim.Core
 {
     public class Solver
     {
-        private readonly List<Volume0D> _volumes = new();
-        private readonly List<Pipe1D> _pipes = new();
-        private readonly List<Connection> _connections = new();
-
+        private readonly List<IComponent> _components = new();
+        private PipeSystem _pipeSystem;
+        private BoundarySystem _boundarySystem;
         private double _dt;
 
-        public void AddVolume(Volume0D v) => _volumes.Add(v);
-        public void AddPipe(Pipe1D p) => _pipes.Add(p);
-        public void AddConnection(Connection c) => _connections.Add(c);
+        public int SubStepCount { get; set; } = 4;
+        public bool EnableProfiling { get; set; } = false;
+
+        private long _stepCount;
+        private long _ticksOrifice, _ticksOpenEnd, _ticksPipe, _ticksUpdate;
+
         public void SetTimeStep(double dt) => _dt = dt;
+        public void AddComponent(IComponent component) => _components.Add(component);
 
-        /// <summary>
-        /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right).
-        /// </summary>
-        public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0)
+        public void SetPipeSystem(PipeSystem pipeSystem)
         {
-            if (isA)
-            {
-                pipe.SetABoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetAAmbientPressure(ambientPressure);
+            _pipeSystem = pipeSystem;
         }
-            else
+        public void SetBoundarySystem(BoundarySystem boundarySystem)
         {
-                pipe.SetBBoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetBAmbientPressure(ambientPressure);
-            }
+            _boundarySystem = boundarySystem;
         }
 
-        public float Step()
+        public void Step()
         {
-            // 1. Volumes publish state
-            foreach (var v in _volumes)
-                v.PushStateToPort();
+            if (_pipeSystem == null || _boundarySystem == null) return;
 
-            // 2. Set volume BCs for volume‑coupled ends
-            foreach (var conn in _connections)
-            {
-                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                {
-                    var pipe = GetPipe(conn.PortA);
-                    bool isA = pipe.PortA == conn.PortA;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortA, conn.PortB);
-                }
-                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                {
-                    var pipe = GetPipe(conn.PortB);
-                    bool isA = pipe.PortB == conn.PortB;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortB, conn.PortA);
-                }
-            }
-
-            // 3. Sub‑steps
-            int nSub = 1;
-            foreach (var p in _pipes)
-                nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
-            double dtSub = _dt / nSub;
+            int nSub = _pipeSystem.GetRequiredSubSteps((float)_dt, 0.8f);
+            nSub = Math.Max(nSub, SubStepCount);   // never go below fixed minimum
+            float dtSub = (float)(_dt / nSub);
 
             for (int sub = 0; sub < nSub; sub++)
             {
-                foreach (var p in _pipes)
-                    p.SimulateSingleStep(dtSub);
+                long t0;
 
-                foreach (var conn in _connections)
-                {
-                    if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortA);
-                        bool isA = pipe.PortA == conn.PortA;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortA, conn.PortB, dtSub);
-                    }
-                    else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortB);
-                        bool isA = pipe.PortB == conn.PortB;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
-                    }
+                t0 = Stopwatch.GetTimestamp();
+                _boundarySystem.ResolveOrifices(dtSub);
+                _ticksOrifice += Stopwatch.GetTimestamp() - t0;
+
+                t0 = Stopwatch.GetTimestamp();
+                _boundarySystem.ResolveOpenEnds(dtSub);
+                _ticksOpenEnd += Stopwatch.GetTimestamp() - t0;
+
+                t0 = Stopwatch.GetTimestamp();
+                _pipeSystem.SimulateStep(dtSub);
+                _ticksPipe += Stopwatch.GetTimestamp() - t0;
             }
 
-                if (sub < nSub - 1)
+            long tUS = Stopwatch.GetTimestamp();
+            foreach (var comp in _components)
+                comp.UpdateState((float)_dt);
+            _ticksUpdate += Stopwatch.GetTimestamp() - tUS;
+
+            _stepCount++;
+            if (_stepCount % 5000 == 0 && EnableProfiling)
             {
-                    foreach (var v in _volumes)
-                        v.PushStateToPort();
+                double freq = Stopwatch.Frequency;
+                double total = _ticksOrifice + _ticksOpenEnd + _ticksPipe + _ticksUpdate;
+                double avgStepUs = (total / freq) * 1e6 / 5000.0;
 
-                    foreach (var conn in _connections)
+                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)
                 {
-                        if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortA);
-                            bool isA = pipe.PortA == conn.PortA;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortA, conn.PortB);
-                        }
-                        else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortB);
-                            bool isA = pipe.PortB == conn.PortB;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortB, conn.PortA);
-                        }
-                    }
-                }
+                    Console.WriteLine(_pipeSystem.GetProfileReport());
                 }
 
-            // 5. Audio samples (none for now, but placeholder)
-            var audioSamples = new List<float>();
-            foreach (var conn in _connections)
-            {
-                if (conn is SoundConnection sc)
-                    audioSamples.Add(sc.GetAudioSample());
+                _ticksOrifice = _ticksOpenEnd = _ticksPipe = _ticksUpdate = 0;
             }
-
-            // 6. Clear BC flags
-            foreach (var p in _pipes)
-                p.ClearBC();
-
-            return SoundProcessor.MixAndClip(audioSamples.ToArray());
-        }
-
-        private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
-        private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
-        private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
-        private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p);
-
-        private void SetVolumeBC(Port pipePort, Port volPort)
-        {
-            var pipe = GetPipe(pipePort);
-            if (pipe == null) return;
-            bool isA = pipe.PortA == pipePort;
-            if (isA)
-                pipe.SetAVolumeState(volPort.Density, volPort.Pressure);
-            else
-                pipe.SetBVolumeState(volPort.Density, volPort.Pressure);
-        }
-
-        private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
-        {
-            double mdot = pipePort.MassFlowRate;
-            volPort.MassFlowRate = -mdot;
-
-            if (mdot < 0)  // pipe → volume
-            {
-                volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
-            }
-            // else volume’s own enthalpy (from PushStateToPort) is used
-
-            GetVolume(volPort)?.Integrate(dtSub);
         }
     }
 }

Core/SoundEngine.cs

@@ -1,131 +0,0 @@
-using SFML.Audio;
-using SFML.System;
-
-namespace FluidSim;
-
-#region Lock‑free ring buffer (unchanged)
-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;
-    }
-}
-#endregion
-
-#region Stereo stream that consumes the ring buffer
-internal class RingBufferStream : SoundStream
-{
-    private readonly RingBuffer ringBuffer;
-
-    public RingBufferStream(RingBuffer buffer)
-    {
-        ringBuffer = buffer;
-        // 2 channels, 44.1 kHz, standard stereo mapping
-        Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
-    }
-
-    protected override bool OnGetData(out short[] samples)
-    {
-        const int monoBlockSize = 512;          // number of mono samples we'll read
-        float[] temp = new float[monoBlockSize];
-        int read = ringBuffer.Read(temp, monoBlockSize);
-        samples = new short[monoBlockSize * 2];
-
-        if (read > 0)
-        {
-            for (int i = 0; i < read; i++)
-            {
-                float clamped = Math.Clamp(temp[i], -1f, 1f);
-                short final = (short)(clamped * short.MaxValue);
-                samples[i * 2] = final;   // left
-                samples[i * 2 + 1] = final;   // right
-            }
-        }
-        for (int i = read * 2; i < samples.Length; i++)
-            samples[i] = 0;
-
-        return true;
-    }
-
-    protected override void OnSeek(Time timeOffset) =>
-        throw new NotSupportedException();
-}
-#endregion
-
-#region Public sound engine API (unchanged)
-public class SoundEngine : IDisposable
-{
-    private readonly RingBuffer ringBuffer;
-    private readonly RingBufferStream stream;
-    private bool isPlaying;
-
-    public SoundEngine(int bufferCapacity = 16384)
-    {
-        ringBuffer = new RingBuffer(bufferCapacity);
-        stream = new RingBufferStream(ringBuffer);
-    }
-
-    public void Start()
-    {
-        if (isPlaying) return;
-        stream.Play();
-        isPlaying = true;
-    }
-
-    public void Stop()
-    {
-        if (!isPlaying) return;
-        stream.Stop();
-        isPlaying = false;
-        float[] drain = new float[ringBuffer.Count];
-        ringBuffer.Read(drain, drain.Length);
-    }
-
-    public int WriteSamples(float[] data, int count) =>
-        ringBuffer.Write(data, count);
-
-    public float Volume
-    {
-        get => stream.Volume;
-        set => stream.Volume = value;
-    }
-
-    public void Dispose()
-    {
-        Stop();
-        stream.Dispose();
-    }
-}
-#endregion

Core/SoundProcessor.cs

@@ -1,23 +1,34 @@
-namespace FluidSim.Core
-{
-    /// <summary>
-    /// Mixes multiple audio samples and applies a soft‑clipping tanh.
-    /// </summary>
-    public static class SoundProcessor
-    {
-        /// <summary>Overall gain applied after mixing (before tanh).</summary>
-        public static float MasterGain { get; set; } = 0.01f;
+using System;
 
-        /// <summary>
-        /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1].
-        /// </summary>
-        public static float MixAndClip(params float[] samples)
+namespace FluidSim.Core
+{
+    public class SoundProcessor
     {
-            float sum = 0f;
-            foreach (float s in samples)
-                sum += s;
-            sum *= MasterGain;
-            return sum;
+        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/Connection.cs

@@ -1,15 +0,0 @@
-namespace FluidSim.Interfaces
-{
-    /// <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);
-    }
-}

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

Interfaces/SoundConnection.cs

@@ -1,25 +0,0 @@
-namespace FluidSim.Interfaces
-{
-    /// <summary>
-    /// A Connection that also produces an audio sample from the pressure drop across it.
-    /// </summary>
-    public class SoundConnection : Connection
-    {
-        /// <summary>Gain applied to the normalised pressure difference.</summary>
-        public float Gain { get; set; } = 1.0f;
-
-        /// <summary>Reference pressure used for normalisation (Pa). Default: 1 atm.</summary>
-        public double ReferencePressure { get; set; } = 101325.0;
-
-        public SoundConnection(Port a, Port b) : base(a, b) { }
-
-        /// <summary>
-        /// Returns a normalised audio sample proportional to the pressure difference.
-        /// </summary>
-        public float GetAudioSample()
-        {
-            double dp = PortA.Pressure - PortB.Pressure;
-            return (float)(dp / ReferencePressure) * Gain;
-        }
-    }
-}

Program.cs

@@ -1,8 +1,13 @@
-using SFML.Graphics;
-using SFML.Window;
-using SFML.System;
-using System.Diagnostics;
+using FluidSim.Audio;
 using FluidSim.Core;
+using FluidSim.Tests;
+using SFML.Graphics;
+using SFML.System;
+using SFML.Window;
+using System;
+using System.Diagnostics;
+using System.Threading;
+using System.Threading.Tasks;
 
 namespace FluidSim;
 
@@ -10,175 +15,202 @@ public class Program
 {
     private const int SampleRate = 44100;
     private const double DrawFrequency = 60.0;
-    private static Scenario scenario;
 
-    // Speed control
-    //private static double desiredSpeed = 1.0;
-    private static double desiredSpeed = 0.0001;
-    private static double currentSpeed = desiredSpeed;
+    // 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;
 
-    // Space‑toggle state
-    private static double lastDesiredSpeed = 0.1;   // remembers the last non‑1.0 scroll speed
-    private static bool isRealTime = true;          // true when desiredSpeed == 1.0
+    private static volatile bool _timeWarpActive;
 
-    private static volatile bool running = true;
+    // 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 = new RenderWindow(mode, "Pipe Resonator");
-        window.SetVerticalSyncEnabled(true);
-        window.Closed += (_, _) => { running = false; window.Close(); };
-        window.MouseWheelScrolled += OnMouseWheel;
-        window.KeyPressed += OnKeyPressed;
+        var window = CreateWindow();
+        LoadFont();
+        _scenario = new SingleCylScenario();
+        _scenario.Initialize(SampleRate);
+        _lastThrottleUpdateTime = 0.0f;
 
-        var soundEngine = new SoundEngine(bufferCapacity: 16384);
-        soundEngine.Volume = 70;
-        soundEngine.Start();
+        _simRingBuffer = new SimulationRingBuffer(8192);
+        _soundEngine = new SoundEngine(_simRingBuffer) { Volume = 100 };
+        _soundEngine.Start();
 
-        //scenario = new PipeResonatorScenario();
-        //scenario = new HelmholtzResonatorScenario();
-        scenario = new SodShockTubeScenario();
-
-        scenario.Initialize(SampleRate);
+        var cts = new CancellationTokenSource();
+        Task.Run(() => SimulationLoop(cts.Token), cts.Token);
 
         var stopwatch = Stopwatch.StartNew();
         double lastDrawTime = 0.0;
-        double drawInterval = 1.0 / DrawFrequency;
-        double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
-
-        // Resampling buffer
-        List<float> simBuffer = new List<float>(4096);
-        double readIndex = 0.0;
-
-        for (int i = 0; i < 4; i++)
-            simBuffer.Add(scenario.Process());
-
-        long totalSimSteps = simBuffer.Count;
-        long totalOutputSamples = 0;
-
-        double lastRealTime = stopwatch.Elapsed.TotalSeconds;
-        const int outputChunk = 256;
-        float[] outputBuf = new float[outputChunk];
 
         while (window.IsOpen)
         {
             window.DispatchEvents();
 
-            double currentRealTime = stopwatch.Elapsed.TotalSeconds;
-            double dtSpeed = currentRealTime - lastSpeedUpdateTime;
-            lastSpeedUpdateTime = currentRealTime;
+            double now = stopwatch.Elapsed.TotalSeconds;
 
-            // Smoothly transition currentSpeed → desiredSpeed
-            // When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp
-            double smoothingRate = 8.0;   // higher = faster catch‑up
-            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
+            // ---- Playback speed smoothing ----
+            double targetSpeed = _timeWarpActive ? 1.0 : _desiredSpeed;
+            _currentDisplaySpeed += (targetSpeed - _currentDisplaySpeed) *
+                (1.0 - Math.Exp(-8.0 * (now - lastDrawTime)));
+            _soundEngine.Speed = _currentDisplaySpeed;
 
-            // ---------- Generate audio ----------
-            double targetAudioClock = currentRealTime + 0.05;
+            // ---- Throttle update ----
+            float dtThrottle = (float)now - _lastThrottleUpdateTime;
+            _lastThrottleUpdateTime = (float)now;
 
-            while (totalOutputSamples < targetAudioClock * SampleRate && running)
+            float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
+
+            // Snap to zero instantly when target is zero (key released)
+            if (throttleDesiredFraction == 0.0)
             {
-                int toGenerate = (int)Math.Min(
-                    (long)outputChunk,
-                    (long)(targetAudioClock * SampleRate) - totalOutputSamples
-                );
-                if (toGenerate <= 0) break;
-
-                double maxIndex = readIndex + (toGenerate - 1) * currentSpeed + 2;
-                int requiredSimIndex = (int)Math.Ceiling(maxIndex);
-                while (simBuffer.Count - 1 < requiredSimIndex)
+                _throttleCurrent = 0.0f;
+            }
+            else
             {
-                    simBuffer.Add(scenario.Process());
-                    totalSimSteps++;
+                float smoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
+                _throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
             }
 
-                for (int i = 0; i < toGenerate; i++)
+            _scenario.Throttle = _throttleCurrent;
+
+            // ---- Drawing ----
+            if (now - lastDrawTime >= 1.0 / DrawFrequency)
             {
-                    int i0 = (int)readIndex;
-                    int i1 = i0 + 1;
-                    double frac = readIndex - i0;
-
-                    float y0 = simBuffer[Math.Clamp(i0, 0, simBuffer.Count - 1)];
-                    float y1 = simBuffer[Math.Clamp(i1, 0, simBuffer.Count - 1)];
-                    outputBuf[i] = (float)(y0 + (y1 - y0) * frac);
-
-                    readIndex += currentSpeed;
-
-                    while (readIndex >= 1.0 && simBuffer.Count > 2)
+                if (_overlayText != null)
                 {
-                        simBuffer.RemoveAt(0);
-                        readIndex -= 1.0;
+                    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" : "---")}";
                 }
-                }
-
-                int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
-                totalOutputSamples += accepted;
-
-                if (accepted < toGenerate)
-                    break;
-            }
-
-            // ---------- Drawing & title ----------
-            if (currentRealTime - lastDrawTime >= drawInterval)
-            {
-                double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate);
-                double simTime = totalSimSteps / (double)SampleRate;
-                string toggleHint = isRealTime ? "[Space] slow mo" : "[Space] real time";
-                window.SetTitle(
-                    $"{toggleHint}  Sim: {simTime:F2}s | " +
-                    $"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
-                    $"Actual: {actualSpeed:F2}x"
-                );
 
                 window.Clear(Color.Black);
-                scenario.Draw(window);
+                _scenario.Draw(window);
+                if (_overlayText != null) window.Draw(_overlayText);
                 window.Display();
-                lastDrawTime = currentRealTime;
+                lastDrawTime = now;
             }
         }
 
-        soundEngine.Dispose();
+        cts.Cancel();
+        _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)
     {
-        bool wasRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-
-        if (e.Delta > 0)
-            desiredSpeed *= ScrollFactor;
-        else if (e.Delta < 0)
-            desiredSpeed /= ScrollFactor;
-
-        desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed);
-
-        // Update the remembered slow-mo speed (unless we are exactly at 1.0)
-        if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6)
-            lastDesiredSpeed = desiredSpeed;
-
-        // Update isRealTime flag
-        isRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
+        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)
     {
-        if (e.Code == Keyboard.Key.Space)
+        switch (e.Code)
         {
-            if (isRealTime)
+            case Keyboard.Key.Space:
+                _timeWarpActive = !_timeWarpActive;
+                if (!_timeWarpActive)
                 {
-                // Switch to the remembered slow speed
-                desiredSpeed = lastDesiredSpeed;
+                    _desiredSpeed = _lastNormalSpeed;
+                    _isRealTime = false;
                 }
-            else
+                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)
     {
-                // Switch back to real time
-                desiredSpeed = 1.0;
-            }
-            isRealTime = !isRealTime;
-        }
+        if (e.Code == Keyboard.Key.W)
+            _wKeyHeld = false;
     }
 }

Scenarios/HelmholtzResonatorScenario.cs

@@ -1,133 +0,0 @@
-using System;
-using FluidSim.Components;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-
-namespace FluidSim.Core
-{
-    public class HelmholtzResonatorScenario : Scenario
-    {
-        private Solver solver;
-        private Volume0D cavity;
-        private Pipe1D neck;
-        private Connection coupling;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            // 1‑litre cavity, 10% over‑pressure
-            double cavityVolume = 1e-3;
-            double initialCavityPressure = 1.1 * ambientPressure;
-            cavity = new Volume0D(cavityVolume, initialCavityPressure, 300.0, sampleRate)
-            {
-                Gamma = 1.4,
-                GasConstant = 287.0
-            };
-
-            // Neck: length 10 cm, radius 1 cm
-            double neckLength = 0.1;
-            double neckRadius = 0.01;
-            double neckArea = Math.PI * neckRadius * neckRadius;
-            neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40);
-            neck.SetUniformState(1.225, 0.0, ambientPressure);
-
-            coupling = new Connection(neck.PortA, cavity.Port)
-            {
-                Area = neckArea,
-                DischargeCoefficient = 0.62,
-                Gamma = 1.4
-            };
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddVolume(cavity);
-            solver.AddPipe(neck);
-            solver.AddConnection(coupling);
-
-            // Port A (left) = volume coupling, Port B (right) = open end
-            solver.SetPipeBoundary(neck, isA: true, BoundaryType.VolumeCoupling);
-            solver.SetPipeBoundary(neck, isA: false, BoundaryType.OpenEnd, ambientPressure);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pOpen = neck.GetCellPressure(neck.GetCellCount() - 1);
-            float audio = (float)((pOpen - ambientPressure) / ambientPressure);
-
-            if (stepCount % 20 == 0)
-            {
-                double pCav = cavity.Pressure;
-                double mdotA = neck.PortA.MassFlowRate;   // positive = into pipe (leaving cavity)
-                Console.WriteLine(
-                    $"t={time * 1e3:F2} ms step={stepCount} " +
-                    $"P_cav={pCav:F1} Pa, P_open={pOpen:F1} Pa, " +
-                    $"mdot_A={mdotA * 1e3:F4} g/s, audio={audio:F4}");
-            }
-
-            return audio;
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winW = target.GetView().Size.X;
-            float winH = target.GetView().Size.Y;
-            float centerY = winH / 2f;
-
-            // Cavity rectangle
-            float cavityWidth = 120f;
-            float cavityHeight = 180f;
-            var cavityRect = new RectangleShape(new Vector2f(cavityWidth, cavityHeight));
-            cavityRect.Position = new Vector2f(40f, centerY - cavityHeight / 2f);
-            cavityRect.FillColor = PressureColor(cavity.Pressure);
-            target.Draw(cavityRect);
-
-            // Neck drawn as tapered pipe
-            int n = neck.GetCellCount();
-            float neckStartX = 40f + cavityWidth + 10f;
-            float neckEndX = winW - 60f;
-            float neckLenPx = neckEndX - neckStartX;
-            float dx = neckLenPx / (n - 1);
-            float baseRadius = 20f;
-
-            Vertex[] vertices = new Vertex[n * 2];
-            for (int i = 0; i < n; i++)
-            {
-                float x = neckStartX + i * dx;
-                double p = neck.GetCellPressure(i);
-                float r = baseRadius * (float)(0.5 + 0.5 * Math.Tanh((p - ambientPressure) / (ambientPressure * 0.2)));
-                if (r < 4f) r = 4f;
-                Color col = PressureColor(p);
-                vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
-                vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
-            }
-            target.Draw(vertices, PrimitiveType.TriangleStrip);
-
-            // Open end indicator
-            var arrow = new CircleShape(8f);
-            arrow.Position = new Vector2f(neckEndX - 4f, centerY - 4f);
-            arrow.FillColor = Color.White;
-            target.Draw(arrow);
-        }
-
-        private Color PressureColor(double pressure)
-        {
-            double range = ambientPressure * 0.1;
-            double t = Math.Clamp((pressure - ambientPressure) / range, -1.0, 1.0);
-            byte r = (byte)(t > 0 ? 255 * t : 0);
-            byte b = (byte)(t < 0 ? -255 * t : 0);
-            byte g = (byte)(255 * (1 - Math.Abs(t)));
-            return new Color(r, g, b);
-        }
-    }
-}

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/PipeResonatorScenario.cs

@@ -1,184 +0,0 @@
-using FluidSim.Components;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-using System;
-
-namespace FluidSim.Core
-{
-    public class PipeResonatorScenario : Scenario
-    {
-        private Solver solver;
-        private Pipe1D pipe;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-        private bool enableLogging = true;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            double length = 2;
-            double radius = 50 * Units.mm;
-            double area = Units.AreaFromDiameter(radius);
-
-            pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: 80);
-            pipe.SetUniformState(1.225, 0.0, ambientPressure);
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddPipe(pipe);
-            // Open end at port A (left), closed end at port B (right)
-            solver.SetPipeBoundary(pipe, isA: true, BoundaryType.OpenEnd, ambientPressure);
-            solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd);
-
-            // Initial pressure pulse
-            int pulseCells = 5;
-            double pulsePressure = 2 * ambientPressure;
-            for (int i = 0; i < pulseCells; i++)
-                pipe.SetCellState(i, 1.225, 0.0, pulsePressure);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pMid = pipe.GetPressureAtFraction(0.5);
-            sample = (float)((pMid - ambientPressure) / ambientPressure);
-
-            Log(sample);
-            return sample;
-        }
-
-        private void Log(float sample)
-        {
-            if (!enableLogging) return;
-            if (stepCount % 10 == 0 && stepCount < 1000)
-            {
-                double pMid = pipe.GetPressureAtFraction(0.5);
-                double pOpen = pipe.GetCellPressure(0);
-                double pClosed = pipe.GetCellPressure(pipe.GetCellCount() - 1);
-                Console.WriteLine(
-                    $"t = {time * 1e3:F3} ms  Step {stepCount:D4}: " +
-                    $"sample = {sample:F3}, " +
-                    $"P_mid = {pMid:F2} Pa ({pMid / ambientPressure:F4} atm), " +
-                    $"P_open = {pOpen:F2} Pa, P_closed = {pClosed:F2} Pa");
-            }
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winWidth = target.GetView().Size.X;
-            float winHeight = target.GetView().Size.Y;
-
-            float pipeCenterY = winHeight / 2f;
-            float margin = 60f;
-            float pipeStartX = margin;
-            float pipeEndX = winWidth - margin;
-            float pipeLengthPx = pipeEndX - pipeStartX;
-            int n = pipe.GetCellCount();
-            float dx = pipeLengthPx / (n - 1);   // spacing between cell centres
-
-            float baseRadius = 25f;
-            float rangeFactor = 1f;
-            float scaleFactor = 5f;
-
-            // ----- smoothstep helper -----
-            static float SmoothStep(float edge0, float edge1, float x)
-            {
-                float t = Math.Clamp((x - edge0) / (edge1 - edge0), 0f, 1f);
-                return t * t * (3f - 2f * t);
-            }
-
-            // ----- Pre‑compute cell positions and radii -----
-            var centers = new float[n];
-            var radii = new float[n];
-            for (int i = 0; i < n; i++)
-            {
-                double p = pipe.GetCellPressure(i);
-                float deviation = (float)Math.Tanh((p - ambientPressure) / ambientPressure / rangeFactor);
-                radii[i] = baseRadius * (1f + deviation * scaleFactor);
-                if (radii[i] < 2f) radii[i] = 2f;
-                centers[i] = pipeStartX + i * dx;
-            }
-
-            // ----- Build triangle‑strip vertices -----
-            int segmentsPerCell = 8;   // smoothness
-            int totalPoints = n + (n - 1) * segmentsPerCell;
-            Vertex[] stripVertices = new Vertex[totalPoints * 2];   // top + bottom for each point
-            int idx = 0;
-
-            for (int i = 0; i < n; i++)
-            {
-                // ---- Cell centre ----
-                float x = centers[i];
-                float r = radii[i];
-                double p = pipe.GetCellPressure(i);
-                Color col = PressureColor(p);
-
-                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col);
-                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col);
-
-                // ---- Intermediate segments after this cell (if not last) ----
-                if (i < n - 1)
-                {
-                    for (int s = 1; s <= segmentsPerCell; s++)
-                    {
-                        float t = s / (float)segmentsPerCell;
-                        float st = SmoothStep(0f, 1f, t);
-                        float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
-                        float ri = radii[i] + (radii[i + 1] - radii[i]) * st;
-                        double pi = pipe.GetCellPressure(i) * (1 - t) + pipe.GetCellPressure(i + 1) * t;
-                        Color coli = PressureColor(pi);
-
-                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli);
-                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), coli);
-                    }
-                }
-            }
-
-            // Draw the pipe as a triangle strip
-            var pipeMesh = new VertexArray(PrimitiveType.TriangleStrip, (uint)stripVertices.Length);
-            for (int i = 0; i < stripVertices.Length; i++)
-                pipeMesh[(uint)i] = stripVertices[i];
-            target.Draw(pipeMesh);
-
-            // ----- Closed end indicator (right) -----
-            float wallThickness = 8f;
-            var wall = new RectangleShape(new Vector2f(wallThickness, winHeight * 0.6f));
-            wall.Position = new Vector2f(pipeEndX, pipeCenterY - winHeight * 0.6f / 2f);
-            wall.FillColor = new Color(180, 180, 180);
-            target.Draw(wall);
-        }
-
-        /// <summary>Blue (low) → Green (ambient) → Red (high).</summary>
-        private Color PressureColor(double pressure)
-        {
-            double range = ambientPressure * 0.05;  // ±5% gives full colour swing
-            double t = (pressure - ambientPressure) / range;
-            t = Math.Clamp(t, -1.0, 1.0);
-
-            byte r, g, b;
-            if (t < 0)
-            {
-                double factor = -t;
-                r = 0;
-                g = (byte)(255 * (1 - factor));
-                b = (byte)(255 * factor);
-            }
-            else
-            {
-                double factor = t;
-                r = (byte)(255 * factor);
-                g = (byte)(255 * (1 - factor));
-                b = 0;
-            }
-            return new Color(r, g, b);
-        }
-    }
-}

Scenarios/Scenario.cs

@@ -1,23 +1,161 @@
 using SFML.Graphics;
+using SFML.System;
+using FluidSim.Core;
+using FluidSim.Components;
 
-namespace FluidSim.Core
+namespace FluidSim.Tests
 {
     public abstract class Scenario
     {
-        /// <summary>
-        /// Initialize the scenario with a given audio sample rate.
-        /// </summary>
+        protected const float AmbientPressure = 101325f;
+        protected const float AmbientTemperature = 300f;
+        public float Throttle { get; set; }
+
         public abstract void Initialize(int sampleRate);
-
-        /// <summary>
-        /// Advance one simulation step and return an audio sample.
-        /// The step size is 1 / sampleRate seconds.
-        /// </summary>
         public abstract float Process();
-
-        /// <summary>
-        /// Draw the current simulation state onto the given SFML render target.
-        /// </summary>
         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/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/SodShockTubeScenario.cs

@@ -1,158 +0,0 @@
-using System;
-using FluidSim.Components;
-using FluidSim.Utils;
-using SFML.Graphics;
-using SFML.System;
-
-namespace FluidSim.Core
-{
-    public class SodShockTubeScenario : Scenario
-    {
-        private Solver solver;
-        private Pipe1D pipe;
-        private int stepCount;
-        private double time;
-        private double dt;
-        private double ambientPressure = 1.0 * Units.atm;
-        private const double GasConstant = 287.0;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-            double length = 1.0;
-            double area = 1.0;
-            int nCells = 200;
-
-            pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: nCells);
-            pipe.SetUniformState(0.125, 0.0, 0.1 * ambientPressure);  // right state
-
-            // Left half high pressure
-            for (int i = 0; i < nCells / 2; i++)
-                pipe.SetCellState(i, 1.0, 0.0, ambientPressure);
-
-            solver = new Solver();
-            solver.SetTimeStep(dt);
-            solver.AddPipe(pipe);
-            solver.SetPipeBoundary(pipe, isA: true, BoundaryType.ClosedEnd);
-            solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd);
-        }
-
-        public override float Process()
-        {
-            float sample = solver.Step();
-            time += dt;
-            stepCount++;
-
-            double pMid = pipe.GetPressureAtFraction(0.5);
-            float audio = (float)((pMid - ambientPressure) / ambientPressure);
-
-            bool log = true;
-
-            if (log)
-            {
-                int n = pipe.GetCellCount();
-                Console.WriteLine($"step {stepCount}:");
-                Console.WriteLine("i    rho (kg/m³)   p (Pa)       T (K)     u (m/s)");
-                for (int i = 0; i < n; i++)
-                {
-                    if (i % 10 == 0)
-                    {
-                        double rho = pipe.GetCellDensity(i);
-                        double p = pipe.GetCellPressure(i);
-                        double u = pipe.GetCellVelocity(i);
-                        double T = p / (rho * GasConstant);   // GasConstant = 287.0
-                        Console.WriteLine($"{i,-4} {rho,10:F4}   {p,10:F1}   {T,8:F2}   {u,10:F4}");
-                    }
-                }
-                Console.WriteLine();
-            }
-
-            return audio;
-        }
-
-        public override void Draw(RenderWindow target)
-        {
-            float winW = target.GetView().Size.X;
-            float winH = target.GetView().Size.Y;
-            float centerY = winH / 2f;
-            float margin = 40f;
-            float pipeStartX = margin;
-            float pipeEndX = winW - margin;
-            float pipeLenPx = pipeEndX - pipeStartX;
-            int n = pipe.GetCellCount();
-            float dx = pipeLenPx / (n - 1);
-            float baseRadius = 60f;
-
-            Vertex[] vertices = new Vertex[n * 2];
-            for (int i = 0; i < n; i++)
-            {
-                float x = pipeStartX + i * dx;
-
-                double p = pipe.GetCellPressure(i);
-                double rho = pipe.GetCellDensity(i);
-                double T = p / (rho * GasConstant);           // temperature in Kelvin
-
-                // Radius from pressure (exaggerated deviation)
-                float r = baseRadius * (float)(p / ambientPressure * 2);
-                if (r < 4f) r = 4f;
-
-                // Colour from temperature
-                Color col = TemperatureColor(T);
-
-                vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);
-                vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col);
-            }
-            target.Draw(vertices, PrimitiveType.TriangleStrip);
-
-            // Diaphragm marker (faint white line at initial interface)
-            float diaphragmX = pipeStartX + (n / 2) * dx;
-            var line = new RectangleShape(new Vector2f(2f, winH * 0.5f));
-            line.Position = new Vector2f(diaphragmX - 1f, centerY - winH * 0.25f);
-            line.FillColor = new Color(255, 255, 255, 80);
-            target.Draw(line);
-        }
-
-        /// <summary>
-        /// Custom temperature‑to‑hue mapping that matches the given Sod‑tube hue values:
-        ///   250 K → 176, 300 K → 122, 350 K → 120?, 450 K → 71.
-        /// Interpolates piecewise linearly, clamping outside [250,450].
-        /// </summary>
-        private Color TemperatureColor(double T)
-        {
-            // 1. Map temperature to hue (0–360)
-            double[] Tknots = { 250, 282, 353, 450 };
-            double[] Hknots = { 176, 179, 122, 71 };
-            double hue;
-            if (T <= Tknots[0]) hue = Hknots[0];
-            else if (T >= Tknots[^1]) hue = Hknots[^1];
-            else
-            {
-                int i = 0;
-                while (i < Tknots.Length - 1 && T > Tknots[i + 1]) i++;
-                double frac = (T - Tknots[i]) / (Tknots[i + 1] - Tknots[i]);
-                hue = Hknots[i] + frac * (Hknots[i + 1] - Hknots[i]);
-            }
-
-            // 2. Convert hue to RGB (S = 1, V = 1)
-            double h = hue / 60.0;
-            int sector = (int)Math.Floor(h);
-            double f = h - sector;
-            byte p = 0;
-            byte q = (byte)(255 * (1 - f));
-            byte tByte = (byte)(255 * f);
-            byte v = 255;
-
-            byte r, g, b;
-            switch (sector % 6)
-            {
-                case 0: r = v; g = tByte; b = p; break;
-                case 1: r = q; g = v; b = p; break;
-                case 2: r = p; g = v; b = tByte; break;
-                case 3: r = p; g = q; b = v; break;
-                case 4: r = tByte; g = p; b = v; break;
-                default: r = v; g = p; b = q; break;
-            }
-            return new Color(r, g, b);
-        }
-    }
-}

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

Sources/EffortSource.cs

@@ -1,10 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Text;
-
-namespace FluidSim.Sources
-{
-    internal class EffortSource
-    {
-    }
-}

Sources/FlowSource.cs

@@ -1,10 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Text;
-
-namespace FluidSim.Sources
-{
-    internal class FlowSource
-    {
-    }
-}