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;
 
@@ -50,4 +40,4 @@ namespace FluidSim.Audio
         protected override void OnSeek(Time timeOffset) =>
             throw new NotSupportedException();
     }
-}
+}

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

Components/Cylinder.cs

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

Components/Pipe1D.cs

@@ -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/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,370 +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) – inertance only ---
-            // If 0 when UseInertance is true, a stable default is auto‑computed at runtime.
-            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) ----------
-        // Simple isentropic nozzle – no built‑in loss. For dissipation use pipe damping
-        // or the inertance model if you need a damped resonator.
-        public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
-                               int areaIndex, float dischargeCoeff = 1f)
-        {
-            _orifices[OrificeCount] = new OrificeDesc
-            {
-                VolumePort = volumePort,
-                PipeIndex = pipeIndex,
-                IsLeftEnd = isLeftEnd,
-                AreaIndex = areaIndex,
-                DischargeCoeff = dischargeCoeff,
-                UseInertance = false,
-                EffectiveLength = 0f,
-                CurrentMdot = 0f,
-                LossCoefficient = 0f
-            };
-            OrificeCount++;
-        }
-
-        // ---------- Add orifice with inertance ----------
-        // effectiveLength – length of the inertial slug (m), typically the physical neck length.
-        // lossCoefficient – linear resistance (N·s/m⁵). If 0 (or omitted) an automatic stable
-        //                   value will be computed from the pipe's characteristic impedance.
-        public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
-                                            int areaIndex, float dischargeCoeff,
-                                            float effectiveLength, float lossCoefficient = 0f)
-        {
-            AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff);
-            ref var d = ref _orifices[OrificeCount - 1];
-            d.UseInertance = true;
-            d.EffectiveLength = effectiveLength;
-            d.LossCoefficient = lossCoefficient;
-        }
-
-        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 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 (for reference) ----
-                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 ghost state ----
-                float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
-
-                if (d.UseInertance)
-                {
-                    // ---- Inertance ODE with (possibly automatic) linear loss ----
-                    float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
-                    float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
-                    float dp = volP - pipeP;
-
-                    // If loss coefficient not provided, use a tiny fraction of the pipe's characteristic impedance
-                    float Rlin = d.LossCoefficient;
-                    if (Rlin <= 0f)
-                    {
-                        // Auto‑sized linear drag: 0.5 % of Z_char
-                        float rhoRef = d.CurrentMdot >= 0 ? volRho : pipeRho;
-                        float cRef = d.CurrentMdot >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
-                        float Z_char = rhoRef * cRef / MathF.Max(area, 1e-12f);
-                        Rlin = 0.005f * Z_char;
-                    }
-
-                    float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
-                    float mdotNew = d.CurrentMdot + dmdot_dt * dt;
-
-                    // Symmetric flow limiters
-                    if (d.VolumePort.Owner is Volume0D vol0)
-                    {
-                        float maxOut = vol0.Mass / dt;
-                        if (mdotNew > maxOut) mdotNew = maxOut;
-                        if (mdotNew < -maxOut) mdotNew = -maxOut;
-                    }
-
-                    int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
-                                              : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
-                    float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
-                    float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
-                    float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
-                    float maxFromPipe = pipeCellMass / dt;
-                    if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
-
-                    // Velocity clamp 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 (purely isentropic) ----
-                    mdotFinal = mdotEst;
-                    rhoFace = rhoFaceEst;
-                    uFace = uFaceEst;
-                    pFace = pFaceEst;
-
-                    // Limit outflow from cavity
-                    if (d.VolumePort.Owner is Volume0D vol0)
-                    {
-                        float maxOut = vol0.Mass / dt;
-                        if (mdotFinal > maxOut) mdotFinal = maxOut;
-                    }
-
-                    // Safety velocity clamp (Mach 0.9)
-                    float cLocal = mdotFinal >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
-                    float maxULocal = 0.9f * cLocal;
-                    float uCheck = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
-                    if (uCheck > maxULocal)
-                    {
-                        uFace = maxULocal;
-                        mdotFinal = rhoFace * uFace * area * (mdotFinal >= 0 ? 1f : -1f);
-                    }
-
-                    d.CurrentMdot = mdotFinal;
-                }
-
-                // ---- Determine air fraction for ghost ----
-                if (mdotFinal >= 0)
-                {
-                    airFracGhost = volAF;
-                }
-                else
-                {
-                    airFracGhost = pipeAF;
-                    if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
-                }
-
-                // ---- Apply sign convention for velocity ----
-                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 flowing into the volume
-                    {
-                        float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
-                        d.VolumePort.SpecificEnthalpy = pipeH;
-                    }
-                    else   // mass flowing out of the volume
-                    {
-                        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;
-
-                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);
-
-                bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
-                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;
-                }
-
-                if (d.IsLeftEnd)
-                    _pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
-                else
-                    _pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
-
-                float area = d.PipeArea;
-                float mdot = rhoGhost * uGhost * area;
-                if (d.IsLeftEnd) mdot = -mdot;
-                d.LastMassFlowRate = mdot;
-                d.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,561 +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, _c is gone
-        private float[] _p;
-
-        // Flux arrays (size = _allCells + 1)
-        private float[] _fluxM, _fluxP, _fluxE, _fluxY;
-
-        // Damping and relaxation (computed on‑the‑fly only if used)
-        private float[] _dampingFactors;
-        private float[] _relaxFactors;
-        private bool _applyDamping;
-        private bool _applyRelax;
-
-        // Ghost buffer
-        private readonly GhostBuffer _ghost;
-
-        // Wall mask – precomputed once
-        private readonly bool[] _isWallFace;
-
-        // ---------- 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];
-
-            _dampingFactors = new float[_allCells];
-            _relaxFactors = new float[_allCells];
-            _applyDamping = _coeffBase != 0f;
-            _applyRelax = _relaxRate != 0f;
-
-            _ghost = new GhostBuffer(_pipeCount);
-            _ambientEnergyRef = initialP * Gm1Inv;
-
-            // Pre‑compute wall face flags: each face that sits between two different pipes is a wall
-            _isWallFace = new bool[faceCount];
-            for (int f = 1; f < _totalCells; f++)
-            {
-                for (int p = 0; p < _pipeCount; p++)
-                {
-                    if (f == _pipeEnd[p] && f < _totalCells)
-                    {
-                        _isWallFace[f] = true;
-                        break;
-                    }
-                }
-            }
-
-            // 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 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: fuses primitive calculation and flux evaluation ----------
-        private void ComputeFluxes(float dt)
-        {
-            float fm, fp, fe;
-            int vecSize = Vector<float>.Count;
-
-            // ---- 1. Left ghost boundaries ----
-            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];
-
-                // store pressure for cell idx
-                _p[idx] = pR;
-
-                LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
-                _fluxM[idx] = fm; _fluxP[idx] = fp; _fluxE[idx] = fe;
-
-                float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
-                ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
-                _fluxY[idx] = fy;
-            }
-
-            // ---- 2. Right ghost boundaries ----
-            for (int p = 0; p < _pipeCount; p++)
-            {
-                int idx = _pipeEnd[p] - 1;
-                int face = idx + 1;
-                int ghostIdx = p * 2 + 1;
-                float rR = _ghost.Rho[ghostIdx];
-                float uR = _ghost.U[ghostIdx];
-                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];
-
-                // store pressure for cell idx
-                _p[idx] = pL;
-
-                LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
-                _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
-
-                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;
-            }
-
-            // ---- 3. Interior faces – vectorised SIMD ----
-            for (int face = 1; face < _totalCells; face++)
-            {
-                // Handle walls (rare) with scalar code
-                if (_isWallFace[face])
-                {
-                    int iL = face - 1;
-                    float rL = _rho[iL], rhouL = _rhou[iL];
-                    float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
-                    float uL = rhouL * invRhoL;
-                    float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
-                    float cL = MathF.Sqrt(Gamma * pL * invRhoL);
-                    _p[iL] = pL;
-
-                    LaxFlux(rL, uL, pL, cL, rL, -uL, pL, cL, out fm, out fp, out fe);
-                    _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
-                    _fluxY[face] = 0f;
-                    continue;
-                }
-
-                // If the next vecSize faces contain a wall, fall back to scalar for this block
-                if (face + vecSize - 1 < _totalCells)
-                {
-                    bool hasWall = false;
-                    for (int f = face; f < face + vecSize; f++)
-                        if (_isWallFace[f]) { hasWall = true; break; }
-
-                    if (!hasWall)
-                    {
-                        // --- 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);
-
-                        // Store pressures for visualisation (left cell of each face)
-                        pL.CopyTo(_p, face - 1);
-
-                        // Lax‑Friedrichs fluxes
-                        var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass
-                        var 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, so we advance vecSize faces
-                        continue;
-                    }
-                }
-
-                // --- Scalar interior face (fallback) ---
-                {
-                    int iLf = face - 1, iRf = face;
-                    float rLf = _rho[iLf], rhouLf = _rhou[iLf];
-                    float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f));
-                    float uLf = rhouLf * invRhoLf;
-                    float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf);
-                    float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf);
-                    float YLf = _Y[iLf];
-                    _p[iLf] = pLf;
-
-                    float rRf = _rho[iRf], rhouRf = _rhou[iRf];
-                    float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f));
-                    float uRf = rhouRf * invRhoRf;
-                    float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf);
-                    float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf);
-                    float YRf = _Y[iRf];
-
-                    LaxFlux(rLf, uLf, pLf, cLf, rRf, uRf, pRf, cRf, out fm, out fp, out fe);
-                    _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
-
-                    float alpha = MathF.Max(MathF.Abs(uLf) + cLf, MathF.Abs(uRf) + cRf);
-                    ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy);
-                    _fluxY[face] = fy;
-                }
-            }
-
-            // If damping/relaxation are active, compute the factors here (re-uses _dampingFactors/_relaxFactors arrays,
-            // but we no longer have a separate precompute pass). We compute them on demand in UpdateCells anyway?
-            // Actually UpdateCells multiplies by these factors; we can compute them there if needed.
-        }
-
-        // ---------- Cell update (unchanged core, but skips relaxation/damping when not needed) ----------
-        private void UpdateCells(float dt)
-        {
-            int vecSize = Vector<float>.Count;
-            float dtRelax = -_relaxRate * dt;
-
-            // Compute damping and relaxation factors if needed
-            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)
-            {
-                var relaxVal = MathF.Exp(dtRelax);
-                for (int i = 0; i < _totalCells; i++)
-                    _relaxFactors[i] = relaxVal;
-            }
-
-            int iCell = 0;
-            for (; iCell <= _totalCells - 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 (only a few cells)
-            for (; iCell < _totalCells; 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);
-            }
-        }
-
-        // ---------- Scalar flux helpers (used in boundaries and scalar fallback) ----------
-        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);
-        }
-
-        // ---------- 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,93 +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 = SubStepCount;
+            // 2. Handle direct volume‑to‑volume connections
-            float dtSub = (float)(_dt / nSub);
+            foreach (var conn in _connections)
-
-            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,216 +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!;   // cast to access ThrottleArea
-    private static Font? _overlayFont;
-    private static Text? _overlayText;
-
-    // Throttle control
-    private static float _throttleTarget = 1.0f;          // 0‑1, set by arrow keys
-    private static float _throttleCurrent = 0.0f;         // actual current fraction (lerped)
-    private const float ThrottleLerpRate = 10.0f;          // times per second (speed of movement)
-    private static bool _wKeyHeld = false;
-    private static float _lastThrottleUpdateTime;
-
-    private const int TargetMaxFill = (int)(SampleRate * 0.2);
 
     public static void Main()
     {
-        var window = CreateWindow();
+        var mode = new VideoMode(new Vector2u(1280, 720));
-        LoadFont();
+        var window = new RenderWindow(mode, "Fluid Simulation");
-        _scenario = new SingleCylScenario();
+        window.SetVerticalSyncEnabled(true);
-        _scenario.Initialize(SampleRate);
+        window.Closed += (_, _) => { running = false; window.Close(); };
-        _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 update ----
-            float dtThrottle = (float)now - _lastThrottleUpdateTime;
-            _lastThrottleUpdateTime = (float)now;
-
-            float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
-
-            // Snap to zero instantly when target is zero (key released)
-            if (throttleDesiredFraction == 0.0)
             {
-                _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;
-            }
-
-            _scenario.Throttle = _throttleCurrent;
-
-            // ---- 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" +
-                        $"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;
-        }
-    }
-
-    private static void OnKeyReleased(object? sender, KeyEventArgs e)
-    {
-        if (e.Code == Keyboard.Key.W)
-            _wKeyHeld = 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,161 +0,0 @@
-using SFML.Graphics;
-using SFML.System;
-using FluidSim.Core;
-using FluidSim.Components;
-
-namespace FluidSim.Tests
-{
-    public abstract class Scenario
-    {
-        protected const float AmbientPressure = 101325f;
-        protected const float AmbientTemperature = 300f;
-        public float Throttle { get; set; }
-
-        public abstract void Initialize(int sampleRate);
-        public abstract float Process();
-        public abstract void Draw(RenderWindow target);
-
-        protected Color PressureColor(float pressurePa)
-        {
-            float bar = pressurePa / 1e5f;
-            byte r, g, b;
-            if (bar < 1f)
-            {
-                float f = Math.Clamp(bar, 0f, 1f);
-                r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
-            }
-            else
-            {
-                float f = Math.Min((bar - 1f) / 9f, 1f);
-                r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
-            }
-            return new Color(r, g, b);
-        }
-
-        protected Color TemperatureColor(float t)
-        {
-            t = Math.Clamp(t, 0f, 2000f);
-            byte r, g, b;
-            if (t < AmbientTemperature)
-            {
-                float f = t / AmbientTemperature;
-                r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
-            }
-            else
-            {
-                float f = (t - AmbientTemperature) / (2000f - AmbientTemperature);
-                r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
-            }
-            return new Color(r, g, b);
-        }
-
-        protected void DrawVolume(RenderWindow target, Volume0D volume,
-                                  float centerX, float topY, float width, float height)
-        {
-            var rect = new RectangleShape(new Vector2f(width, height))
-            {
-                FillColor = PressureColor(volume.Pressure),
-                Position = new Vector2f(centerX - width / 2f, topY)
-            };
-            target.Draw(rect);
-            var border = new RectangleShape(new Vector2f(width, height))
-            {
-                FillColor = Color.Transparent,
-                OutlineColor = Color.White,
-                OutlineThickness = 1f,
-                Position = rect.Position
-            };
-            target.Draw(border);
-        }
-
-        protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
-                                    float centerX, float topY, float width, float maxHeight)
-        {
-            float fraction = cylinder.PistonFraction;
-            float currentHeight = maxHeight * fraction;
-            var wall = new RectangleShape(new Vector2f(width, maxHeight))
-            {
-                FillColor = new Color(60, 60, 60),
-                Position = new Vector2f(centerX - width / 2f, topY)
-            };
-            target.Draw(wall);
-            var gas = new RectangleShape(new Vector2f(width, currentHeight))
-            {
-                FillColor = PressureColor(cylinder.Pressure),
-                Position = new Vector2f(centerX - width / 2f, topY)
-            };
-            target.Draw(gas);
-            var piston = new RectangleShape(new Vector2f(width, 4f))
-            {
-                FillColor = Color.White,
-                Position = new Vector2f(centerX - width / 2f, topY + currentHeight)
-            };
-            target.Draw(piston);
-            float valveW = 6f, valveH = 10f, valveY = topY + 4f;
-            var iv = new RectangleShape(new Vector2f(valveW, valveH))
-            {
-                FillColor = cylinder.IntakeValveArea > 0f ? Color.Green : Color.Red,
-                Position = new Vector2f(centerX - width / 2f - valveW - 2f, valveY)
-            };
-            target.Draw(iv);
-            var ev = new RectangleShape(new Vector2f(valveW, valveH))
-            {
-                FillColor = cylinder.ExhaustValveArea > 0f ? Color.Green : Color.Red,
-                Position = new Vector2f(centerX + width / 2f + 2f, valveY)
-            };
-            target.Draw(ev);
-        }
-
-        protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex,
-                                float pipeCenterY, float pipeStartX, float pipeEndX)
-        {
-            int start = pipeSystem.GetPipeStart(pipeIndex);
-            int end = pipeSystem.GetPipeEnd(pipeIndex);
-            int n = end - start;
-            if (n < 2) return;
-
-            float pipeLen = pipeEndX - pipeStartX;
-            float dx = pipeLen / (n - 1);
-            float baseRadius = 25f;
-
-            var centers = new float[n];
-            var radii = new float[n];
-            var temps = new float[n];
-            for (int i = 0; i < n; i++)
-            {
-                int cell = start + i;
-                float p = pipeSystem.GetCellPressure(cell);
-                float rho = pipeSystem.GetCellDensity(cell);
-                temps[i] = p / MathF.Max(rho * 287f, 1e-12f);
-                float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
-                radii[i] = baseRadius * (1f + dev * 2f);
-                if (radii[i] < 2f) radii[i] = 2f;
-                centers[i] = pipeStartX + i * dx;
-            }
-
-            int segments = 8;
-            var va = new VertexArray(PrimitiveType.TriangleStrip);
-            for (int i = 0; i < n; i++)
-            {
-                float x = centers[i], r = radii[i];
-                Color col = TemperatureColor(temps[i]);
-                va.Append(new Vertex(new Vector2f(x, pipeCenterY - r), col));
-                va.Append(new Vertex(new Vector2f(x, pipeCenterY + r), col));
-                if (i < n - 1)
-                {
-                    for (int s = 1; s <= segments; s++)
-                    {
-                        float t = s / (float)segments;
-                        float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
-                        float ri = radii[i] + (radii[i + 1] - radii[i]) * t;
-                        float Ti = temps[i] + (temps[i + 1] - temps[i]) * t;
-                        Color colS = TemperatureColor(Ti);
-                        va.Append(new Vertex(new Vector2f(xi, pipeCenterY - ri), colS));
-                        va.Append(new Vertex(new Vector2f(xi, pipeCenterY + ri), colS));
-                    }
-                }
-            }
-            target.Draw(va);
-        }
-    }
-}

Scenarios/Simulation.cs

@@ -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,272 +0,0 @@
-using FluidSim.Components;
-using FluidSim.Core;
-using FluidSim.Interfaces;
-using SFML.Graphics;
-using SFML.System;
-using System;
-
-namespace FluidSim.Tests
-{
-    public class SingleCylScenario : Scenario
-    {
-        // ---------- Engine components ----------
-        private Crankshaft crankshaft;
-        private Cylinder cylinder;
-
-        // ---------- Fluid network ----------
-        private PipeSystem pipeSystem;
-        private BoundarySystem boundaries;
-        private Solver solver;
-
-        // Volumes
-        private Volume0D intakePlenum;
-
-        // Ports
-        private Port plenumInlet, plenumOutlet;
-
-        // Orifice / open‑end indices
-        private int throttleAreaIdx, plenumRunnerIdx, intakeValveIdx, exhaustValveIdx;
-        private int intakeOpenIdx, exhaustOpenIdx;
-        private float[] orificeAreas;
-
-        // Sound
-        private SoundProcessor exhaustSound, intakeSound;
-        private OutdoorExhaustReverb reverb;
-
-        // ---------- Simulation state ----------
-        private double dt;
-        private int stepCount;
-
-        // ---------- Geometry (Lifan YX140) ----------
-        // Bore 56 mm, Stroke 57 mm, CR 9.5
-        private const float Bore = 0.056f;
-        private const float Stroke = 0.057f;
-        private const float ConRod = 0.110f;           // typical for 57 mm stroke
-        private const float CompressionRatio = 9.5f;
-
-        // Valve diameters (intake 27 mm, exhaust 23 mm)
-        private const float IntakeValveDiam = 0.027f;
-        private const float ExhaustValveDiam = 0.023f;
-        private const float ValveLift = 0.006f;         // 6 mm peak lift
-
-        // Valve timings (degrees, 720° four‑stroke)
-        // Intake: 15° BTDC → 45° ABDC
-        private const float IVO = 345f;   // 15° BTDC
-        private const float IVC = 585f;   // 45° ABDC (180°+45°)
-        // Exhaust: 45° BBDC → 15° ATDC
-        private const float EVO = 135f;   // 45° BBDC (180°-45°)
-        private const float EVC = 375f;   // 15° ATDC (360°+15°)
-
-        // Spark advance: 30° BTDC
-        private const float SparkAdv = 30f;
-
-        // Pipe / plenum sizes
-        private const float PipeDiam = 0.025f;          // 25 mm intake / exhaust
-        private const float PipeArea = 0.00049087f;      // π*D²/4
-        private const float PlenumVolume = 0.0005f;      // 500 mL
-        private const float MaxThrottleArea = 1e-4f;     // ~1 cm² (fully open)
-
-        // Pipe lengths and cell counts
-        private const float IntakeLenBefore = 0.15f;     // 15 cm before throttle
-        private const float RunnerLen = 0.25f;           // 25 cm runner
-        private const float ExhaustLen = 0.60f;          // 60 cm exhaust
-        private const int CellsBefore = 6;
-        private const int CellsRunner = 10;
-        private const int CellsExhaust = 24;
-
-        public override void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            // ---------- Crankshaft ----------
-            crankshaft = new Crankshaft(600);            // start at ~600 RPM
-            crankshaft.Inertia = 0.2f;
-            crankshaft.FrictionConstant = 2.0f;
-            crankshaft.FrictionViscous = 0.04f;
-
-            // ---------- Cylinder ----------
-            cylinder = new Cylinder(Bore, Stroke, ConRod, CompressionRatio,
-                                    IVO, IVC, EVO, EVC, crankshaft)
-            {
-                IntakeValveDiameter = IntakeValveDiam,
-                ExhaustValveDiameter = ExhaustValveDiam,
-                IntakeValveLift = ValveLift,
-                ExhaustValveLift = ValveLift,
-                SparkAdvance = SparkAdv,
-                EnergyVariationFraction = 0.03f,   // small cycle‑to‑cycle variation
-                MisfireProbability = 0.0f
-            };
-
-            // ---------- Pipe system ----------
-            int totalCells = CellsBefore + CellsRunner + CellsExhaust;
-            int[] pipeStart = { 0, CellsBefore, CellsBefore + CellsRunner };
-            int[] pipeEnd   = { CellsBefore, CellsBefore + CellsRunner, totalCells };
-
-            float[] areas = new float[totalCells];
-            float[] dxs   = new float[totalCells];
-            float dxBefore = IntakeLenBefore / CellsBefore;
-            float dxRunner = RunnerLen / CellsRunner;
-            float dxExh    = ExhaustLen / CellsExhaust;
-
-            for (int i = 0; i < totalCells; i++)
-            {
-                areas[i] = PipeArea;
-                if (i < CellsBefore)
-                    dxs[i] = dxBefore;
-                else if (i < CellsBefore + CellsRunner)
-                    dxs[i] = dxRunner;
-                else
-                    dxs[i] = dxExh;
-            }
-
-            float rho0 = 101325f / (287f * 300f);
-            pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, areas, dxs,
-                                        rho0, 0f, 101325f);
-            pipeSystem.DampingMultiplier = 0.5f;
-            pipeSystem.EnergyRelaxationRate = 0f;        // adiabatic pipes
-            pipeSystem.AmbientPressure = 101325f;
-
-            // ---------- Volumes ----------
-            intakePlenum = new Volume0D(PlenumVolume, 101325f, 300f);
-            plenumInlet  = intakePlenum.CreatePort();
-            plenumOutlet = intakePlenum.CreatePort();
-
-            // ---------- Boundary system ----------
-            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
-
-            throttleAreaIdx  = 0;
-            plenumRunnerIdx  = 1;
-            intakeValveIdx   = 2;
-            exhaustValveIdx  = 3;
-
-            // Open ends
-            boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, PipeArea);
-            intakeOpenIdx = 0;
-            boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, PipeArea);
-            exhaustOpenIdx = 1;
-
-            // Orifices
-            // throttle – variable area, low discharge for restriction
-            boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false,
-                                  throttleAreaIdx, dischargeCoeff: 0.8f);
-            // plenum → runner
-            boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true,
-                                  plenumRunnerIdx, dischargeCoeff: 1.0f);
-            // intake valve
-            boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false,
-                                  intakeValveIdx, dischargeCoeff: 1.0f);
-            // exhaust valve
-            boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true,
-                                  exhaustValveIdx, dischargeCoeff: 1.0f);
-
-            orificeAreas = new float[4];
-            orificeAreas[plenumRunnerIdx] = PipeArea;   // fixed full‑bore
-
-            // ---------- Solver ----------
-            solver = new Solver { SubStepCount = 5, EnableProfiling = false };
-            solver.SetTimeStep(dt);
-            solver.SetPipeSystem(pipeSystem);
-            solver.SetBoundarySystem(boundaries);
-            solver.AddComponent(cylinder);
-            solver.AddComponent(intakePlenum);
-
-            // ---------- Sound ----------
-            exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
-            intakeSound  = new SoundProcessor(sampleRate, 1f) { Gain = 0.2f };
-            reverb = new OutdoorExhaustReverb(sampleRate);
-
-            stepCount = 0;
-            Console.WriteLine("Single‑cylinder engine (YX140) ready.");
-        }
-
-        public override float Process()
-        {
-            // ---- Crank and cylinder pre‑step ----
-            crankshaft.Step((float)dt);
-            cylinder.PreStep((float)dt);
-
-            // ---- Update variable areas ----
-            float throttledArea = MaxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1.0f);
-            orificeAreas[throttleAreaIdx] = throttledArea;
-            orificeAreas[intakeValveIdx]  = cylinder.IntakeValveArea;
-            orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
-            boundaries.SetOrificeAreas(orificeAreas);
-
-            // ---- Fluids step ----
-            solver.Step();
-            stepCount++;
-
-            // ---- Sound ----
-            float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
-            float intakeFlow  = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
-
-            float exhaustDry = exhaustSound.Process(exhaustFlow);
-            float intakeDry  = intakeSound.Process(intakeFlow);
-
-            if (stepCount % 2000 == 0)
-            {
-                float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
-                Console.WriteLine($"Step {stepCount}, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar, " +
-                                  $"Throttle={Throttle * 100:F0}%");
-            }
-
-            return reverb.Process(exhaustDry + intakeDry);
-        }
-
-        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 leftX = 40f;
-
-            // Intake open end marker
-            var om = new CircleShape(5f) { FillColor = Color.Cyan };
-            om.Position = new Vector2f(leftX - 5f, intakeY - 5f);
-            target.Draw(om);
-
-            // Pipe 0 – before throttle
-            float p0EndX = leftX + 80f;
-            DrawPipe(target, pipeSystem, 0, intakeY, leftX, p0EndX);
-
-            // Throttle symbol
-            float thrX = p0EndX + 5f;
-            var thr = new RectangleShape(new Vector2f(8f, 30f))
-            {
-                FillColor = Color.Yellow,
-                Position = new Vector2f(thrX, intakeY - 15f)
-            };
-            target.Draw(thr);
-
-            // Plenum volume
-            float plenW = 60f, plenH = 50f;
-            float plenLeftX = thrX + 12f;
-            float plenCenterX = plenLeftX + plenW / 2f;
-            float plenTopY = intakeY - plenH / 2f;
-            DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
-
-            // Pipe 1 – runner
-            float rStartX = plenLeftX + plenW + 10f;
-            float rEndX = rStartX + 100f;
-            DrawPipe(target, pipeSystem, 1, intakeY, rStartX, rEndX);
-
-            // Cylinder
-            float cylCX = rEndX + 50f;
-            float cylTopY = intakeY - 120f;
-            float cylW = 80f, cylMaxH = 240f;
-            DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
-
-            // Pipe 2 – exhaust
-            float exhStartX = cylCX + cylW / 2f + 20f;
-            float exhEndX = winW - 60f;
-            DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
-
-            // Exhaust open end
-            var em = new CircleShape(5f) { FillColor = Color.Magenta };
-            em.Position = new Vector2f(exhEndX - 5f, exhaustY - 5f);
-            target.Draw(em);
-        }
-    }
-}

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

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