orifice confirmed working

Components/Volume0D.cs

@@ -4,23 +4,17 @@ using FluidSim.Interfaces;
 
 namespace FluidSim.Components
 {
-    /// <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 List<Port> Ports { get; } = new List<Port>();
 
-        public double Mass { get; private set; }
+        public double Mass { get; set; }           // made public setter
-        public double InternalEnergy { get; private set; }
+        public double InternalEnergy { get; set; }  // made public setter
         public double Volume { get; set; }
         public double Dvdt { get; set; }
         public double Gamma { get; set; } = 1.4;
         public double GasConstant { get; set; } = 287.0;
 
-        // Ambient pressure used for emergency refill – default 101325 Pa
         public double AmbientPressure { get; set; } = 101325.0;
 
         // Derived quantities
@@ -42,11 +36,9 @@ namespace FluidSim.Components
             InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
         }
 
-        /// <summary>Add a new port and initialise it to the volume's current state.</summary>
         public Port CreatePort()
         {
             var port = new Port { Owner = this };
-            // Set the port state immediately to avoid a mismatch before the first integration
             port.Pressure = Pressure;
             port.Density = Density;
             port.Temperature = Temperature;
@@ -56,9 +48,18 @@ namespace FluidSim.Components
         }
 
         /// <summary>
-        /// Integrate over dt using the MassFlowRate and SpecificEnthalpy
+        /// Set the pressure to a specific value while keeping the current temperature constant.
-        /// that have been set on each port during the coupling resolution phase.
+        /// Updates Mass and InternalEnergy accordingly.
         /// </summary>
+        public void SetPressure(double pressure)
+        {
+            double V = Math.Max(Volume, 1e-12);
+            double currentT = Temperature;   // current temperature before changes
+            double rho = pressure / (GasConstant * currentT);
+            Mass = rho * V;
+            InternalEnergy = pressure * V / (Gamma - 1.0);
+        }
+
         public void UpdateState(double dt)
         {
             double totalMdot = 0.0;
@@ -67,17 +68,15 @@ namespace FluidSim.Components
             foreach (var port in Ports)
             {
                 totalMdot += port.MassFlowRate;
-                // mdot * h gives energy flow: positive mdot = inflow, negative = outflow
                 totalEdot += port.MassFlowRate * port.SpecificEnthalpy;
             }
 
             double dm = totalMdot * dt;
-            double dE = totalEdot * dt - Pressure * Dvdt * dt;   // piston work
+            double dE = totalEdot * dt - Pressure * Dvdt * dt;
 
             Mass += dm;
             InternalEnergy += dE;
 
-            // Safeguards: keep at least 1 µg of gas at a small pressure
             double V = Math.Max(Volume, 1e-12);
             if (Mass < 1e-9)
             {
@@ -89,18 +88,16 @@ namespace FluidSim.Components
                 InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
             }
 
-            // 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
+                port.SpecificEnthalpy = h;
             }
         }
 

Core/OpenEndLink.cs

@@ -3,14 +3,19 @@ using FluidSim.Components;
 
 namespace FluidSim.Core
 {
+    /// <summary>
+    /// Characteristic open‑end boundary condition after Jones (1978).
+    /// For all subsonic flow (outflow and inflow), the ghost state is derived
+    /// from the isentropic expansion to ambient pressure, using the pipe's entropy,
+    /// and the outgoing Riemann invariant. This avoids a density jump at flow reversal.
+    /// Supersonic outflow extrapolates the interior state.
+    /// </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;
-        public double GasConstant { get; set; } = 287.0;
-        public double AmbientTemperature { get; set; } = 300.0;
 
         public double LastMassFlowRate { get; private set; }
         public double LastFaceDensity { get; private set; }
@@ -33,61 +38,63 @@ namespace FluidSim.Core
             double gm1 = gamma - 1.0;
             double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12));
             double pAmb = AmbientPressure;
-            double rhoAmb = pAmb / (GasConstant * AmbientTemperature);
+
-            double aAmb = Math.Sqrt(gamma * pAmb / rhoAmb);
+            // Riemann invariants
+            double J_plus  = uInt + 2.0 * cInt / gm1;
+            double J_minus = uInt - 2.0 * cInt / gm1;
 
             double rhoGhost, uGhost, pGhost;
 
-            // ----- Supersonic outflow: extrapolate interior -----
+            // ---- Subsonic branch (used for both outflow and inflow) ----
-            bool supersonicOut = IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
+            // Isentropic expansion to ambient pressure using pipe's entropy
-            if (supersonicOut)
+            double s = pInt / Math.Pow(rhoInt, gamma);        // entropy constant
+            double rhoIso = Math.Pow(pAmb / s, 1.0 / gamma);
+            double cIso = Math.Sqrt(gamma * pAmb / Math.Max(rhoIso, 1e-12));
+            double uIso = IsLeftEnd
+                ? (J_minus + 2.0 * cIso / gm1)
+                : (J_plus  - 2.0 * cIso / gm1);
+
+            // Check for supersonic outflow: if the isentropic velocity exceeds the speed of sound,
+            // the flow is supersonic and we extrapolate the interior state.
+            bool supersonic = IsLeftEnd
+                ? (uInt <= -cInt)   // left end: outflow is when u < -c
+                : (uInt >= cInt);   // right end: outflow is when u > c
+
+            // Extra check: if the isentropic velocity is supersonic in the outflow direction,
+            // also treat as supersonic (this can happen when the interior pressure is very high).
+            if (!supersonic)
             {
+                if (IsLeftEnd)
+                    supersonic = uIso <= -cIso;
+                else
+                    supersonic = uIso >= cIso;
+            }
+
+            if (supersonic)
+            {
+                // Supersonic outflow – extrapolate interior
                 rhoGhost = rhoInt;
                 uGhost   = uInt;
                 pGhost   = pInt;
             }
             else
             {
-                // Riemann invariants
+                // Subsonic flow – use the isentropic state
-                double J_plus  = uInt + 2.0 * cInt / gm1;
+                rhoGhost = rhoIso;
-                double J_minus = uInt - 2.0 * cInt / gm1;
+                uGhost   = uIso;
-
+                pGhost   = pAmb;
-                // Trial subsonic outflow ghost state
-                double s = pInt / Math.Pow(rhoInt, gamma);
-                double rhoOut = Math.Pow(pAmb / s, 1.0 / gamma);
-                double cOut = Math.Sqrt(gamma * pAmb / rhoOut);
-                double uOut = IsLeftEnd
-                    ? (J_minus + 2.0 * cOut / gm1)
-                    : (J_plus  - 2.0 * cOut / gm1);
-
-                bool outflowPossible = IsLeftEnd ? (uOut <= 0) : (uOut >= 0);
-
-                if (outflowPossible)
-                {
-                    // Subsonic outflow
-                    pGhost = pAmb;
-                    rhoGhost = rhoOut;
-                    uGhost = uOut;
-                }
-                else
-                {
-                    // Subsonic inflow (ambient reservoir model)
-                    pGhost = pAmb;
-                    rhoGhost = rhoAmb;
-                    uGhost = IsLeftEnd
-                        ? (J_minus + 2.0 * aAmb / gm1)
-                        : (J_plus  - 2.0 * aAmb / gm1);
-                }
             }
 
+            // Apply ghost to pipe
             if (IsLeftEnd)
                 Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost);
             else
                 Pipe.SetGhostRight(rhoGhost, uGhost, pGhost);
 
+            // Mass flow out of the pipe (positive = leaving)
             double area = Pipe.Area;
             double mdot = rhoGhost * uGhost * area;
-            if (IsLeftEnd) mdot = -mdot;
+            if (IsLeftEnd) mdot = -mdot;   // left end: positive u is into pipe, so out is -u
             LastMassFlowRate = mdot;
             LastFaceDensity = rhoGhost;
             LastFaceVelocity = uGhost;

Core/Solver.cs

@@ -6,10 +6,6 @@ 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<IComponent> _components = new();
@@ -19,6 +15,9 @@ namespace FluidSim.Core
 
         private double _dt;
 
+        /// <summary>CFL target for sub‑stepping (0.3‑0.8). Lower values are safer for shocks.</summary>
+        public double CflTarget { get; set; } = 0.8;
+
         public void SetTimeStep(double dt) => _dt = dt;
 
         public void AddComponent(IComponent component) => _components.Add(component);
@@ -26,57 +25,38 @@ namespace FluidSim.Core
         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>
-        /// Advance the whole system by one global time step.
-        /// </summary>
         public void Step()
         {
             var pipes = _components.OfType<Pipe1D>().ToList();
             if (pipes.Count == 0) return;
 
-            // 1. Determine sub‑step count (max CFL over all pipes)
             int nSub = 1;
             foreach (var p in pipes)
-                nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
+                nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt, CflTarget));
             double dtSub = _dt / nSub;
 
-            // 2. Sub‑step loop
+            const int maxSubSteps = 10000;
+            if (nSub > maxSubSteps)
+            {
+                Console.WriteLine($"Warning: required sub‑steps {nSub} exceeds limit. Simulation stopped.");
+                return;
+            }
+
             for (int sub = 0; sub < nSub; sub++)
             {
-                // 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);
             }
 
-            // 3. Clear ghost flags
             foreach (var p in pipes)
                 p.ClearGhostFlags();
 
-            // 4. Integrate non‑pipe components (volumes, atmosphere, etc.)
             foreach (var comp in _components)
                 comp.UpdateState(_dt);
         }

Program.cs

@@ -14,7 +14,7 @@ public class Program
     private static Scenario scenario;
 
     // Speed control (existing + new throttle)
-    private static double desiredSpeed = 0.001;
+    private static double desiredSpeed = 0.01;
     //private static double desiredSpeed = 1.0;
     private static double currentSpeed = desiredSpeed;
     private const double MinSpeed = 0.0001;

Scenarios/TestScenario.cs

@@ -3,60 +3,74 @@ using SFML.Graphics;
 using SFML.System;
 using FluidSim.Components;
 using FluidSim.Core;
+using FluidSim.Utils;
 
 namespace FluidSim.Tests
 {
     public class TestScenario : Scenario
     {
         private Solver solver;
+        private Volume0D volume;
         private Pipe1D pipe;
-        private OrificeLink closedEnd;    // left end – closed wall
+        private OrificeLink orifice;       // volume → pipe left
-        private OpenEndLink openEndLink;  // right end – atmosphere
+        private OpenEndLink openEnd;       // pipe right → atmosphere
         private SoundProcessor soundProcessor;
-        private OutdoorExhaustReverb reverb;
         private int stepCount;
-        private double simTime;            // elapsed simulation time (seconds)
+
-        private double pulseInterval = 0.1; // seconds between pulses
+        // Pressure reset scheduling
-        private double nextPulseTime;
+        private double simTime;
         private double dt;
+        private double resetInterval = 0.2;          // seconds between resets
+        private double nextResetTime;
+        private double targetPressure = 10 * Units.atm;
+        private double rampDuration = 0.002;          // 2 ms ramp
+        private double rampStartTime;
+        private double rampStartPressure;             // pressure when ramp started
+        private bool ramping;
 
         public override void Initialize(int sampleRate)
         {
             dt = 1.0 / sampleRate;
 
             soundProcessor = new SoundProcessor(sampleRate, 1);
-            soundProcessor.Gain = 10;
+            soundProcessor.Gain = 2.0f;   // lower gain to avoid clipping
-
-            reverb = new OutdoorExhaustReverb(sampleRate);
 
             solver = new Solver();
             solver.SetTimeStep(dt);
+            solver.CflTarget = 0.4;
 
-            // Pipe: 2 m long, 1 cm² area, 200 cells
+            volume = new Volume0D(1e-3, targetPressure, 300.0);
-            pipe = new Pipe1D(length: 2, area: 1e-4, cellCount: 100);
+            solver.AddComponent(volume);
+
+            pipe = new Pipe1D(1.0, 1e-4, 200);
+            pipe.EnergyRelaxationRate = 10;
             solver.AddComponent(pipe);
 
-            // Initially pipe at ambient conditions
+            var volPort = volume.CreatePort();
-            pipe.SetUniformState(1.225, 0.0, 101325.0);
+            double orificeArea = 1e-5;
+            orifice = new OrificeLink(volPort, pipe, isPipeLeftEnd: true,
+                                      areaProvider: () => orificeArea)
+            {
+                DischargeCoefficient = 0.62,
+                UseInertance = true,
+                EffectiveLength = 0.001
+            };
+            solver.AddOrificeLink(orifice);
 
-            // Left end: closed wall (area = 0 → reflective)
+            openEnd = new OpenEndLink(pipe, isLeftEnd: false)
-            closedEnd = new OrificeLink(null, pipe, isPipeLeftEnd: true, areaProvider: () => 0.0);
-            solver.AddOrificeLink(closedEnd);
-
-            // Right end: open to atmosphere
-            openEndLink = new OpenEndLink(pipe, isLeftEnd: false)
             {
                 AmbientPressure = 101325.0,
                 Gamma = 1.4
             };
-            solver.AddOpenEndLink(openEndLink);
+            solver.AddOpenEndLink(openEnd);
 
             stepCount = 0;
             simTime = 0.0;
-            nextPulseTime = pulseInterval;   // first pulse at 100 ms
+            nextResetTime = resetInterval;
+            ramping = false;
 
-            Console.WriteLine("Pulse reflection test – closed left, open right");
+            Console.WriteLine("Pressure reset test with smooth ramp");
-            Console.WriteLine("Pulse injected every 100 ms at left end (cell 0)");
+            Console.WriteLine($"Volume 1L, reset to {targetPressure} Pa every {resetInterval*1000} ms, ramp {rampDuration*1000} ms");
         }
 
         public override float Process()
@@ -65,59 +79,54 @@ namespace FluidSim.Tests
             stepCount++;
             simTime += dt;
 
-            // ---- Inject a pressure pulse at the closed end every 100 ms ----
+            // ---- Smooth pressure ramp ----
-            if (simTime >= nextPulseTime)
+            if (ramping)
             {
-                // Apply a Gaussian pulse to the first few cells
+                double elapsed = simTime - rampStartTime;
-                double ambientPressure = 101325.0;
+                if (elapsed >= rampDuration)
-                double pulseAmplitude = 20 * ambientPressure;   // 0.5 atm overpressure
-                double pulseWidth = 0.05;   // m (spread over a few cells)
-                int n = pipe.CellCount;
-                double dx = 2.0 / n;
-
-                // Only modify cells within 2*pulseWidth from the left end
-                int maxCell = Math.Min(5, n - 1);   // at most the first 5 cells
-                for (int i = 0; i <= maxCell; i++)
                 {
-                    double x = (i + 0.5) * dx;
+                    // Ramp finished, set exactly to target
-                    double P = pulseAmplitude * Math.Exp(-x * x / (pulseWidth * pulseWidth));
+                    volume.SetPressure(targetPressure);
-                    double currentP = pipe.GetCellPressure(i);
+                    ramping = false;
-                    double newP = P;
+                }
-                    // Update pressure, keeping density and velocity unchanged
+                else
-                    // We recompute total energy accordingly
+                {
-                    double rho = pipe.GetCellDensity(i);
+                    double frac = elapsed / rampDuration;
-                    double u = pipe.GetCellVelocity(i);
+                    double currentTarget = rampStartPressure + (targetPressure - rampStartPressure) * frac;
-                    double e = newP / ((1.4 - 1.0) * rho);
+                    volume.SetPressure(currentTarget);
-                    double E = rho * e + 0.5 * rho * u * u;
-                    pipe.SetCellState(i, rho, u, newP);
                 }
-                Console.WriteLine($"Pulse injected at t = {simTime:F3} s");
-                nextPulseTime += pulseInterval;
             }
 
-            // Audio from open‑end mass flow
+            // ---- Trigger a new reset ----
-            float sample = soundProcessor.Process(openEndLink);
+            if (!ramping && simTime >= nextResetTime)
-
-            // Log every 200 steps
-            if (stepCount % 1000 == 0)
             {
-                int leftIdx = 0;
+                rampStartPressure = volume.Pressure;   // current pressure before reset
-                int midIdx = pipe.CellCount / 2;
+                rampStartTime = simTime;
-                int rightIdx = pipe.CellCount - 1;
+                ramping = true;
-                double pL = pipe.GetCellPressure(leftIdx);
+                nextResetTime += resetInterval;
-                double pM = pipe.GetCellPressure(midIdx);
+            }
-                double pR = pipe.GetCellPressure(rightIdx);
 
-                Console.WriteLine($"Step {stepCount}: P_left={pL:F1} Pa, P_mid={pM:F1} Pa, P_right={pR:F1} Pa");
+            // Log every 500 steps
+            if (stepCount % 500 == 0)
+            {
+                double volP = volume.Pressure;
+                double pipeL = pipe.GetCellPressure(0);
+                double pipeR = pipe.GetCellPressure(pipe.CellCount - 1);
+                double mdotOrif = orifice.LastMassFlowRate;
+                double mdotOpen = openEnd.LastMassFlowRate;
+
+                Console.WriteLine($"Step {stepCount}: " +
+                    $"VolP={volP:F1} Pa, PipeL={pipeL:F1}, PipeR={pipeR:F1}, " +
+                    $"mdot_orif={mdotOrif:E4} kg/s, mdot_open={mdotOpen:E4} kg/s");
             }
 
             if (double.IsNaN(pipe.GetCellPressure(0)))
             {
-                Console.WriteLine("NaN detected – stopping simulation.");
+                Console.WriteLine("NaN detected – stopping.");
                 return 0f;
             }
 
-            return reverb.Process(sample);
+            return soundProcessor.Process(openEnd);
         }
 
         public override void Draw(RenderWindow target)