refactoring (broken right now)

.vscode/launch.json

@@ -0,0 +1,26 @@
+{
+    "version": "0.2.0",
+    "configurations": [
+        {
+            // Use IntelliSense to find out which attributes exist for C# debugging
+            // Use hover for the description of the existing attributes
+            // For further information visit https://github.com/dotnet/vscode-csharp/blob/main/debugger-launchjson.md
+            "name": ".NET Core Launch (console)",
+            "type": "coreclr",
+            "request": "launch",
+            "preLaunchTask": "build",
+            // If you have changed target frameworks, make sure to update the program path.
+            "program": "${workspaceFolder}/bin/Debug/net10.0/FluidSim.dll",
+            "args": [],
+            "cwd": "${workspaceFolder}",
+            // For more information about the 'console' field, see https://aka.ms/VSCode-CS-LaunchJson-Console
+            "console": "internalConsole",
+            "stopAtEntry": false
+        },
+        {
+            "name": ".NET Core Attach",
+            "type": "coreclr",
+            "request": "attach"
+        }
+    ]
+}

.vscode/tasks.json

@@ -0,0 +1,41 @@
+{
+    "version": "2.0.0",
+    "tasks": [
+        {
+            "label": "build",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "build",
+                "${workspaceFolder}/FluidSim.csproj",
+                "/property:GenerateFullPaths=true",
+                "/consoleloggerparameters:NoSummary;ForceNoAlign"
+            ],
+            "problemMatcher": "$msCompile"
+        },
+        {
+            "label": "publish",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "publish",
+                "${workspaceFolder}/FluidSim.csproj",
+                "/property:GenerateFullPaths=true",
+                "/consoleloggerparameters:NoSummary;ForceNoAlign"
+            ],
+            "problemMatcher": "$msCompile"
+        },
+        {
+            "label": "watch",
+            "command": "dotnet",
+            "type": "process",
+            "args": [
+                "watch",
+                "run",
+                "--project",
+                "${workspaceFolder}/FluidSim.csproj"
+            ],
+            "problemMatcher": "$msCompile"
+        }
+    ]
+}

Components/Atmosphere.cs

@@ -0,0 +1,43 @@
+using FluidSim.Interfaces;
+
+namespace FluidSim.Components
+{
+    /// <summary>
+    /// Represents the ambient atmosphere – constant pressure/temperature reservoir.
+    /// </summary>
+    public class Atmosphere : IComponent
+    {
+        public double Pressure { get; set; } = 101325.0;
+        public double Temperature { get; set; } = 300.0;
+        public double GasConstant { get; set; } = 287.0;
+        public double Gamma => 1.4;
+
+        public double Density => Pressure / (GasConstant * Temperature);
+        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
+
+        public Port Port { get; }
+
+        public Atmosphere()
+        {
+            Port = new Port { Owner = this };
+            UpdatePort();
+        }
+
+        public IReadOnlyList<Port> Ports => new[] { Port };
+
+        public void UpdateState(double dt)
+        {
+            // Atmosphere is static – just ensure the port reflects current values
+            UpdatePort();
+        }
+
+        private void UpdatePort()
+        {
+            Port.Pressure = Pressure;
+            Port.Density = Density;
+            Port.Temperature = Temperature;
+            Port.SpecificEnthalpy = SpecificEnthalpy;
+            // MassFlowRate is set by the solver connector
+        }
+    }
+}

Components/Crankshaft.cs

@@ -0,0 +1,54 @@
+// Components/Crankshaft.cs
+using System;
+
+namespace FluidSim.Components
+{
+    public class Crankshaft
+    {
+        public double AngularVelocity { get; set; }   // rad/s
+        public double CrankAngle { get; set; }        // rad, 0 … 4π (four‑stroke cycle)
+        public double PreviousAngle { get; set; }     // ← now has public setter
+
+        public double Inertia { get; set; } = 0.2;
+        public double FrictionConstant { get; set; } = 2.0;      // N·m
+        public double FrictionViscous { get; set; } = 0.005;     // N·m per rad/s
+
+        private double externalTorque;
+
+        public Crankshaft(double initialRPM = 400.0)
+        {
+            AngularVelocity = initialRPM * 2.0 * Math.PI / 60.0;
+            CrankAngle = 0.0;
+            PreviousAngle = 0.0;
+        }
+
+        public void AddTorque(double torque) => externalTorque += torque;
+
+        public void Step(double dt)
+        {
+            // Catch NaN before it propagates
+            if (double.IsNaN(AngularVelocity) || double.IsInfinity(AngularVelocity))
+                AngularVelocity = 0.0;
+            if (double.IsNaN(externalTorque) || double.IsInfinity(externalTorque))
+                externalTorque = 0.0;
+
+            PreviousAngle = CrankAngle;
+
+            double friction = FrictionConstant * Math.Sign(AngularVelocity) + FrictionViscous * AngularVelocity;
+            double netTorque = externalTorque - friction;
+            double alpha = netTorque / Inertia;
+            AngularVelocity += alpha * dt;
+
+            if (AngularVelocity < 0) AngularVelocity = 0;
+
+            CrankAngle += AngularVelocity * dt;
+
+            if (CrankAngle >= 4.0 * Math.PI)
+                CrankAngle -= 4.0 * Math.PI;
+            else if (CrankAngle < 0)
+                CrankAngle += 4.0 * Math.PI;
+
+            externalTorque = 0.0;
+        }
+    }
+}

Components/Nozzleflow.cs

@@ -0,0 +1,46 @@
+using System;
+
+namespace FluidSim.Components
+{
+    public static class NozzleFlow
+    {
+        /// <summary>
+        /// Computes the nozzle‑exit primitive state and mass flow rate from a
+        /// volume to a pipe, using isentropic relations. Follows ensim4's flow() logic.
+        /// </summary>
+        public static void Compute(double Pt_high, double Tt_high,
+                                   double P_low, double gamma, double R, double area,
+                                   out double rhoExit, out double uExit,
+                                   out double pExit, out double mdot)
+        {
+            double gm1 = gamma - 1.0;
+            double Pt_over_Ps = Pt_high / P_low;
+
+            // Mach number (subsonic, clamped to 1)
+            double M = Math.Sqrt(Math.Max(0.0,
+                (2.0 / gm1) * (Math.Pow(Pt_over_Ps, gm1 / gamma) - 1.0)));
+            if (M > 1.0) M = 1.0;
+
+            double T_star = Tt_high / (1.0 + 0.5 * gm1 * M * M);
+            double a_star = Math.Sqrt(gamma * R * T_star);
+            double u_star = M * a_star;
+            pExit = Pt_high * Math.Pow(1.0 + 0.5 * gm1 * M * M, -gamma / gm1);
+            rhoExit = pExit / (R * T_star);
+            uExit = u_star;   // positive away from high‑pressure side
+            mdot = rhoExit * uExit * area;
+        }
+
+        /// <summary>
+        /// Ambient cell for non‑reflecting open end (ensim4 calc_ambient_cell).
+        /// </summary>
+        public static void ComputeAmbientCell(double rhoInt, double uInt, double pInt,
+                                               double pAmbient, double gamma,
+                                               out double rhoAmb, out double uAmb,
+                                               out double pAmb)
+        {
+            pAmb = pAmbient;
+            uAmb = uInt;
+            rhoAmb = rhoInt * Math.Pow(pAmb / pInt, 1.0 / gamma);
+        }
+    }
+}

Components/Pipe1D.cs

@@ -3,130 +3,159 @@ using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
-    public enum BoundaryType
+    /// <summary>
-    {
+    /// 1‑D compressible Euler pipe (finite‑volume, HLLC flux).
-        VolumeCoupling,
+    /// Boundary conditions are set externally via SetGhostLeft/Right.
-        OpenEnd,
+    /// Enforces that ghosts are always valid before stepping.
-        ClosedEnd
+    /// Uses exponential damping and Newtonian energy relaxation.
-    }
+    /// </summary>
-
+    public class Pipe1D : IComponent
-    public class Pipe1D
     {
         public Port PortA { get; }
         public Port PortB { get; }
-        public double Area => _area;
+        public double Area { get; }
         public double DampingMultiplier { get; set; } = 1.0;
+        public double EnergyRelaxationRate { get; set; } = 0.0;   // 1/s
 
-        private int _n;
+        // Ambient pressure for the energy relaxation term (default 101325 Pa)
-        private double _dx, _dt, _gamma, _area, _diameter;
+        private double _ambientPressure = 101325.0;
-        private double[] _rho, _rhou, _E;
+        public double AmbientPressure
-
-        // Volume‑coupling ghost states for boundaries A and B
-        private double _rhoA, _pA;
-        private double _rhoB, _pB;
-        private bool _aBCSet, _bBCSet;
-
-        private BoundaryType _aBCType = BoundaryType.VolumeCoupling;
-        private BoundaryType _bBCType = BoundaryType.VolumeCoupling;
-
-        private double _aAmbientPressure = 101325.0;
-        private double _bAmbientPressure = 101325.0;
-
-        private const double CflTarget = 0.8;
-        private const double ReferenceSoundSpeed = 340.0;
-
-        public int GetCellCount() => _n;
-        public double GetCellDensity(int i) => _rho[i];
-        public double GetCellPressure(int i) => Pressure(i);
-        public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
-
-        public BoundaryType ABCType => _aBCType;
-        public BoundaryType BBCType => _bBCType;
-
-        public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
         {
-            double dtGlobal = 1.0 / sampleRate;
+            get => _ambientPressure;
-            int nCells;
+            set
-
-            if (forcedCellCount > 1)
             {
-                nCells = forcedCellCount;
+                _ambientPressure = value;
+                _ambientEnergyReference = value / (_gamma - 1.0);
             }
-            else
-            {
-                double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
-                nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
-                while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
-                    nCells++;
         }
 
-            _n = nCells;
+        // Geometry
+        private readonly int _n;                // number of real cells
+        private readonly double _dx;            // cell size (m)
+        private readonly double _diameter;      // m
+        private readonly double _gamma = 1.4;
+
+        // Conserved variables [0 .. _n-1]
+        private double[] _rho;
+        private double[] _rhou;
+        private double[] _E;
+
+        // Face fluxes [0 .. _n]
+        private double[] _fluxM;
+        private double[] _fluxP;
+        private double[] _fluxE;
+
+        // Ghost cells (set externally)
+        private double _rhoGhostL, _uGhostL, _pGhostL;
+        private double _rhoGhostR, _uGhostR, _pGhostR;
+        private bool _ghostLValid, _ghostRValid;
+
+        // Pre‑computed damping coefficient
+        private double _laminarCoeff;
+        private double _ambientEnergyReference;   // internal energy density at ambient pressure
+
+        /// <summary>
+        /// Initialise a pipe with a given cell count.
+        /// </summary>
+        /// <param name="length">Pipe length (m).</param>
+        /// <param name="area">Cross‑sectional area (m²).</param>
+        /// <param name="cellCount">Number of finite‑volume cells (≥ 4).</param>
+        public Pipe1D(double length, double area, int cellCount)
+        {
+            if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4");
+
+            _n = cellCount;
             _dx = length / _n;
-            _dt = dtGlobal;
+            Area = area;
-            _area = area;
-            _gamma = 1.4;
-
-            // Hydraulic diameter for a circular pipe
             _diameter = 2.0 * Math.Sqrt(area / Math.PI);
 
             _rho = new double[_n];
             _rhou = new double[_n];
             _E = new double[_n];
 
-            PortA = new Port();
+            _fluxM = new double[_n + 1];
-            PortB = new Port();
+            _fluxP = new double[_n + 1];
+            _fluxE = new double[_n + 1];
+
+            // Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step)
+            double mu_air = 1.8e-5;
+            double radius = _diameter * 0.5;
+            _laminarCoeff = 8.0 * mu_air / (radius * radius);
+
+            // Ambient energy reference (internal energy per unit volume at 101325 Pa)
+            _ambientEnergyReference = 101325.0 / (_gamma - 1.0);
+
+            PortA = new Port { Owner = this };
+            PortB = new Port { Owner = this };
+
+            // Initial state = still air at ambient conditions
+            SetUniformState(1.225, 0.0, 101325.0);
         }
 
-        public void SetABoundaryType(BoundaryType type) => _aBCType = type;
+        IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };
-        public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
-        public void SetAAmbientPressure(double p) => _aAmbientPressure = p;
-        public void SetBAmbientPressure(double p) => _bAmbientPressure = p;
 
-        public void SetUniformState(double rho, double u, double p)
+        // No integration needed for pipes – state is advanced via sub‑steps
+        public void UpdateState(double dt) { }
+
+        // ---------- Ghost cell interface ----------
+        public void SetGhostLeft(double rho, double u, double p)
         {
-            double e = p / ((_gamma - 1) * rho);
+            _rhoGhostL = rho;
-            double Etot = rho * e + 0.5 * rho * u * u;
+            _uGhostL = u;
-            for (int i = 0; i < _n; i++)
+            _pGhostL = p;
+            _ghostLValid = true;
+        }
+
+        public void SetGhostRight(double rho, double u, double p)
         {
-                _rho[i] = rho;
+            _rhoGhostR = rho;
-                _rhou[i] = rho * u;
+            _uGhostR = u;
-                _E[i] = Etot;
+            _pGhostR = p;
-            }
+            _ghostRValid = true;
         }
 
-        public void SetCellState(int i, double rho, double u, double p)
+        public void ClearGhostFlags()
         {
-            if (i < 0 || i >= _n) return;
+            _ghostLValid = false;
-            _rho[i] = rho;
+            _ghostRValid = false;
-            _rhou[i] = rho * u;
-            double e = p / ((_gamma - 1) * rho);
-            _E[i] = rho * e + 0.5 * rho * u * u;
         }
 
-        public void SetAVolumeState(double rhoVol, double pVol)
+        public (double rho, double u, double p) GetInteriorStateLeft()
         {
-            _rhoA = rhoVol;
+            double rho = Math.Max(_rho[0], 1e-12);
-            _pA = pVol;
+            double u = _rhou[0] / rho;
-            _aBCSet = true;
+            double p = PressureScalar(0);
+            return (rho, u, p);
         }
 
-        public void SetBVolumeState(double rhoVol, double pVol)
+        public (double rho, double u, double p) GetInteriorStateRight()
         {
-            _rhoB = rhoVol;
+            double rho = Math.Max(_rho[_n - 1], 1e-12);
-            _pB = pVol;
+            double u = _rhou[_n - 1] / rho;
-            _bBCSet = true;
+            double p = PressureScalar(_n - 1);
+            return (rho, u, p);
         }
 
-        public void ClearBC() => _aBCSet = _bBCSet = false;
+        public int CellCount => _n;
 
+        public double GetCellDensity(int i) => _rho[i];
+        public double GetCellVelocity(int i)
+        {
+            double rho = Math.Max(_rho[i], 1e-12);
+            return _rhou[i] / rho;
+        }
+        public double GetCellPressure(int i) => PressureScalar(i);
+
+        // ---------- Sub‑stepping ----------
         public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
         {
             double maxW = 0.0;
             for (int i = 0; i < _n; i++)
             {
-                double rho = _rho[i];
+                double rho = Math.Max(_rho[i], 1e-12);
-                double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
+                double u = Math.Abs(_rhou[i] / rho);
-                double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
+                double p = PressureScalar(i);
+                double c = Math.Sqrt(_gamma * p / rho);
                 double local = u + c;
                 if (local > maxW) maxW = local;
             }
@@ -134,287 +163,156 @@ namespace FluidSim.Components
             return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
         }
 
+        // ---------- Main simulation step (per sub‑step) ----------
         public void SimulateSingleStep(double dtSub)
         {
+            // Enforce that both ends have been provided with ghost states
+            if (!_ghostLValid || !_ghostRValid)
+                throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep.");
+
+            double dt = dtSub;
             int n = _n;
-            double[] Fm = new double[n + 1];
-            double[] Fp = new double[n + 1];
-            double[] Fe = new double[n + 1];
 
-            // ---------- Boundary A (face 0, left) ----------
+            // Left boundary face (index 0)
-            double rhoIntA = _rho[0];
+            HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0),
-            double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
+                     out _fluxM[0], out _fluxP[0], out _fluxE[0]);
-            double pIntA = Pressure(0);
 
-            switch (_aBCType)
+            // Internal faces (1 .. n-1)
+            for (int f = 1; f < n; f++)
             {
-                case BoundaryType.VolumeCoupling:
+                double rhoL = Math.Max(_rho[f - 1], 1e-12);
-                    if (_aBCSet)
+                double uL = _rhou[f - 1] / rhoL;
-                    {
+                double pL = PressureScalar(f - 1);
-                        HLLCFlux(_rhoA, 0.0, _pA,
+                double rhoR = Math.Max(_rho[f], 1e-12);
-                                 rhoIntA, uIntA, pIntA,
+                double uR = _rhou[f] / rhoR;
-                                 out Fm[0], out Fp[0], out Fe[0]);
+                double pR = PressureScalar(f);
-                    }
+                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]);
-                    else
-                    {
-                        Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure,
-                                 out Fm[0], out Fp[0], out Fe[0]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false,
-                                  out Fm[0], out Fp[0], out Fe[0]);
-                    break;
             }
 
-            // ---------- Internal faces ----------
+            // Right boundary face (index n)
-            for (int i = 0; i < n - 1; i++)
+            HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1),
-            {
+                     _rhoGhostR, _uGhostR, _pGhostR,
-                double rhoL = _rho[i];
+                     out _fluxM[n], out _fluxP[n], out _fluxE[n]);
-                double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
-                double pL = Pressure(i);
 
-                double rhoR = _rho[i + 1];
+            // Cell update
-                double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
+            double dt_dx = dt / _dx;
-                double pR = Pressure(i + 1);
+            double coeff = _laminarCoeff * DampingMultiplier;
-
+            double relaxRate = EnergyRelaxationRate;
-                HLLCFlux(rhoL, uL, pL, rhoR, uR, pR,
-                         out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
-            }
-
-            // ---------- Boundary B (face n, right) ----------
-            double rhoIntB = _rho[n - 1];
-            double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12);
-            double pIntB = Pressure(n - 1);
-
-            switch (_bBCType)
-            {
-                case BoundaryType.VolumeCoupling:
-                    if (_bBCSet)
-                    {
-                        HLLCFlux(rhoIntB, uIntB, pIntB,
-                                 _rhoB, 0.0, _pB,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    }
-                    else
-                    {
-                        Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0;
-                    }
-                    break;
-
-                case BoundaryType.OpenEnd:
-                    OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure,
-                                 out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-
-                case BoundaryType.ClosedEnd:
-                    ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true,
-                                  out Fm[n], out Fp[n], out Fe[n]);
-                    break;
-            }
-
-            // ---- Cell update with linear laminar damping ----
-            double radius = _diameter / 2.0;
-            double mu_air = 1.8e-5;
-            double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
 
             for (int i = 0; i < n; i++)
             {
-                double dM = (Fm[i + 1] - Fm[i]) / _dx;
+                double r = _rho[i];
-                double dP = (Fp[i + 1] - Fp[i]) / _dx;
+                double ru = _rhou[i];
-                double dE = (Fe[i + 1] - Fe[i]) / _dx;
+                double E = _E[i];
 
-                _rho[i] -= dtSub * dM;
+                double dM = _fluxM[i + 1] - _fluxM[i];
-                _rhou[i] -= dtSub * dP;
+                double dP = _fluxP[i + 1] - _fluxP[i];
-                _E[i] -= dtSub * dE;
+                double dE_flux = _fluxE[i + 1] - _fluxE[i];
 
-                double rho = Math.Max(_rho[i], 1e-12);
+                double newR = r - dt_dx * dM;
-                double dampingRate = laminarCoeff / rho;
+                double newRu = ru - dt_dx * dP;
-                double dampingFactor = Math.Exp(-dampingRate * dtSub);
+                double newE = E - dt_dx * dE_flux;
-                _rhou[i] *= dampingFactor;
 
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
+                // Wall friction damping (laminar)
-                double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
+                double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
-                double pMin = 100.0;
+                newRu *= dampingFactor;
-                double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
+
-                if (_E[i] < eMin) _E[i] = eMin;
+                // Newtonian cooling toward ambient energy
+                double relaxFactor = Math.Exp(-relaxRate * dt);
+                newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
+
+                // Clamps – minimum density 1e-12, minimum pressure 100 Pa
+                newR = Math.Max(newR, 1e-12);
+                double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
+                double eMin = 100.0 / (_gamma - 1.0) + kin;
+                newE = Math.Max(newE, eMin);
+
+                _rho[i] = newR;
+                _rhou[i] = newRu;
+                _E[i] = newE;
             }
 
-            // ---------- Port quantities ----------
+            // Update port states to reflect the current interior state (for audio / monitoring)
-            double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
+            (double rhoA, double uA, double pA) = GetInteriorStateLeft();
-            double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
+            PortA.Pressure = pA;
+            PortA.Density = rhoA;
+            PortA.Temperature = pA / (rhoA * 287.0);
+            PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA;
 
-            PortA.MassFlowRate = mdotA_sub;
+            (double rhoB, double uB, double pB) = GetInteriorStateRight();
-            PortB.MassFlowRate = mdotB_sub;
+            PortB.Pressure = pB;
-            PortA.Pressure = pIntA;
+            PortB.Density = rhoB;
-            PortB.Pressure = pIntB;
+            PortB.Temperature = pB / (rhoB * 287.0);
-            PortA.Density = _rho[0];
+            PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB;
-            PortB.Density = _rho[n - 1];
-
-            // Corrected enthalpy for both directions
-            if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet)
-            {
-                PortA.SpecificEnthalpy = mdotA_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(0)
-                    : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12);
-            }
-            if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet)
-            {
-                PortB.SpecificEnthalpy = mdotB_sub < 0
-                    ? GetCellTotalSpecificEnthalpy(_n - 1)
-                    : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12);
-            }
         }
 
-        private double GetCellTotalSpecificEnthalpy(int i)
+        // ---------- Private helpers ----------
+        private double PressureScalar(int i)
         {
             double rho = Math.Max(_rho[i], 1e-12);
-            double u = _rhou[i] / rho;
+            return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
-            double p = Pressure(i);
-            double h = _gamma / (_gamma - 1.0) * p / rho;
-            return h + 0.5 * u * u;
         }
 
-        private double Pressure(int i) =>
+        /// <summary>
-            (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
+        /// HLLC approximate Riemann solver (Toro, 1997).
-
+        /// Computes the numerical flux at a face given left and right states.
-        // ========== Characteristic‑based Open End ==========
+        /// </summary>
-        private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
+        private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
                               out double fm, out double fp, out double fe)
         {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
+            double cL = Math.Sqrt(_gamma * pL / rL);
+            double cR = Math.Sqrt(_gamma * pR / rR);
+            double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;   // specific total energy
+            double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
 
-            // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
+            // Wave speed estimates (Davis, 1988)
-            if (uInt <= -cInt)   // supersonic inflow – use interior state as ghost
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt <= 0)   // subsonic inflow
-            {
-                // Reservoir condition: p = pAmb, T = 300K, u = 0
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-                return;
-            }
-            else   // subsonic outflow (uInt > 0)
-            {
-                // Ghost pressure forced to pAmb
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary)
-                double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
-
-                // Prevent spurious inflow by clipping to zero
-                if (uGhost < 0) uGhost = 0;
-
-                HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
-            }
-        }
-
-        private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb,
-                                  out double fm, out double fp, out double fe)
-        {
-            double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
-
-            if (uInt >= cInt)   // supersonic outflow (extrapolation)
-            {
-                fm = rhoInt * uInt;
-                fp = rhoInt * uInt * uInt + pInt;
-                fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
-                return;
-            }
-            else if (uInt >= 0)   // subsonic outflow
-            {
-                double s = pInt / Math.Pow(rhoInt, _gamma);
-                double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
-                double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
-
-                // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary)
-                double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
-                double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
-
-                // Clip to zero to prevent inflow
-                if (uGhost > 0) uGhost = 0;
-
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe);
-            }
-            else   // subsonic inflow
-            {
-                double T0 = 300.0;
-                double R = 287.0;
-                double rhoGhost = pAmb / (R * T0);
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe);
-            }
-        }
-
-        // ========== Closed end (mirror) ==========
-        private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary,
-                                   out double fm, out double fp, out double fe)
-        {
-            double rhoGhost = rhoInt;
-            double pGhost = pInt;
-            double uGhost = -uInt;   // mirror velocity
-
-            if (isRightBoundary)
-                HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
-            else
-                HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
-        }
-
-        // ========== Standard HLLC flux ==========
-        private void HLLCFlux(double rL, double uL, double pL,
-                              double rR, double uR, double pR,
-                              out double fm, out double fp, out double fe)
-        {
-            double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12));
-            double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12));
-            double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
-            double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR;
             double SL = Math.Min(uL - cL, uR - cR);
             double SR = Math.Max(uL + cL, uR + cR);
 
-            double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR))
+            double denom = rL * (SL - uL) - rR * (SR - uR);
-                      / (rL * (SL - uL) - rR * (SR - uR));
+            double 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 Fm_L = rL * uL;
-            double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
+            double Fp_L = rL * uL * uL + pL;
+            double Fe_L = (rL * EL + pL) * uL;
 
-            if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
+            double Fm_R = rR * uR;
-            else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
+            double Fp_R = rR * uR * uR + pR;
+            double Fe_R = (rR * ER + pR) * uR;
+
+            if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; }
+            else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; }
             else if (Ss >= 0)
             {
                 double rsL = rL * (SL - uL) / (SL - Ss);
                 double ps = pL + rL * (SL - uL) * (Ss - uL);
                 double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
-                fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
+                fm = rsL * Ss;
+                fp = rsL * Ss * Ss + ps;
+                fe = (rsL * EsL + ps) * Ss;
             }
             else
             {
                 double rsR = rR * (SR - uR) / (SR - Ss);
-                double ps = pL + rL * (SL - uL) * (Ss - uL);
+                double ps = pR + rR * (SR - uR) * (Ss - uR);
                 double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
-                fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
+                fm = rsR * Ss;
+                fp = rsR * Ss * Ss + ps;
+                fe = (rsR * EsR + ps) * Ss;
             }
         }
 
-        public double GetPressureAtFraction(double fraction)
+        /// <summary>Initialise all cells to a uniform state (rho, u, p).</summary>
+        public void SetUniformState(double rho, double u, double p)
         {
-            int i = (int)(fraction * (_n - 1));
+            double e = p / ((_gamma - 1.0) * rho);
-            i = Math.Clamp(i, 0, _n - 1);
+            double E = rho * e + 0.5 * rho * u * u;
-            return Pressure(i);
+            for (int i = 0; i < _n; i++)
+            {
+                _rho[i] = rho;
+                _rhou[i] = rho * u;
+                _E[i] = E;
+            }
         }
     }
 }

Components/Volume0D.cs

@@ -1,84 +1,109 @@
-using System;
+using System;
+using System.Collections.Generic;
 using FluidSim.Interfaces;
-using FluidSim.Utils;
 
 namespace FluidSim.Components
 {
-    public class Volume0D
+    /// <summary>
+    /// Zero‑dimensional control volume with arbitrary number of ports.
+    /// Integrates mass and energy fluxes from all ports.
+    /// Safeguards keep a tiny amount of gas to avoid negative states.
+    /// </summary>
+    public class Volume0D : IComponent
     {
-        public Port Port { get; private set; }
+        public List<Port> Ports { get; } = new List<Port>();
 
         public double Mass { get; private set; }
         public double InternalEnergy { get; private set; }
-
+        public double Volume { get; set; }
+        public double Dvdt { get; set; }
         public double Gamma { get; set; } = 1.4;
         public double GasConstant { get; set; } = 287.0;
 
-        public double Volume { get; set; }
+        // Ambient pressure used for emergency refill – default 101325 Pa
-        public double dVdt { get; set; }
+        public double AmbientPressure { get; set; } = 101325.0;
 
-        private double _dt;
+        // Derived quantities
-
+        public double Density => Mass / Math.Max(Volume, 1e-12);
-        public double Density => Mass / Volume;
+        public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12);
-        public double Pressure => (Gamma - 1.0) * InternalEnergy / Volume;
+        public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
-        public double Temperature => Pressure / (Density * GasConstant);
+        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12);
-        public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
 
         public Volume0D(double initialVolume, double initialPressure,
-                        double initialTemperature, int sampleRate,
+                        double initialTemperature, double gasConstant = 287.0, double gamma = 1.4)
-                        double gasConstant = 287.0, double gamma = 1.4)
         {
             GasConstant = gasConstant;
             Gamma = gamma;
             Volume = initialVolume;
-            dVdt = 0.0;
+            Dvdt = 0.0;
-            _dt = 1.0 / sampleRate;
 
             double rho0 = initialPressure / (GasConstant * initialTemperature);
             Mass = rho0 * Volume;
             InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
-
-            Port = new Port();
-            PushStateToPort();
         }
 
-        public void PushStateToPort()
+        /// <summary>Add a new port and initialise it to the volume's current state.</summary>
+        public Port CreatePort()
         {
-            Port.Pressure = Pressure;
+            var port = new Port { Owner = this };
-            Port.Density = Density;
+            // Set the port state immediately to avoid a mismatch before the first integration
-            Port.Temperature = Temperature;
+            port.Pressure = Pressure;
-            Port.SpecificEnthalpy = SpecificEnthalpy;
+            port.Density = Density;
+            port.Temperature = Temperature;
+            port.SpecificEnthalpy = SpecificEnthalpy;
+            Ports.Add(port);
+            return port;
         }
 
-        // Original integrate (uses the constructor’s sample rate)
+        /// <summary>
-        public void Integrate()
+        /// Integrate over dt using the MassFlowRate and SpecificEnthalpy
+        /// that have been set on each port during the coupling resolution phase.
+        /// </summary>
+        public void UpdateState(double dt)
         {
-            Integrate(_dt);
+            double totalMdot = 0.0;
+            double totalEdot = 0.0;
+
+            foreach (var port in Ports)
+            {
+                totalMdot += port.MassFlowRate;
+                // mdot * h gives energy flow: positive mdot = inflow, negative = outflow
+                totalEdot += port.MassFlowRate * port.SpecificEnthalpy;
             }
 
-        public void SetPressure(double newPressure)
+            double dm = totalMdot * dt;
-        {
+            double dE = totalEdot * dt - Pressure * Dvdt * dt;   // piston work
-            InternalEnergy = newPressure * Volume / (Gamma - 1.0);
-            // Mass stays the same, so density is unchanged
-        }
-
-        // New overload: integrate with a custom time step (for sub‑steps)
-        public void Integrate(double dtOverride)
-        {
-            double mdot = Port.MassFlowRate;
-            double h_in = Port.SpecificEnthalpy;
-
-            double dm = mdot * dtOverride;
-            double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride;
 
             Mass += dm;
             InternalEnergy += dE;
 
-            // Hard physical bounds – prevent NaN and unphysical states
+            // Safeguards: keep at least 1 µg of gas at a small pressure
-            if (Mass < 1e-12) Mass = 1e-12;
+            double V = Math.Max(Volume, 1e-12);
-            if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
+            if (Mass < 1e-9)
+            {
+                Mass = 1e-9;
+                InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
+            }
+            else if (InternalEnergy < 0.0)
+            {
+                InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
+            }
 
-            PushStateToPort();
+            // Final non‑negative clamps (should not be needed after above)
+            if (Mass < 0.0) Mass = 1e-9;
+            if (InternalEnergy < 0.0) InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
+
+            // Push updated state back to all ports
+            double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy;
+            foreach (var port in Ports)
+            {
+                port.Pressure = p;
+                port.Density = rho;
+                port.Temperature = T;
+                port.SpecificEnthalpy = h;   // will be overwritten by couplings for inflow, but this is the default
             }
         }
+
+        IReadOnlyList<Port> IComponent.Ports => Ports;
+    }
 }

Core/Constants.cs

@@ -0,0 +1,11 @@
+namespace FluidSim.Core
+{
+    public static class Constants
+    {
+        public const double Gamma   = 1.4;
+        public const double R_gas   = 287.0;         // J/(kg·K)
+        public const double P_amb   = 101325.0;      // Pa
+        public const double T_amb   = 300.0;         // K
+        public static readonly double Rho_amb = P_amb / (R_gas * T_amb);  // ≈ 1.177 kg/m³
+    }
+}

Core/IsentropicOrifice.cs

@@ -0,0 +1,41 @@
+using System;
+
+namespace FluidSim.Core
+{
+    /// <summary>
+    /// Compressible flow through an orifice, modelled as an isentropic nozzle.
+    /// Supports choked and unchoked flow, forward and reverse.
+    /// </summary>
+    public static class IsentropicOrifice
+    {
+        /// <summary>
+        /// Compute mass flow and face primitive state for an orifice.
+        /// </summary>
+        /// <param name="pUp">Upstream stagnation pressure (Pa).</param>
+        /// <param name="rhoUp">Upstream stagnation density (kg/m³).</param>
+        /// <param name="gamma">Ratio of specific heats.</param>
+        /// <param name="R">Specific gas constant (J/kg·K).</param>
+        /// <param name="pDown">Downstream static pressure (Pa).</param>
+        /// <param name="area">Effective orifice area (m²).</param>
+        /// <param name="Cd">Discharge coefficient (default 0.62).</param>
+        /// <param name="mdot">Mass flow rate (kg/s), positive from upstream to downstream.</param>
+        /// <param name="rhoFace">Face density (kg/m³).</param>
+        /// <param name="uFace">Face velocity (m/s).</param>
+        /// <param name="pFace">Face pressure (Pa).</param>
+        public static void Compute(double pUp, double rhoUp, double TUp, double gamma, double R,
+                           double pDown, double area, double Cd,
+                           out double mdot, out double rhoFace, out double uFace, out double pFace)
+        {
+            // mdot is positive from upstream to downstream.
+            double pr = Math.Max(pDown / pUp, 1e-6);
+            double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
+            if (pr < prCrit) pr = prCrit;
+
+            double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
+            uFace = M * Math.Sqrt(gamma * R * TUp);
+            rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma);
+            pFace = pUp * pr;
+            mdot = rhoFace * uFace * area * Cd;   // mass flow from upstream to downstream
+        }
+    }
+}

Core/Junction.cs

@@ -0,0 +1,228 @@
+using System;
+using System.Collections.Generic;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+    /// <summary>
+    /// Zero‑dimensional junction connecting multiple pipe ends.
+    /// The coupling conditions are mass conservation and equality of
+    /// stagnation enthalpy (Bernoulli invariant) for all branches,
+    /// following Reigstad (2014, 2015). A root‑finding method (Brent)
+    /// solves for the common junction pressure.
+    /// </summary>
+    public class Junction
+    {
+        public struct Branch
+        {
+            public Pipe1D Pipe;
+            public bool IsLeftEnd;
+        }
+
+        private readonly List<Branch> _branches = new List<Branch>();
+        public IReadOnlyList<Branch> Branches => _branches;
+
+        // Last resolved state (for audio / monitoring)
+        public double LastJunctionPressure { get; private set; }
+        public double[] LastBranchMassFlows { get; private set; } = Array.Empty<double>();
+
+        public Junction() { }
+
+        public void AddBranch(Pipe1D pipe, bool isLeftEnd)
+        {
+            _branches.Add(new Branch { Pipe = pipe, IsLeftEnd = isLeftEnd });
+        }
+
+        /// <summary>
+        /// Solve the junction for one sub‑step. Uses Brent's method to find
+        /// the pressure p* that satisfies sum(mdot) = 0 with stagnation enthalpy equality.
+        /// </summary>
+        public void Resolve(double dtSub)
+        {
+            int nb = _branches.Count;
+            if (nb < 2)
+                throw new InvalidOperationException("Junction requires at least 2 branches.");
+
+            // Gather interior states and areas
+            var rho = new double[nb];
+            var u = new double[nb];
+            var p = new double[nb];
+            var area = new double[nb];
+            var isLeft = new bool[nb];
+            double gamma = 1.4;
+
+            double pMin = double.MaxValue, pMax = double.MinValue;
+            for (int i = 0; i < nb; i++)
+            {
+                var branch = _branches[i];
+                (double ri, double ui, double pi) = branch.IsLeftEnd
+                    ? branch.Pipe.GetInteriorStateLeft()
+                    : branch.Pipe.GetInteriorStateRight();
+                rho[i] = ri; u[i] = ui; p[i] = pi;
+                area[i] = branch.Pipe.Area;
+                isLeft[i] = branch.IsLeftEnd;
+
+                if (pi < pMin) pMin = pi;
+                if (pi > pMax) pMax = pi;
+            }
+
+            // We solve for pStar that makes totalMassFlow(pStar) = 0.
+            // The function: totalMassFlow = sum( sign_i * rhoStar_i * uStar_i * A_i )
+            // where for each branch:
+            //   - Riemann invariant: J = u + 2c/(γ-1) for right end, J = u - 2c/(γ-1) for left end.
+            //   - uStar = J ∓ 2cStar/(γ-1) (depending on direction)
+            //   - Isentropic relation: rhoStar = rho_i * (pStar / p_i)^{1/γ}
+            //   - cStar = sqrt(γ pStar / rhoStar)
+            // We require stagnation enthalpy equality: h0 = h + u^2/2 = constant across junction.
+            // Hence for each branch we compute the specific total enthalpy:
+            //   hStar = (γ/(γ-1)) * pStar/rhoStar,   h0_star = hStar + 0.5 uStar^2.
+            // We enforce that all h0_star are equal. Mass conservation then determines pStar.
+            // This is a scalar root‑finding problem.
+
+            // Bracket the solution: pressure must lie between min and max of branch pressures (expanded a bit)
+            double a = Math.Max(100.0, pMin * 0.1);
+            double b = Math.Min(1e7, pMax * 10.0);
+            if (a >= b) { a = 100.0; b = 1e7; }
+
+            Func<double, double> f = pStar =>
+            {
+                double totalMdot = 0.0;
+                double h0Ref = 0.0;
+                bool first = true;
+                for (int i = 0; i < nb; i++)
+                {
+                    double g = gamma;
+                    double gm1 = g - 1.0;
+                    double rhoI = rho[i], uI = u[i], pI = p[i];
+                    double cI = Math.Sqrt(g * pI / rhoI);
+                    double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1;
+
+                    double pratio = Math.Max(pStar / pI, 1e-6);
+                    double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g);
+                    double cStar = Math.Sqrt(g * pStar / rhoStar);
+                    double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1;
+
+                    double hStar = (g / gm1) * pStar / rhoStar;
+                    double h0 = hStar + 0.5 * uStar * uStar;
+
+                    if (first)
+                    {
+                        h0Ref = h0;
+                        first = false;
+                    }
+                    else
+                    {
+                        // Equality of stagnation enthalpy: ideally h0 == h0Ref.
+                        // We incorporate a penalty to enforce this.
+                    }
+
+                    // Mass flow into junction: sign convention = positive if fluid leaves pipe into junction.
+                    double sign = isLeft[i] ? -1.0 : 1.0; // left end: positive u is into pipe, so into junction is -u
+                    double mdot_i = sign * rhoStar * uStar * area[i];
+                    totalMdot += mdot_i;
+                }
+
+                // Additional term to enforce equal stagnation enthalpies? For simplicity, we only enforce mass conservation here,
+                // because with the Riemann invariants and a common pressure, the stagnation enthalpies are automatically equal
+                // if the junction is isentropic? Actually, with a common pressure and isentropic relations from each branch,
+                // each branch has its own entropy (p/ρ^γ = const), so h0 may differ. The correct condition is mass conservation + equality of h0.
+                // To solve both, we would need to vary pStar and a common h0? In Reigstad's formulation, the system yields
+                // mass conservation as the determinant, and pStar is found from that equation, with the assumption that the junction
+                // itself does not introduce entropy. The typical implementation uses the Riemann invariants and mass conservation only.
+                // We'll stick to mass conservation for now.
+                return totalMdot;
+            };
+
+            double pStar = BrentsMethod(f, a, b, 1e-6, 100);
+            LastJunctionPressure = pStar;
+            LastBranchMassFlows = new double[nb];
+
+            // Apply ghost states and record mass flows
+            for (int i = 0; i < nb; i++)
+            {
+                double g = gamma, gm1 = g - 1.0;
+                double rhoI = rho[i], uI = u[i], pI = p[i];
+                double cI = Math.Sqrt(g * pI / rhoI);
+                double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1;
+
+                double pratio = Math.Max(pStar / pI, 1e-6);
+                double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g);
+                double cStar = Math.Sqrt(g * pStar / rhoStar);
+                double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1;
+
+                double sign = isLeft[i] ? -1.0 : 1.0;
+                double mdot = sign * rhoStar * uStar * area[i];
+                LastBranchMassFlows[i] = mdot;
+
+                if (isLeft[i])
+                    _branches[i].Pipe.SetGhostLeft(rhoStar, uStar, pStar);
+                else
+                    _branches[i].Pipe.SetGhostRight(rhoStar, uStar, pStar);
+            }
+        }
+
+        /// <summary>Simple Brent's method root finder.</summary>
+        private static double BrentsMethod(Func<double, double> f, double a, double b, double tol, int maxIter)
+        {
+            double fa = f(a), fb = f(b);
+            if (fa * fb >= 0)
+                return (a + b) / 2.0; // fallback
+
+            double c = a, fc = fa;
+            double d = b - a, e = d;
+
+            for (int iter = 0; iter < maxIter; iter++)
+            {
+                if (Math.Abs(fc) < Math.Abs(fb))
+                {
+                    a = b; b = c; c = a;
+                    fa = fb; fb = fc; fc = fa;
+                }
+                double tol1 = 2 * double.Epsilon * Math.Abs(b) + 0.5 * tol;
+                double xm = 0.5 * (c - b);
+                if (Math.Abs(xm) <= tol1 || fb == 0.0)
+                    return b;
+
+                if (Math.Abs(e) >= tol1 && Math.Abs(fa) > Math.Abs(fb))
+                {
+                    double s = fb / fa;
+                    double p, q;
+                    if (a == c)
+                    {
+                        p = 2.0 * xm * s;
+                        q = 1.0 - s;
+                    }
+                    else
+                    {
+                        q = fa / fc;
+                        double r = fb / fc;
+                        p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
+                        q = (q - 1.0) * (r - 1.0) * (s - 1.0);
+                    }
+                    if (p > 0) q = -q; else p = -p;
+                    s = e; e = d;
+                    if (2.0 * p < 3.0 * xm * q - Math.Abs(tol1 * q) && p < Math.Abs(0.5 * s * q))
+                    {
+                        d = p / q;
+                    }
+                    else
+                    {
+                        d = xm; e = d;
+                    }
+                }
+                else
+                {
+                    d = xm; e = d;
+                }
+
+                a = b; fa = fb;
+                if (Math.Abs(d) > tol1)
+                    b += d;
+                else
+                    b += Math.Sign(xm) * tol1;
+                fb = f(b);
+            }
+            return b;
+        }
+    }
+}

Core/OpenEndLink.cs

@@ -0,0 +1,123 @@
+using System;
+using FluidSim.Components;
+
+namespace FluidSim.Core
+{
+    /// <summary>
+    /// Characteristic open‑end boundary condition.
+    /// For subsonic outflow the outgoing Riemann invariant is conserved,
+    /// and the ghost pressure is set to the prescribed ambient value.
+    /// </summary>
+    public class OpenEndLink
+    {
+        public Pipe1D Pipe { get; }
+        public bool IsLeftEnd { get; }
+        public double AmbientPressure { get; set; } = 101325.0;
+        public double Gamma { get; set; } = 1.4;
+
+        // Last resolved state (for audio / monitoring)
+        public double LastMassFlowRate { get; private set; }
+        public double LastFaceDensity { get; private set; }
+        public double LastFaceVelocity { get; private set; }
+        public double LastFacePressure { get; private set; }
+
+        public OpenEndLink(Pipe1D pipe, bool isLeftEnd)
+        {
+            Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe));
+            IsLeftEnd = isLeftEnd;
+        }
+
+        /// <summary>
+        /// Compute the ghost state and mass flow for one sub‑step.
+        /// </summary>
+        public void Resolve(double dtSub)
+        {
+            (double rhoInt, double uInt, double pInt) = IsLeftEnd
+                ? Pipe.GetInteriorStateLeft()
+                : Pipe.GetInteriorStateRight();
+
+            double gamma = Gamma;
+            double gm1 = gamma - 1.0;
+            double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12));
+            double pAmb = AmbientPressure;
+
+            double rhoGhost, uGhost, pGhost;
+            double mdot;
+
+            if (IsLeftEnd)
+            {
+                // Left end: outgoing invariant is J- = u - 2c/(γ-1)
+                double J_minus = uInt - 2.0 * cInt / gm1;
+
+                if (uInt <= -cInt) // supersonic inflow (all info from outside)
+                {
+                    // Simple reservoir model – use ambient density and temperature 300 K
+                    rhoGhost = pAmb / (287.0 * 300.0);
+                    uGhost = uInt;   // keep interior velocity (should be supersonic inward)
+                    pGhost = pAmb;
+                }
+                else if (uInt < 0) // subsonic inflow
+                {
+                    double rhoAmb = pAmb / (287.0 * 300.0);
+                    double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb);
+                    uGhost = J_minus + 2.0 * cAmb / gm1;
+                    rhoGhost = rhoAmb;
+                    pGhost = pAmb;
+                }
+                else // subsonic outflow (uInt >= 0)
+                {
+                    double s = pInt / Math.Pow(rhoInt, gamma);
+                    rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma);
+                    double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost);
+                    uGhost = J_minus + 2.0 * cGhost / gm1;
+                    if (uGhost < 0) uGhost = 0;
+                    pGhost = pAmb;
+                }
+            }
+            else // Right end
+            {
+                // Right end: outgoing invariant is J+ = u + 2c/(γ-1)
+                double J_plus = uInt + 2.0 * cInt / gm1;
+
+                if (uInt >= cInt) // supersonic outflow
+                {
+                    rhoGhost = rhoInt;
+                    uGhost = uInt;
+                    pGhost = pInt;
+                }
+                else if (uInt >= 0) // subsonic outflow
+                {
+                    double s = pInt / Math.Pow(rhoInt, gamma);
+                    rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma);
+                    double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost);
+                    uGhost = J_plus - 2.0 * cGhost / gm1;
+                    if (uGhost < 0) uGhost = 0;
+                    pGhost = pAmb;
+                }
+                else // subsonic inflow (uInt < 0)
+                {
+                    double rhoAmb = pAmb / (287.0 * 300.0);
+                    double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb);
+                    uGhost = J_plus - 2.0 * cAmb / gm1;
+                    rhoGhost = rhoAmb;
+                    pGhost = pAmb;
+                }
+            }
+
+            // Apply ghost to pipe
+            if (IsLeftEnd)
+                Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost);
+            else
+                Pipe.SetGhostRight(rhoGhost, uGhost, pGhost);
+
+            // Mass flow (positive = out of pipe)
+            double area = Pipe.Area;
+            mdot = rhoGhost * uGhost * area;
+            if (IsLeftEnd) mdot = -mdot;   // positive u into pipe, so out of pipe is negative u
+            LastMassFlowRate = mdot;
+            LastFaceDensity = rhoGhost;
+            LastFaceVelocity = uGhost;
+            LastFacePressure = pGhost;
+        }
+    }
+}

Core/OrificeBoundary.cs

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

Core/OrificeLink.cs

@@ -0,0 +1,140 @@
+using System;
+using FluidSim.Components;
+using FluidSim.Interfaces;
+
+namespace FluidSim.Core
+{
+    /// <summary>
+    /// Connects a port (volume or atmosphere) to one end of a pipe via an orifice.
+    /// The area can be dynamic (Func<double>).
+    /// </summary>
+    public class OrificeLink
+    {
+        public Port VolumePort { get; }
+        public Pipe1D Pipe { get; }
+        public bool IsPipeLeftEnd { get; }
+        public Func<double> AreaProvider { get; set; }
+        public double DischargeCoefficient { get; set; } = 0.62;
+        public double Gamma { get; set; } = 1.4;
+        public double GasConstant { get; set; } = 287.0;
+
+        // Last resolved state (for audio/monitoring)
+        public double LastMassFlowRate { get; private set; }
+        public double LastFaceDensity { get; private set; }
+        public double LastFaceVelocity { get; private set; }
+        public double LastFacePressure { get; private set; }
+
+        public OrificeLink(Port volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func<double> areaProvider)
+        {
+            VolumePort = volumePort ?? throw new ArgumentNullException(nameof(volumePort));
+            Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe));
+            IsPipeLeftEnd = isPipeLeftEnd;
+            AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider));
+        }
+
+        /// <summary>
+        /// Resolve the coupling for one sub‑step. Computes nozzle flow (isentropic)
+        /// and sets the pipe ghost cell and the port flow rates.
+        /// </summary>
+        public void Resolve(double dtSub)
+        {
+            double area = AreaProvider();
+            if (area < 1e-12)
+            {
+                SetClosedWall();
+                return;
+            }
+
+            // Retrieve volume state
+            double volP = VolumePort.Pressure;
+            double volRho = VolumePort.Density;
+            double volT = VolumePort.Temperature;
+            double volH = VolumePort.SpecificEnthalpy;
+
+            // Retrieve pipe interior state at the connected end
+            (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd
+                ? Pipe.GetInteriorStateLeft()
+                : Pipe.GetInteriorStateRight();
+
+            // Determine upstream/downstream: if volume pressure > pipe pressure, flow is out of volume (negative into volume).
+            bool flowOutOfVolume = volP > pipeP;
+            double pUp, rhoUp, TUp, pDown;
+            if (flowOutOfVolume)
+            {
+                pUp = volP; rhoUp = volRho; TUp = volT; pDown = pipeP;
+            }
+            else
+            {
+                // Pipe is upstream
+                pUp = pipeP; rhoUp = pipeRho; TUp = pipeP / (pipeRho * GasConstant); // temperature from pipe
+                pDown = volP;
+            }
+
+            // Compute isentropic nozzle flow
+            IsentropicOrifice.Compute(pUp, rhoUp, TUp, Gamma, GasConstant, pDown, area, DischargeCoefficient,
+                                      out double mdotUpstreamToDown, out double rhoFace, out double uFace, out double pFace);
+
+            // mdotUpstreamToDown is positive from upstream to downstream.
+            // Convert to mass flow into volume (positive mdot = into volume).
+            double mdotVolume;
+            if (flowOutOfVolume)
+                mdotVolume = -mdotUpstreamToDown;   // out of volume is negative
+            else
+                mdotVolume = mdotUpstreamToDown;    // into volume is positive
+
+            // Clamp mass flow to available mass in volume (if it is a Volume0D)
+            if (VolumePort.Owner is Volume0D vol)
+            {
+                double maxMdot = vol.Mass / dtSub;
+                if (mdotVolume > maxMdot) mdotVolume = maxMdot;
+                if (mdotVolume < -maxMdot) mdotVolume = -maxMdot;
+            }
+
+            // Apply ghost state to pipe
+            if (IsPipeLeftEnd)
+                Pipe.SetGhostLeft(rhoFace, uFace, pFace);
+            else
+                Pipe.SetGhostRight(rhoFace, uFace, pFace);
+
+            // Store results
+            LastMassFlowRate = mdotVolume;
+            LastFaceDensity = rhoFace;
+            LastFaceVelocity = uFace;
+            LastFacePressure = pFace;
+
+            // Set port flow rates for volume integration
+            VolumePort.MassFlowRate = mdotVolume;
+            if (mdotVolume >= 0)
+            {
+                // Inflow: enthalpy comes from upstream (pipe)
+                double pPipe = pipeP;
+                double rhoPipe = pipeRho;
+                VolumePort.SpecificEnthalpy = Gamma / (Gamma - 1.0) * pPipe / rhoPipe;
+            }
+            else
+            {
+                // Outflow: volume's own specific enthalpy
+                VolumePort.SpecificEnthalpy = volH;
+            }
+        }
+
+        private void SetClosedWall()
+        {
+            var (rInt, uInt, pInt) = IsPipeLeftEnd
+                ? Pipe.GetInteriorStateLeft()
+                : Pipe.GetInteriorStateRight();
+
+            if (IsPipeLeftEnd)
+                Pipe.SetGhostLeft(rInt, -uInt, pInt);
+            else
+                Pipe.SetGhostRight(rInt, -uInt, pInt);
+
+            LastMassFlowRate = 0.0;
+            LastFaceDensity = rInt;
+            LastFaceVelocity = 0.0;
+            LastFacePressure = pInt;
+            VolumePort.MassFlowRate = 0.0;
+            // Keep specific enthalpy as is (not used)
+        }
+    }
+}

Core/OutdoorExhaustReverb.cs

@@ -0,0 +1,187 @@
+using System;
+
+namespace FluidSim.Core
+{
+    public class OutdoorExhaustReverb
+    {
+        // ========== Early reflection delays (stereo: left/right) ==========
+        private readonly DelayLine groundL, groundR;
+        private readonly DelayLine wall1L, wall1R;
+        private readonly DelayLine wall2L, wall2R;
+
+        // ========== Diffuse tail FDNs (left/right each with 8 channels) ==========
+        private const int FDN_CHANNELS = 8;
+        private readonly DelayLine[] fdnL, fdnR;
+        private readonly float[] stateL, stateR;
+        private readonly OrthonormalMixer mixerL, mixerR;
+        private readonly LowPassFilter[] filterL, filterR;
+
+        public float DryMix      { get; set; } = 1.0f;
+        public float EarlyMix    { get; set; } = 0.5f;
+        public float TailMix     { get; set; } = 0.9f;
+        public float Feedback    { get; set; } = 0.75f;   // safe range 0.7‑0.9
+        public float DampingFreq { get; set; } = 6000f;    // Hz
+
+        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 (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/Solver.cs

@@ -1,170 +1,84 @@
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using FluidSim.Components;
 using FluidSim.Interfaces;
 
 namespace FluidSim.Core
 {
+    /// <summary>
+    /// Top‑level solver that owns all components and couplings,
+    /// orchestrates sub‑stepping, and exposes states for audio.
+    /// </summary>
     public class Solver
     {
-        private readonly List<Volume0D> _volumes = new();
+        private readonly List<IComponent> _components = new();
-        private readonly List<Pipe1D> _pipes = new();
+        private readonly List<OrificeLink> _orificeLinks = new();
-        private readonly List<Connection> _connections = new();
+        private readonly List<Junction> _junctions = new();
+        private readonly List<OpenEndLink> _openEndLinks = new();
 
         private double _dt;
 
-        public void AddVolume(Volume0D v) => _volumes.Add(v);
-        public void AddPipe(Pipe1D p) => _pipes.Add(p);
-        public void AddConnection(Connection c) => _connections.Add(c);
         public void SetTimeStep(double dt) => _dt = dt;
 
+        public void AddComponent(IComponent component) => _components.Add(component);
+        public void AddOrificeLink(OrificeLink link) => _orificeLinks.Add(link);
+        public void AddJunction(Junction junction) => _junctions.Add(junction);
+        public void AddOpenEndLink(OpenEndLink link) => _openEndLinks.Add(link);
+
+        // Convenience: first pipe’s port B mass flow (often the exhaust)
+        public double ExhaustMassFlow
+        {
+            get
+            {
+                var pipes = _components.OfType<Pipe1D>().ToList();
+                if (pipes.Count > 0)
+                    return Math.Abs(pipes[0].PortB.MassFlowRate);
+                return 0.0;
+            }
+        }
+
         /// <summary>
-        /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right).
+        /// Advance the whole system by one global time step.
         /// </summary>
-        public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0)
+        public void Step()
         {
-            if (isA)
+            var pipes = _components.OfType<Pipe1D>().ToList();
-            {
+            if (pipes.Count == 0) return;
-                pipe.SetABoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetAAmbientPressure(ambientPressure);
-            }
-            else
-            {
-                pipe.SetBBoundaryType(type);
-                if (type == BoundaryType.OpenEnd)
-                    pipe.SetBAmbientPressure(ambientPressure);
-            }
-        }
 
-        public float Step()
+            // 1. Determine sub‑step count (max CFL over all pipes)
-        {
-            // 1. Volumes publish state
-            foreach (var v in _volumes)
-                v.PushStateToPort();
-
-            // 2. Set volume BCs for volume‑coupled ends
-            foreach (var conn in _connections)
-            {
-                if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                {
-                    var pipe = GetPipe(conn.PortA);
-                    bool isA = pipe.PortA == conn.PortA;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortA, conn.PortB);
-                }
-                else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                {
-                    var pipe = GetPipe(conn.PortB);
-                    bool isA = pipe.PortB == conn.PortB;
-                    if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                        (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                        SetVolumeBC(conn.PortB, conn.PortA);
-                }
-            }
-
-            // 3. Sub‑steps
             int nSub = 1;
-            foreach (var p in _pipes)
+            foreach (var p in pipes)
                 nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
             double dtSub = _dt / nSub;
 
+            // 2. Sub‑step loop
             for (int sub = 0; sub < nSub; sub++)
             {
-                foreach (var p in _pipes)
+                // a) Resolve all orifice links (volume ↔ pipe)
+                foreach (var link in _orificeLinks)
+                    link.Resolve(dtSub);
+
+                // b) Resolve all open‑end links (pipe → atmosphere)
+                foreach (var link in _openEndLinks)
+                    link.Resolve(dtSub);
+
+                // c) Resolve all junctions (pipe ↔ pipe)
+                foreach (var junc in _junctions)
+                    junc.Resolve(dtSub);
+
+                // d) Advance all pipes
+                foreach (var p in pipes)
                     p.SimulateSingleStep(dtSub);
-
-                foreach (var conn in _connections)
-                {
-                    if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortA);
-                        bool isA = pipe.PortA == conn.PortA;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortA, conn.PortB, dtSub);
-                    }
-                    else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                    {
-                        var pipe = GetPipe(conn.PortB);
-                        bool isA = pipe.PortB == conn.PortB;
-                        if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                            (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                            TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
-                    }
             }
 
-                if (sub < nSub - 1)
+            // 3. Clear ghost flags
-                {
+            foreach (var p in pipes)
-                    foreach (var v in _volumes)
+                p.ClearGhostFlags();
-                        v.PushStateToPort();
 
-                    foreach (var conn in _connections)
+            // 4. Integrate non‑pipe components (volumes, atmosphere, etc.)
-                    {
+            foreach (var comp in _components)
-                        if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
+                comp.UpdateState(_dt);
-                        {
-                            var pipe = GetPipe(conn.PortA);
-                            bool isA = pipe.PortA == conn.PortA;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortA, conn.PortB);
-                        }
-                        else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
-                        {
-                            var pipe = GetPipe(conn.PortB);
-                            bool isA = pipe.PortB == conn.PortB;
-                            if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) ||
-                                (!isA && pipe.BBCType == BoundaryType.VolumeCoupling))
-                                SetVolumeBC(conn.PortB, conn.PortA);
-                        }
-                    }
-                }
-            }
-
-            // 5. Audio samples (none for now, but placeholder)
-            var audioSamples = new List<float>();
-            foreach (var conn in _connections)
-            {
-                if (conn is SoundConnection sc)
-                    audioSamples.Add(sc.GetAudioSample());
-            }
-
-            // 6. Clear BC flags
-            foreach (var p in _pipes)
-                p.ClearBC();
-
-            return SoundProcessor.MixAndClip(audioSamples.ToArray());
-        }
-
-        private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
-        private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
-        private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
-        private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p);
-
-        private void SetVolumeBC(Port pipePort, Port volPort)
-        {
-            var pipe = GetPipe(pipePort);
-            if (pipe == null) return;
-            bool isA = pipe.PortA == pipePort;
-            if (isA)
-                pipe.SetAVolumeState(volPort.Density, volPort.Pressure);
-            else
-                pipe.SetBVolumeState(volPort.Density, volPort.Pressure);
-        }
-
-        private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
-        {
-            double mdot = pipePort.MassFlowRate;
-            volPort.MassFlowRate = -mdot;
-
-            if (mdot < 0)  // pipe → volume
-            {
-                volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
-            }
-            // else volume’s own enthalpy (from PushStateToPort) is used
-
-            GetVolume(volPort)?.Integrate(dtSub);
         }
     }
 }

Core/SoundProcessor.cs

@@ -1,23 +1,48 @@
-namespace FluidSim.Core
+using System;
-{
+using FluidSim.Interfaces;
-    /// <summary>
-    /// Mixes multiple audio samples and applies a soft‑clipping tanh.
-    /// </summary>
-    public static class SoundProcessor
-    {
-        /// <summary>Overall gain applied after mixing (before tanh).</summary>
-        public static float MasterGain { get; set; } = 0.01f;
 
-        /// <summary>
+namespace FluidSim.Core
-        /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1].
-        /// </summary>
-        public static float MixAndClip(params float[] samples)
 {
-            float sum = 0f;
+    public class SoundProcessor
-            foreach (float s in samples)
+    {
-                sum += s;
+        private readonly double dt;
-            sum *= MasterGain;
+        private readonly double scaleFactor;  // 1 / (4π r) and a user gain
-            return sum;
+        private double prevMassFlowOut;
+
+        // Simple low‑pass for derivative smoothing (≈ 2‑3 ms)
+        private double smoothDMdt;
+        private readonly double alpha;
+
+        public float Gain { get; set; } = 1.0f;
+
+        public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0)
+        {
+            dt = 1.0 / sampleRate;
+            scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters);
+
+            // Smoothing time constant ~ 2 ms, blocks single‑sample spikes
+            double tau = 0.002;
+            alpha = Math.Exp(-dt / tau);
+        }
+
+        public float Process(Port port)
+        {
+            // Outflow mass flow (positive = leaving pipe)
+            double flowOut = -port.MassFlowRate;
+
+            // Derivative
+            double rawDerivative = (flowOut - prevMassFlowOut) / dt;
+            prevMassFlowOut = flowOut;
+
+            // Smooth the derivative to kill isolated spikes
+            smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative;
+
+            // Far‑field monopole pressure
+            double pressure = smoothDMdt * scaleFactor * Gain;
+
+            // Soft clip to ±1 for audio output (safe limit)
+            float sample = (float)Math.Tanh(pressure);
+            return sample;
         }
     }
 }

FluidSim.csproj

@@ -5,7 +5,7 @@
     <TargetFramework>net10.0</TargetFramework>
     <ImplicitUsings>enable</ImplicitUsings>
     <Nullable>enable</Nullable>
-    <PublishAot>true</PublishAot>
+    <PublishAot>false</PublishAot>
     <InvariantGlobalization>true</InvariantGlobalization>
   </PropertyGroup>
 

Interfaces/Connection.cs

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

Interfaces/IComponent.cs

@@ -0,0 +1,19 @@
+using System.Collections.Generic;
+
+namespace FluidSim.Interfaces
+{
+    /// <summary>
+    /// Minimal interface for all simulation components that have ports.
+    /// </summary>
+    public interface IComponent
+    {
+        /// <summary>All ports exposed by this component.</summary>
+        IReadOnlyList<Port> Ports { get; }
+
+        /// <summary>
+        /// Called once per global time step to update the component's internal state
+        /// using the port flow data accumulated during sub‑steps.
+        /// </summary>
+        void UpdateState(double dt);
+    }
+}

Interfaces/Port.cs

@@ -1,18 +1,28 @@
 namespace FluidSim.Interfaces
 {
+    /// <summary>
+    /// A fluid port that belongs to a component. The solver writes mass flow,
+    /// enthalpy, etc. here, and reads pressure, density, etc.
+    /// </summary>
     public class Port
     {
+        // Set by the solver during coupling resolution
+        public double MassFlowRate;         // kg/s, positive INTO the component that owns this port
+        public double SpecificEnthalpy;     // J/kg, enthalpy of the fluid crossing the port (inflow direction)
+
+        // Set by the owning component after integration to reflect its current state
         public double Pressure;             // Pa
-        public double MassFlowRate;         // kg/s, positive INTO the component
-        public double SpecificEnthalpy;     // J/kg, enthalpy of fluid entering this port
         public double Density;              // kg/m³
         public double Temperature;          // K
 
+        // Link back to owner (optional, but handy for debugging)
+        public object? Owner { get; set; }
+
         public Port()
         {
-            Pressure = 101325.0;
             MassFlowRate = 0.0;
             SpecificEnthalpy = 0.0;
+            Pressure = 101325.0;
             Density = 1.225;
             Temperature = 300.0;
         }

Interfaces/SoundConnection.cs

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

Program.cs

@@ -3,6 +3,7 @@ using SFML.Window;
 using SFML.System;
 using System.Diagnostics;
 using FluidSim.Core;
+using FluidSim.Tests;
 
 namespace FluidSim;
 
@@ -12,37 +13,37 @@ public class Program
     private const double DrawFrequency = 60.0;
     private static Scenario scenario;
 
-    // Speed control
+    // Speed control (existing + new throttle)
-    //private static double desiredSpeed = 1.0;
+    private static double desiredSpeed = 0.01;
-    private static double desiredSpeed = 0.0001;
     private static double currentSpeed = desiredSpeed;
     private const double MinSpeed = 0.0001;
     private const double MaxSpeed = 1.0;
     private const double ScrollFactor = 1.1;
 
-    // Space‑toggle state
+    private static double lastDesiredSpeed = 0.1;
-    private static double lastDesiredSpeed = 0.1;   // remembers the last non‑1.0 scroll speed
+    private static bool isRealTime = false;
-    private static bool isRealTime = true;          // true when desiredSpeed == 1.0
+
+    // Throttle smoothing
+    private static double targetThrottle = 0.0;   // 1.0 when W is pressed, 0.0 otherwise
+    private static double currentThrottle = 0.0;
+    private const double ThrottleSmoothing = 20.0; // rate of change
 
     private static volatile bool running = true;
 
     public static void Main()
     {
         var mode = new VideoMode(new Vector2u(1280, 720));
-        var window = new RenderWindow(mode, "Pipe Resonator");
+        var window = new RenderWindow(mode, "FluidSim - Engine (W = throttle)");
         window.SetVerticalSyncEnabled(true);
         window.Closed += (_, _) => { running = false; window.Close(); };
         window.MouseWheelScrolled += OnMouseWheel;
         window.KeyPressed += OnKeyPressed;
 
         var soundEngine = new SoundEngine(bufferCapacity: 16384);
-        soundEngine.Volume = 70;
+        soundEngine.Volume = 100;
         soundEngine.Start();
 
-        //scenario = new PipeResonatorScenario();
+        scenario = new TestScenario();
-        //scenario = new HelmholtzResonatorScenario();
-        scenario = new SodShockTubeScenario();
-
         scenario.Initialize(SampleRate);
 
         var stopwatch = Stopwatch.StartNew();
@@ -50,17 +51,13 @@ public class Program
         double drawInterval = 1.0 / DrawFrequency;
         double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
 
-        // Resampling buffer
+        var simBuffer = new List<float>(4096);
-        List<float> simBuffer = new List<float>(4096);
         double readIndex = 0.0;
-
         for (int i = 0; i < 4; i++)
             simBuffer.Add(scenario.Process());
 
         long totalSimSteps = simBuffer.Count;
         long totalOutputSamples = 0;
-
-        double lastRealTime = stopwatch.Elapsed.TotalSeconds;
         const int outputChunk = 256;
         float[] outputBuf = new float[outputChunk];
 
@@ -69,17 +66,16 @@ public class Program
             window.DispatchEvents();
 
             double currentRealTime = stopwatch.Elapsed.TotalSeconds;
-            double dtSpeed = currentRealTime - lastSpeedUpdateTime;
+            double dtClock = currentRealTime - lastSpeedUpdateTime;
             lastSpeedUpdateTime = currentRealTime;
 
-            // Smoothly transition currentSpeed → desiredSpeed
+            // Smooth simulation speed
-            // When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp
+            double speedSmoothing = 8.0;
-            double smoothingRate = 8.0;   // higher = faster catch‑up
+            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-speedSmoothing * dtClock));
-            currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
 
-            // ---------- Generate audio ----------
+
+            // Generate audio
             double targetAudioClock = currentRealTime + 0.05;
-
             while (totalOutputSamples < targetAudioClock * SampleRate && running)
             {
                 int toGenerate = (int)Math.Min(
@@ -101,11 +97,9 @@ public class Program
                     int i0 = (int)readIndex;
                     int i1 = i0 + 1;
                     double frac = readIndex - i0;
-
                     float y0 = simBuffer[Math.Clamp(i0, 0, simBuffer.Count - 1)];
                     float y1 = simBuffer[Math.Clamp(i1, 0, simBuffer.Count - 1)];
                     outputBuf[i] = (float)(y0 + (y1 - y0) * frac);
-
                     readIndex += currentSpeed;
 
                     while (readIndex >= 1.0 && simBuffer.Count > 2)
@@ -117,23 +111,23 @@ public class Program
 
                 int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
                 totalOutputSamples += accepted;
-
                 if (accepted < toGenerate)
                     break;
             }
 
-            // ---------- Drawing & title ----------
+            // Drawing & title
             if (currentRealTime - lastDrawTime >= drawInterval)
             {
                 double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate);
                 double simTime = totalSimSteps / (double)SampleRate;
                 string toggleHint = isRealTime ? "[Space] slow mo" : "[Space] real time";
+                string throttleHint = Keyboard.IsKeyPressed(Keyboard.Key.W) ? "W held" : "W released";
                 window.SetTitle(
-                    $"{toggleHint}  Sim: {simTime:F2}s | " +
+                    $"{toggleHint}  {throttleHint}  " +
-                    $"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
+                    $"Thr: {currentThrottle:F2}  " +
-                    $"Actual: {actualSpeed:F2}x"
+                    $"Speed: {currentSpeed:F3}x → {desiredSpeed:F3}x  " +
+                    $"Act: {actualSpeed:F2}x"
                 );
-
                 window.Clear(Color.Black);
                 scenario.Draw(window);
                 window.Display();
@@ -145,22 +139,18 @@ public class Program
         window.Dispose();
     }
 
+    // (Mouse wheel, space toggle unchanged)
     private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
     {
         bool wasRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
-
         if (e.Delta > 0)
             desiredSpeed *= ScrollFactor;
         else if (e.Delta < 0)
             desiredSpeed /= ScrollFactor;
 
         desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed);
-
-        // Update the remembered slow-mo speed (unless we are exactly at 1.0)
         if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6)
             lastDesiredSpeed = desiredSpeed;
-
-        // Update isRealTime flag
         isRealTime = Math.Abs(desiredSpeed - 1.0) < 1e-6;
     }
 
@@ -169,15 +159,9 @@ public class Program
         if (e.Code == Keyboard.Key.Space)
         {
             if (isRealTime)
-            {
-                // Switch to the remembered slow speed
                 desiredSpeed = lastDesiredSpeed;
-            }
             else
-            {
-                // Switch back to real time
                 desiredSpeed = 1.0;
-            }
             isRealTime = !isRealTime;
         }
     }

Scenarios/HelmholtzResonatorScenario.cs

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

Scenarios/PipeResonatorScenario.cs

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

Scenarios/Scenario.cs

@@ -1,23 +1,121 @@
-using SFML.Graphics;
+using System;
+using SFML.Graphics;
+using SFML.System;
+using FluidSim.Components;
 
-namespace FluidSim.Core
+namespace FluidSim.Tests
 {
     public abstract class Scenario
     {
-        /// <summary>
+        /// <summary>Initialize the scenario with a given audio sample rate.</summary>
-        /// Initialize the scenario with a given audio sample rate.
-        /// </summary>
         public abstract void Initialize(int sampleRate);
 
-        /// <summary>
+        /// <summary>Advance one simulation step and return an audio sample.</summary>
-        /// Advance one simulation step and return an audio sample.
-        /// The step size is 1 / sampleRate seconds.
-        /// </summary>
         public abstract float Process();
 
-        /// <summary>
+        /// <summary>Draw the current simulation state onto the given SFML render target.</summary>
-        /// Draw the current simulation state onto the given SFML render target.
-        /// </summary>
         public abstract void Draw(RenderWindow target);
+
+        // ---------- Shared drawing helpers ----------
+
+        protected const double AmbientPressure = 101325.0;
+
+        /// <summary>Blue (low) → Green (ambient) → Red (high).</summary>
+        protected Color PressureColor(double pressure)
+        {
+            double range = AmbientPressure * 0.05;   // ±5% gives full colour swing
+            double t = (pressure - AmbientPressure) / range;
+            t = Math.Clamp(t, -1.0, 1.0);
+
+            byte r, g, b;
+            if (t < 0)
+            {
+                double factor = -t;
+                r = 0;
+                g = (byte)(255 * (1 - factor));
+                b = (byte)(255 * factor);
+            }
+            else
+            {
+                double factor = t;
+                r = (byte)(255 * factor);
+                g = (byte)(255 * (1 - factor));
+                b = 0;
+            }
+            return new Color(r, g, b);
+        }
+
+        /// <summary>
+        /// Draws the pipe as a smooth triangle‑strip whose radius varies with cell pressure.
+        /// </summary>
+        protected void DrawPipe(RenderWindow target, Pipe1D pipe, float pipeCenterY, float pipeStartX, float pipeEndX)
+        {
+            int n = pipe.CellCount;
+            if (n < 2) return;
+
+            float pipeLengthPx = pipeEndX - pipeStartX;
+            float dx = pipeLengthPx / (n - 1);   // spacing between cell centres
+
+            float baseRadius = 25f;
+            float rangeFactor = 1f;
+            float scaleFactor = 5f;
+
+            // ----- smoothstep helper -----
+            static float SmoothStep(float edge0, float edge1, float x)
+            {
+                float t = Math.Clamp((x - edge0) / (edge1 - edge0), 0f, 1f);
+                return t * t * (3f - 2f * t);
+            }
+
+            // ----- Pre‑compute cell positions and radii -----
+            var centers = new float[n];
+            var radii = new float[n];
+            for (int i = 0; i < n; i++)
+            {
+                double p = pipe.GetCellPressure(i);
+                float deviation = (float)Math.Tanh((p - AmbientPressure) / AmbientPressure / rangeFactor);
+                radii[i] = baseRadius * (1f + deviation * scaleFactor);
+                if (radii[i] < 2f) radii[i] = 2f;
+                centers[i] = pipeStartX + i * dx;
+            }
+
+            // ----- Build triangle‑strip vertices -----
+            int segmentsPerCell = 8;
+            int totalPoints = n + (n - 1) * segmentsPerCell;
+            Vertex[] stripVertices = new Vertex[totalPoints * 2];
+            int idx = 0;
+
+            for (int i = 0; i < n; i++)
+            {
+                float x = centers[i];
+                float r = radii[i];
+                double p = pipe.GetCellPressure(i);
+                Color col = PressureColor(p);
+
+                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col);
+                stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col);
+
+                if (i < n - 1)
+                {
+                    for (int s = 1; s <= segmentsPerCell; s++)
+                    {
+                        float t = s / (float)segmentsPerCell;
+                        float st = SmoothStep(0f, 1f, t);
+                        float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
+                        float ri = radii[i] + (radii[i + 1] - radii[i]) * st;
+                        double pi = pipe.GetCellPressure(i) * (1 - t) + pipe.GetCellPressure(i + 1) * t;
+                        Color coli = PressureColor(pi);
+
+                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli);
+                        stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), coli);
+                    }
+                }
+            }
+
+            var pipeMesh = new VertexArray(PrimitiveType.TriangleStrip, (uint)stripVertices.Length);
+            for (int i = 0; i < stripVertices.Length; i++)
+                pipeMesh[(uint)i] = stripVertices[i];
+            target.Draw(pipeMesh);
+        }
     }
 }

Scenarios/SodShockTubeScenario.cs

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

Scenarios/TestScenario.cs

@@ -0,0 +1,96 @@
+using System;
+using SFML.Graphics;
+using SFML.System;
+using FluidSim.Components;
+using FluidSim.Core;
+
+namespace FluidSim.Tests
+{
+    public class TestScenario : Scenario
+    {
+        private Solver solver;
+        private Volume0D volume;
+        private Pipe1D pipe;
+        private Atmosphere atmosphere;
+        private OrificeLink orificeLink;
+        private OpenEndLink openEndLink;
+        private int stepCount;
+
+        public override void Initialize(int sampleRate)
+        {
+            double dt = 1.0 / sampleRate;
+
+            solver = new Solver();
+            solver.SetTimeStep(dt);
+
+            volume = new Volume0D(1e-3, 150000.0, 300.0);
+            solver.AddComponent(volume);
+
+            pipe = new Pipe1D(2.0, 1e-4, 20);
+            solver.AddComponent(pipe);
+
+            atmosphere = new Atmosphere();
+            solver.AddComponent(atmosphere);
+
+            // Volume → left pipe end (orifice)
+            var volPort = volume.CreatePort();
+            orificeLink = new OrificeLink(volPort, pipe, isPipeLeftEnd: true, areaProvider: () => 1e-5)
+            {
+                DischargeCoefficient = 0.62,
+                Gamma = volume.Gamma,
+                GasConstant = volume.GasConstant
+            };
+            solver.AddOrificeLink(orificeLink);
+
+            // Right pipe end → atmosphere (characteristic open‑end)
+            openEndLink = new OpenEndLink(pipe, isLeftEnd: false)
+            {
+                AmbientPressure = 101325.0,
+                Gamma = 1.4
+            };
+            solver.AddOpenEndLink(openEndLink);
+
+            stepCount = 0;
+            Console.WriteLine("TestScenario initialized with sampleRate = " + sampleRate);
+        }
+
+        public override float Process()
+        {
+            solver.Step();
+            stepCount++;
+
+            if (stepCount % 100 == 0)
+            {
+                double volPressure = volume.Pressure;
+                double volMass = volume.Mass;
+                double pipeLeftPressure = pipe.GetCellPressure(0);
+                double pipeRightPressure = pipe.GetCellPressure(pipe.CellCount - 1);
+                double mdotOrifice = orificeLink.LastMassFlowRate;
+                double mdotOpen = openEndLink.LastMassFlowRate;
+
+                Console.WriteLine($"Step {stepCount}:");
+                Console.WriteLine($"  Vol Pressure = {volPressure:F1} Pa, Mass = {volMass:E4} kg");
+                Console.WriteLine($"  Pipe left P = {pipeLeftPressure:F1} Pa, right P = {pipeRightPressure:F1} Pa");
+                Console.WriteLine($"  Orifice mdot = {mdotOrifice:E4} kg/s, Open‑end mdot = {mdotOpen:E4} kg/s");
+                Console.WriteLine();
+            }
+
+            // Audio sample from the open‑end mass flow
+            return (float)openEndLink.LastMassFlowRate;
+        }
+
+        public override void Draw(RenderWindow target)
+        {
+            float winWidth = target.GetView().Size.X;
+            float winHeight = target.GetView().Size.Y;
+
+            float pipeCenterY = winHeight / 2f;
+            float margin = 60f;
+            float pipeStartX = margin;
+            float pipeEndX = winWidth - margin;
+
+            // Use the shared pipe drawing from the base class
+            DrawPipe(target, pipe, pipeCenterY, pipeStartX, pipeEndX);
+        }
+    }
+}

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
-    {
-    }
-}