Attempting to generate sound, and set cells automatically

.vscode/launch.json

@@ -1,26 +0,0 @@
-{
-    "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

@@ -1,41 +0,0 @@
-{
-    "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/RingBuffer.cs

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

Audio/RingBufferStream.cs

@@ -3,31 +3,22 @@ using SFML.System;
 
 namespace FluidSim.Audio
 {
-    public class AudioOutputStream : SoundStream
+    internal class RingBufferStream : SoundStream
     {
-        private readonly SimulationRingBuffer _sourceBuffer;
+        private readonly RingBuffer ringBuffer;
-        private double _speed = 1.0;   // non‑volatile, accessed with Volatile.Read/Write
 
-        public AudioOutputStream(SimulationRingBuffer sourceBuffer)
+        public RingBufferStream(RingBuffer buffer)
         {
-            _sourceBuffer = sourceBuffer;
+            ringBuffer = buffer;
-            // 2 channels, 44.1 kHz, stereo
+            // 2 channels, 44.1 kHz, standard stereo mapping
             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;
+            const int monoBlockSize = 512;          // number of mono samples we'll read
             float[] temp = new float[monoBlockSize];
-
+            int read = ringBuffer.Read(temp, monoBlockSize);
-            int read = _sourceBuffer.ReadInterpolated(temp, monoBlockSize, Speed);
             samples = new short[monoBlockSize * 2];
 
             if (read > 0)
@@ -36,11 +27,10 @@ namespace FluidSim.Audio
                 {
                     float clamped = Math.Clamp(temp[i], -1f, 1f);
                     short final = (short)(clamped * short.MaxValue);
-                    samples[i * 2] = final;      // left
+                    samples[i * 2] = final;   // left
-                    samples[i * 2 + 1] = final;  // right
+                    samples[i * 2 + 1] = final;   // right
                 }
             }
-            // Fill rest with silence
             for (int i = read * 2; i < samples.Length; i++)
                 samples[i] = 0;
 

Audio/SimulationRingBuffer.cs

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

Audio/SoundEngine.cs

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

Audio/SoundProcessor.cs

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

Components/Atmosphere.cs

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

Components/Connection.cs

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

Components/Crankshaft.cs

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

Components/Cylinder.cs

@@ -1,117 +0,0 @@
-using System;
-using FluidSim.Components;   // if needed
-
-namespace FluidSim.Components
-{
-    public class Cylinder : EngineCylinder
-    {
-        public float IVO, IVC, EVO, EVC;   // degrees in a 720° cycle
-
-        protected override float CycleLengthRad => 4f * MathF.PI;
-        protected override float MaxCycleDeg => 720f;
-
-        public override float IntakeValveArea =>
-            MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
-        public override float ExhaustValveArea =>
-            MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
-
-        public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio,
-                        float ivo, float ivc, float evo, float evc, Crankshaft crankshaft)
-            : base(bore, stroke, conRodLength, compressionRatio, crankshaft)
-        {
-            IVO = ivo; IVC = ivc; EVO = evo; EVC = evc;
-        }
-
-        private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
-        {
-            float deg = thetaDeg % 720f;
-            if (deg < 0f) deg += 720f;
-
-            float effectiveOpen = opens;
-            float effectiveClose = closes;
-            if (closes < opens) effectiveClose += 720f;
-            float duration = effectiveClose - effectiveOpen;
-            if (duration <= 0f) return 0f;
-
-            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;
-        }
-
-        protected override void HandleCycleEvents(float prevDeg, float currDeg, float dt)
-        {
-            // Intake closing → fuel injection
-            if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
-            {
-                trappedAirMass = _airMass;
-                fuelMass = trappedAirMass / StoichiometricAFR;
-                fuelInjected = true;
-            }
-
-            // Spark – occurs at TDC (0°) minus advance, every 720°
-            float sparkAngle = (0f - SparkAdvance + 720f) % 720f;
-            bool crossedSpark = false;
-            if (prevDeg < sparkAngle && currDeg >= sparkAngle)
-                crossedSpark = true;
-            else if (prevDeg > sparkAngle && currDeg < sparkAngle)
-                crossedSpark = true;
-
-            if (crossedSpark && !combustionActive && fuelInjected)
-            {
-                if (_random.NextDouble() < MisfireProbability)
-                {
-                    combustionActive = false;
-                }
-                else
-                {
-                    combustionActive = true; burnFraction = 0f;
-                    float range = EnergyVariationFraction;
-                    _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
-                }
-            }
-
-            // Combustion progression
-            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;
-                }
-            }
-        }
-    }
-}

Components/EngineCylinder.cs

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

Components/Pipe1D.cs

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

Components/SoundConnection.cs

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

Components/TwoStrokeCylinder.cs

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

Components/Vehicle.cs

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

Components/Volume0D.cs

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

@@ -1,404 +0,0 @@
-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

@@ -1,11 +0,0 @@
-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

@@ -1,27 +0,0 @@
-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

@@ -1,31 +0,0 @@
-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

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

Core/OutdoorExhaustReverb.cs

@@ -1,187 +0,0 @@
-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

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

Core/SoundProcessor.cs

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

Core/ThreadLoadTracker.cs

@@ -1,34 +0,0 @@
-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,21 +5,12 @@
     <TargetFramework>net10.0</TargetFramework>
     <ImplicitUsings>enable</ImplicitUsings>
     <Nullable>enable</Nullable>
-    <PublishAot>false</PublishAot>
+    <PublishAot>true</PublishAot>
     <InvariantGlobalization>true</InvariantGlobalization>
   </PropertyGroup>
 
   <ItemGroup>
-      <PackageReference Include="SFML.Net" Version="3.0.0" />
+    <PackageReference Include="SFML.Net" Version="3.0.0" />
-  </ItemGroup>
-
-  <ItemGroup>
-      <None Update="fonts\FiraCodeNerdFont-Medium.ttf">
-        <CopyToOutputDirectory>Always</CopyToOutputDirectory>
-      </None>
-      <None Update="fonts\LiberationMono-Regular.ttf">
-          <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
-      </None>
   </ItemGroup>
 
 </Project>

Interfaces/IComponent.cs

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

Interfaces/Port.cs

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

Program.cs

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

Scenarios/HelmholtzScenario.cs

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

Scenarios/Scenario.cs

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

Scenarios/Simulation.cs

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

Scenarios/SingleCylScenario.cs

@@ -1,266 +0,0 @@
-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;
-
-        private float _maxThrottleArea;
-        private float intakePipeArea, exhaustPipeArea;
-        private const float MaxBrakeTorque = 30.0f;   // Nm at full load
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            // Throttle body diameter 44mm (typical for 250cc MX)
-            _maxThrottleArea = (float)Units.AreaFromDiameter(44 * Units.mm);
-
-            // ---- Crankshaft ----
-            crankshaft = new Crankshaft(2000);
-            crankshaft.Inertia = 0.02f;          // kg·m² (crank + flywheel)
-            crankshaft.FrictionConstant = 3.0f;  // Nm – bearings, rings, seals
-            crankshaft.FrictionViscous = 0.002f; // Nm/(rad/s) – oil windage
-
-            // ---- Cylinder (CRF250R) ----
-            float bore = 0.078f;          // 78 mm
-            float stroke = 0.0522f;       // 52.2 mm → 249.4 cc
-            float conRod = 0.1044f;       // 2× stroke
-            float compRatio = 13.5f;      // typical
-
-            // Valve events (high‑performance MX cam)
-            float ivo = 340f, ivc = 600f; // intake opens 20° BTDC (overlap), closes 60° ABDC
-            float evo = 120f, evc = 380f; // exhaust opens 60° BBDC, closes 20° ATDC
-
-            cylinder = new Cylinder(bore, stroke, conRod, compRatio,
-                                    ivo, ivc, evo, evc, crankshaft)
-            {
-                IntakeValveDiameter = 0.036f,   // 36 mm
-                IntakeValveLift = 0.0095f,      // 9.5 mm
-                ExhaustValveDiameter = 0.030f,  // 30 mm
-                ExhaustValveLift = 0.0085f      // 8.5 mm
-            };
-
-            // ---- Pipe system ----
-            int[] pipeStart = { 0, 10, 20 };
-            int[] pipeEnd   = { 10, 20, 70 };
-            int totalCells = pipeEnd[^1];
-            float[] area = new float[totalCells];
-            float[] dx   = new float[totalCells];
-
-            float intakeDia = 0.040f;     // 40 mm intake runner
-            float exhaustDia = 0.038f;    // 38 mm exhaust primary
-            intakePipeArea  = MathF.PI * 0.25f * intakeDia * intakeDia;
-            exhaustPipeArea = MathF.PI * 0.25f * exhaustDia * exhaustDia;
-
-            float intakeLenBefore = 0.15f;  // throttle body to plenum
-            float intakeLenRunner = 0.25f;  // plenum to valve
-            float exhaustLen      = 0.50f;  // exhaust length
-
-            for (int i = 0; i < totalCells; i++)
-            {
-                if (i < 10)
-                {
-                    area[i] = intakePipeArea; dx[i] = intakeLenBefore / 10f;
-                }
-                else if (i < 20)
-                {
-                    area[i] = intakePipeArea; dx[i] = intakeLenRunner / 10f;
-                }
-                else
-                {
-                    area[i] = exhaustPipeArea; 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(1.0e-3f, 101325f, 300f);   // 1 litre airbox
-            plenumInlet = intakePlenum.CreatePort();
-            plenumOutlet = intakePlenum.CreatePort();
-            exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // unused
-            colIn = exhaustCollector.CreatePort();
-            colOut = exhaustCollector.CreatePort();
-
-            // ---- Boundary system ----
-            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
-
-            throttleAreaIdx = 0;
-            plenumRunnerAreaIdx = 1;
-            intakeValveIdx = 2;
-            exhaustValveIdx = 3;
-
-            // Open ends (pipe area = pipe cross‑section)
-            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, intakePipeArea);
-            intakeOpenIdx = 0;
-            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, exhaustPipeArea);
-            exhaustOpenIdx = 1;
-
-            // Orifices
-            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.7f);      // throttle
-            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1.0f); // plenum→runner
-            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1.0f); // intake valve
-            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1.0f); // exhaust valve
-
-            orificeAreas = new float[4];
-            orificeAreas[plenumRunnerAreaIdx] = intakePipeArea; // runner cross‑section (fixed)
-
-            // ---- 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 = 10f };
-            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 10f };
-            reverb = new OutdoorExhaustReverb(sampleRate);
-
-            stepCount = 0;
-            Console.WriteLine("CRF250R engine ready.");
-        }
-
-        public override float Process()
-        {   
-            // Manual brake torque (0..30 Nm)
-            float loadTorque = Load * MaxBrakeTorque;
-            crankshaft.SetLoadTorque(loadTorque);
-
-            crankshaft.Step((float)dt);
-            cylinder.PreStep((float)dt);
-
-            float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.001f, 1f);
-            orificeAreas[throttleAreaIdx] = throttledArea;
-
-            orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
-            orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
-            boundaries.SetOrificeAreas(orificeAreas);
-
-            solver.Step();
-            stepCount++;
-
-            float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
-            float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
-
-            float exhaustDry = exhaustSound.Process(exhaustFlow);
-            float intakeDry  = intakeSound.Process(intakeFlow);
-
-            if (stepCount % 1000 == 0)
-            {
-                float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-                float crankDeg = (crankshaft.CrankAngle + cylinder.PhaseOffset) * 180f / MathF.PI % 720f;
-                Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1}°, RPM={rpm:F0}, CylP={cylinder.Pressure/1e5f:F2} bar");
-                Console.WriteLine($"  intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
-
-                var (r0L, u0L, p0L) = pipeSystem.GetInteriorStateLeft(0);
-                var (r0R, u0R, p0R) = pipeSystem.GetInteriorStateRight(0);
-                Console.WriteLine($"  Pipe0 L: rho={r0L:F4} u={u0L:F3} p={p0L/1e5:F3}bar | R: rho={r0R:F4} u={u0R:F3} p={p0R/1e5:F3}bar");
-
-                var (r1L, u1L, p1L) = pipeSystem.GetInteriorStateLeft(1);
-                var (r1R, u1R, p1R) = pipeSystem.GetInteriorStateRight(1);
-                Console.WriteLine($"  Pipe1 L: rho={r1L:F4} u={u1L:F3} p={p1L/1e5:F3}bar | R: rho={r1R:F4} u={u1R:F3} p={p1R/1e5:F3}bar");
-
-                var (r2L, u2L, p2L) = pipeSystem.GetInteriorStateLeft(2);
-                var (r2R, u2R, p2R) = pipeSystem.GetInteriorStateRight(2);
-                Console.WriteLine($"  Pipe2 L: rho={r2L:F4} u={u2L:F3} p={p2L/1e5:F3}bar | R: rho={r2R:F4} u={u2R:F3} p={p2R/1e5:F3}bar");
-
-                Console.WriteLine($"  Plenum P={intakePlenum.Pressure/1e5:F3}bar, mass={intakePlenum.Mass:E4} kg");
-                Console.WriteLine($"  Cyl mass={cylinder.Mass:E4} kg");
-            }
-
-            return reverb.Process((intakeDry + exhaustDry) * 0.5f);
-        }
-
-        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;
-
-            float pipe1StartX = openEndX;
-            float pipe1EndX = pipe1StartX + 120f;
-            DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
-
-            float throttleX = pipe1EndX + 5f;
-            var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
-            {
-                FillColor = Color.Yellow,
-                Position = new Vector2f(throttleX, intakeY - 15f)
-            };
-            target.Draw(throttleRect);
-
-            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);
-
-            float runnerStartX = plenLeftX + plenW + 5f;
-            float runnerEndX = runnerStartX + 100f;
-            DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
-
-            float cylCX = runnerEndX + 50f;
-            float cylTopY = intakeY - 120f;
-            float cylW = 80f, cylMaxH = 240f;
-            DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
-
-            float exhStartX = cylCX + cylW / 2f + 20f;
-            float exhEndX = winW - 60f;
-            DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
-
-            // --- RPM & Power labels ---
-            float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-            float powerKw = crankshaft.AveragePower * 1e-3f;
-            DrawLabel(target, $"RPM: {rpm:F0}", new Vector2f(20, 90), Color.White, 24);
-            DrawLabel(target, $"Power: {powerKw:F2} kW", new Vector2f(20, 115), Color.White, 24);
-
-            // --- Dyno curve ---
-            float torqueNm = crankshaft.AverageTorque;
-            UpdateDynoCurve(rpm, powerKw, torqueNm);
-
-            float graphX = winW - 410f;
-            float graphY = winH - 260f;
-            float graphW = 400f;
-            float graphH = 250f;
-            DrawDynoCurve(target, graphX, graphY, graphW, graphH, rpm, powerKw);
-        }
-    }
-}

Scenarios/TestScenario.cs

@@ -1,102 +0,0 @@
-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);
-        }
-    }
-}

Scenarios/TwoStrokeScenario.cs

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

Utilities/Units.cs

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