From 1489f278dc31388b73c101c20f73e09e5c3df90a Mon Sep 17 00:00:00 2001
From: max <maxwesterlund2@gmail.com>
Date: Sat, 9 May 2026 02:29:09 +0200
Subject: [PATCH] fixed under damping inertance

---
 Core/BoundarySystem.cs         | 32 +++++++++++++++++++-------------
 Scenarios/HelmholtzScenario.cs |  3 ++-
 2 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/Core/BoundarySystem.cs b/Core/BoundarySystem.cs
index c66876c..5264250 100644
--- a/Core/BoundarySystem.cs
+++ b/Core/BoundarySystem.cs
@@ -184,45 +184,51 @@ namespace FluidSim.Core
 
                 if (d.UseInertance)
                 {
-                    // ---- Inertance ODE with loss ----
                     float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho;
                     float inertance = rhoUp * d.EffectiveLength / MathF.Max(area, 1e-12f);
                     float dp = volP - pipeP;
-                    float Rlin = d.LossCoefficient;   // linear resistance
+                    float Rlin = d.LossCoefficient;
 
                     // Forward Euler
                     float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
                     float mdotNew = d.CurrentMdot + dmdot_dt * dt;
 
-                    // ---------- Symmetric flow limiter ----------
-                    // 1) limit by cavity mass
+                    // ---------- Symmetric flow limiters (existing) ----------
                     if (d.VolumePort.Owner is Volume0D vol0)
                     {
                         float maxOut = vol0.Mass / dt;
-                        if (mdotNew > maxOut) mdotNew = maxOut;   // cavity → pipe
-                        if (mdotNew < -maxOut) mdotNew = -maxOut;   // pipe → cavity
+                        if (mdotNew > maxOut) mdotNew = maxOut;
+                        if (mdotNew < -maxOut) mdotNew = -maxOut;
                     }
 
-                    // 2) limit by mass in the adjacent pipe cell
                     int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
                                               : _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
                     float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
-                    float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);   // uses the new public method
+                    float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
                     float pipeCellMass = pipeRhoAdj * area * pipeDxAdj;
                     float maxFromPipe = pipeCellMass / dt;
-                    if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;   // prevent emptying the pipe cell
+                    if (mdotNew < -maxFromPipe) mdotNew = -maxFromPipe;
+
+                    // ---------- Velocity clamp to Mach 0.9 ----------
+                    float rhoFacePrelim = mdotNew >= 0 ? volRho : pipeRho;
+                    float uFacePrelim = MathF.Abs(mdotNew) / MathF.Max(rhoFacePrelim * area, 1e-12f);
+                    float cUp = mdotNew >= 0 ? MathF.Sqrt(Gamma * Rgas * volT) : MathF.Sqrt(Gamma * Rgas * pipeT);
+                    float maxU = 0.9f * cUp;
+                    if (uFacePrelim > maxU)
+                    {
+                        uFacePrelim = maxU;
+                        mdotNew = rhoFacePrelim * uFacePrelim * area * (mdotNew >= 0 ? 1f : -1f);
+                    }
 
                     // NaN safety
                     if (float.IsNaN(mdotNew)) mdotNew = 0f;
 
-                    // Store for next step
                     d.CurrentMdot = mdotNew;
                     mdotFinal = mdotNew;
 
-                    // Ghost state: use the isentropic face pressure (pFaceEst) as reference,
-                    // but compute velocity from mdotFinal.
+                    // Ghost state
                     rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
-                    pFace = pFaceEst;   // approximate face pressure
+                    pFace = pFaceEst;
                     uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
                 }
                 else
diff --git a/Scenarios/HelmholtzScenario.cs b/Scenarios/HelmholtzScenario.cs
index 18395c4..1c0231c 100644
--- a/Scenarios/HelmholtzScenario.cs
+++ b/Scenarios/HelmholtzScenario.cs
@@ -56,6 +56,7 @@ namespace FluidSim.Tests
 
             float rho0 = 101325f / (287f * 300f);
             pipeSystem = new PipeSystem(neckCells, pipeStart, pipeEnd, areas, dxs, rho0, 0f, 101325f);
+            pipeSystem.DampingMultiplier = 500f;
 
             // --- Boundary system ---
             boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
@@ -72,7 +73,7 @@ namespace FluidSim.Tests
                 areaIndex: cavityOrificeIdx,
                 dischargeCoeff: 0.9f,
                 effectiveLength: neckLength,   // physical length (or L_eff)
-                lossCoefficient: 9000        // start with this, adjust for decay time
+                lossCoefficient: 7000        // start with this, adjust for decay time
             );
 
             // Open end at right side of pipe