Attempting to generate sound, and set cells automatically

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

@@ -0,0 +1,43 @@
+using SFML.Audio;
+using SFML.System;
+
+namespace FluidSim.Audio
+{
+    internal class RingBufferStream : SoundStream
+    {
+        private readonly RingBuffer ringBuffer;
+
+        public RingBufferStream(RingBuffer buffer)
+        {
+            ringBuffer = buffer;
+            // 2 channels, 44.1 kHz, standard stereo mapping
+            Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
+        }
+
+        protected override bool OnGetData(out short[] samples)
+        {
+            const int monoBlockSize = 512;          // number of mono samples we'll read
+            float[] temp = new float[monoBlockSize];
+            int read = ringBuffer.Read(temp, monoBlockSize);
+            samples = new short[monoBlockSize * 2];
+
+            if (read > 0)
+            {
+                for (int i = 0; i < read; i++)
+                {
+                    float clamped = Math.Clamp(temp[i], -1f, 1f);
+                    short final = (short)(clamped * short.MaxValue);
+                    samples[i * 2] = final;   // left
+                    samples[i * 2 + 1] = final;   // right
+                }
+            }
+            for (int i = read * 2; i < samples.Length; i++)
+                samples[i] = 0;
+
+            return true;
+        }
+
+        protected override void OnSeek(Time timeOffset) =>
+            throw new NotSupportedException();
+    }
+}

Audio/SoundEngine.cs

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

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

@@ -12,29 +12,52 @@ namespace FluidSim.Components
         private int _n;
         private double _dx, _dt, _gamma = 1.4, _area;
         private double[] _rho, _rhou, _E;
+        private double _hydraulicDiameter;
 
-        // Volume states at boundaries
         private double _rhoLeft, _pLeft, _rhoRight, _pRight;
         private bool _leftBCSet, _rightBCSet;
 
-        public double FrictionFactor { get; set; } = 0.02;
+        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);
 
-        public Pipe1D(double length, double area, int nCells, int sampleRate)
+        /// <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)
         {
-            _n = nCells;
+            if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
-            _dx = length / nCells;
-            _dt = 1.0 / sampleRate;
             _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();
         }
@@ -56,7 +79,6 @@ namespace FluidSim.Components
         public double GetLeftDensity() => _rho[0];
         public double GetRightDensity() => _rho[_n - 1];
 
-        // ★ New: pass both density and pressure from the volume
         public void SetLeftVolumeState(double rhoVol, double pVol)
         {
             _rhoLeft = rhoVol;
@@ -88,22 +110,13 @@ namespace FluidSim.Components
             // --- Left boundary (face 0) ---
             if (_leftBCSet)
             {
-                // Ghost = actual volume state (ρ_vol, u=0, p_vol)
+                double rhoL = _rhoLeft, uL = 0.0, pL = _pLeft;
-                double rhoL = _rhoLeft;
+                double rhoR = _rho[0], uR = _rhou[0] / Math.Max(rhoR, 1e-12), pR = Pressure(0);
-                double uL = 0.0;
-                double pL = _pLeft;
-
-                double rhoR = _rho[0];
-                double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
-                double pR = Pressure(0);
-
                 HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]);
             }
             else
             {
-                Fm[0] = 0;
+                Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0;
-                Fp[0] = Pressure(0);
-                Fe[0] = 0;
             }
 
             // --- Internal faces ---
@@ -118,25 +131,16 @@ namespace FluidSim.Components
             // --- Right boundary (face n) ---
             if (_rightBCSet)
             {
-                double rhoL = _rho[n - 1];
+                double rhoL = _rho[n - 1], uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12), pL = Pressure(n - 1);
-                double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
+                double rhoR = _rhoRight, uR = 0.0, pR = _pRight;
-                double pL = Pressure(n - 1);
-
-                // Ghost = actual volume state (ρ_vol, u=0, p_vol)
-                double rhoR = _rhoRight;
-                double uR = 0.0;
-                double pR = _pRight;
-
                 HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]);
             }
             else
             {
-                Fm[n] = 0;
+                Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0;
-                Fp[n] = Pressure(n - 1);
-                Fe[n] = 0;
             }
 
-            // --- Cell update ---
+            // --- Cell update (inviscid fluxes) ---
             for (int i = 0; i < n; i++)
             {
                 double dM = (Fm[i + 1] - Fm[i]) / _dx;
@@ -149,12 +153,41 @@ namespace FluidSim.Components
                 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 disabled ---
+            // --- Friction (Darcy–Weisbach, energy‑conserving) ---
-            // if (FrictionFactor > 0) { … }
+            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];
 
-            // --- Port flows ---
             PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0;
             PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0;
 
@@ -170,14 +203,24 @@ namespace FluidSim.Components
         void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
                       out double fm, out double fp, out double fe)
         {
-            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12));
+            const double eps = 1e-12;
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 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 Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
+
-                      / (rL * (SL - uL) - rR * (SR - uR));
+            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;
@@ -186,16 +229,20 @@ namespace FluidSim.Components
             else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
             else if (Ss >= 0)
             {
-                double rsL = rL * (SL - uL) / (SL - Ss);
+                double diffSL = SL - uL;
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
+                if (Math.Abs(diffSL) < eps) diffSL = eps;
-                double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
+                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 rsR = rR * (SR - uR) / (SR - Ss);
+                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 * (SR - uR)));
+                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,6 +1,4 @@
-using System;
+using FluidSim.Interfaces;
-using FluidSim.Interfaces;
-using FluidSim.Utils;
 
 namespace FluidSim.Components
 {

Core/OrificeBoundary.cs

@@ -1,5 +1,4 @@
-using System;
+using FluidSim.Components;
-using FluidSim.Components;
 
 namespace FluidSim.Core
 {

Core/Simulation.cs

@@ -1,67 +0,0 @@
-using System;
-using FluidSim.Components;
-using FluidSim.Utils;
-
-namespace FluidSim.Core
-{
-    public static class Simulation
-    {
-        private static Solver solver;
-        private static Volume0D volA, volB;
-        private static Pipe1D pipe;
-        private static Connection connA, connB;
-        private static int stepCount;
-        private static double time;
-        private static double dt;
-
-        public static void Initialize(int sampleRate)
-        {
-            dt = 1.0 / sampleRate;
-
-            double V = 5.0 * Units.L;
-            volA = new Volume0D(V, 1.1 * Units.atm, Units.Celsius(20), sampleRate);
-            volB = new Volume0D(V, 1.0 * Units.atm, Units.Celsius(20), sampleRate);
-
-            double length = 150 * Units.mm;
-            double diameter = 25 * Units.mm;
-            double area = Units.AreaFromDiameter(25, Units.mm);
-            int nCells = 10;
-
-            pipe = new Pipe1D(length, area, nCells, sampleRate);
-            pipe.SetUniformState(volA.Density, 0.0, volA.Pressure);
-            pipe.FrictionFactor = 0.02;
-
-            // Connections with orifice area equal to pipe area (flange joint)
-            connA = new Connection(volA.Port, pipe.PortA) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
-            connB = new Connection(pipe.PortB, volB.Port) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 };
-
-            solver = new Solver();
-            solver.AddVolume(volA);
-            solver.AddVolume(volB);
-            solver.AddPipe(pipe);
-            solver.AddConnection(connA);
-            solver.AddConnection(connB);
-        }
-
-        public static float Process()
-        {
-            solver.Step();
-            time += dt;
-            stepCount++;
-            Log();
-            return 0f;
-        }
-
-        public static void Log()
-        {
-            if (stepCount <= 50 || stepCount % 200 == 0)
-            {
-                Console.WriteLine(
-                    $"t = {time * 1e3:F3} ms  Step {stepCount:D4}: " +
-                    $"PA = {volA.Pressure / 1e5:F6} bar, " +
-                    $"PB = {volB.Pressure / 1e5:F6} bar, " +
-                    $"FlowA = {pipe.PortA.MassFlowRate * 1e3:F2} g/s");
-            }
-        }
-    }
-}

Core/Solver.cs

@@ -1,4 +1,4 @@
-using System.Collections.Generic;
+using FluidSim.Audio;
 using FluidSim.Components;
 using FluidSim.Interfaces;
 
@@ -10,17 +10,51 @@ namespace FluidSim.Core
         private readonly List<Pipe1D> _pipes = new();
         private readonly List<Connection> _connections = new();
 
+        public float LastSample { get; private set; }
+
         public void AddVolume(Volume0D v) => _volumes.Add(v);
         public void AddPipe(Pipe1D p) => _pipes.Add(p);
         public void AddConnection(Connection c) => _connections.Add(c);
 
         public void Step()
         {
-            // 1. Volumes publish state to their ports
+            // 1. Publish volume states to their own ports
             foreach (var v in _volumes)
                 v.PushStateToPort();
 
-            // 2. Set volume states as boundary conditions on pipes
+            // 2. Handle direct volume‑to‑volume connections
+            foreach (var conn in _connections)
+            {
+                if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB))
+                {
+                    Volume0D volA = _volumes.Find(v => v.Port == conn.PortA);
+                    Volume0D volB = _volumes.Find(v => v.Port == conn.PortB);
+
+                    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))
@@ -29,11 +63,11 @@ namespace FluidSim.Core
                     SetVolumeBC(conn.PortB, conn.PortA);
             }
 
-            // 3. Run pipe simulations
+            // 4. Run pipe simulations
             foreach (var p in _pipes)
                 p.Simulate();
 
-            // 4. Transfer pipe‑port flows to volume ports
+            // 5. Transfer pipe‑to‑volume flows – unchanged
             foreach (var conn in _connections)
             {
                 if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
@@ -42,9 +76,21 @@ namespace FluidSim.Core
                     TransferPipeToVolume(conn.PortB, conn.PortA);
             }
 
-            // 5. Integrate volumes
+            // 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);

Core/SoundEngine.cs

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

Program.cs

@@ -2,7 +2,8 @@
 using SFML.Window;
 using SFML.System;
 using System.Diagnostics;
-using FluidSim.Core;
+using FluidSim.Scenarios;
+using FluidSim.Audio;
 
 namespace FluidSim;
 
@@ -10,6 +11,8 @@ public class Program
 {
     private const int SampleRate = 44100;
     private static volatile bool running = true;
+    // Global step counter – incremented every simulation step
+    private static long stepCount = 0;
 
     public static void Main()
     {
@@ -22,47 +25,71 @@ public class Program
         soundEngine.Volume = 70;
         soundEngine.Start();
 
-        double lastAudioTime = 0.0;
         var stopwatch = Stopwatch.StartNew();
 
-        int warmupSamples = SampleRate / 2;
+        // --- 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);
-        lastAudioTime = stopwatch.Elapsed.TotalSeconds;
+
+        // 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];
 
-        Simulation.Initialize(SampleRate);
+        double lastLogTime = 0.0;                         // for periodic speed printout
 
         while (window.IsOpen)
         {
             window.DispatchEvents();
 
+            // --- Compute how many audio samples are needed since last frame ---
             double currentTime = stopwatch.Elapsed.TotalSeconds;
-            double elapsed = currentTime - lastAudioTime;
+            double elapsed = currentTime;                  // since stopwatch was reset
-            int samplesNeeded = (int)(elapsed * SampleRate);
+            int samplesNeeded = (int)(elapsed * SampleRate) - (int)(stepCount);
+            // (stepCount is total generated samples, so we just need the remainder)
 
+            // --- Generate the required number of simulation steps ---
             while (samplesNeeded > 0 && running)
             {
                 int toGenerate = Math.Min(samplesNeeded, chunkSize);
                 for (int i = 0; i < toGenerate; i++)
                 {
                     buffer[i] = Simulation.Process();
+                    stepCount++;
                 }
                 soundEngine.WriteSamples(buffer, toGenerate);
                 samplesNeeded -= toGenerate;
             }
 
-            lastAudioTime = currentTime;
+            // --- 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();
         }
 
+        // --- Cleanup ---
         soundEngine.Dispose();
         window.Dispose();
     }

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

Sources/EffortSource.cs

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

Sources/FlowSource.cs

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

Utilities/Units.cs

@@ -1,6 +1,4 @@
-using System;
+namespace FluidSim.Utilities
-
-namespace FluidSim.Utils
 {
     public static class Units
     {