diff --git a/Audio/RingBuffer.cs b/Audio/RingBuffer.cs new file mode 100644 index 0000000..efdf28c --- /dev/null +++ b/Audio/RingBuffer.cs @@ -0,0 +1,41 @@ +namespace FluidSim.Audio +{ + internal class RingBuffer + { + private readonly float[] buffer; + private volatile int readPos; + private volatile int writePos; + + public RingBuffer(int capacity) + { + if ((capacity & (capacity - 1)) != 0) + throw new ArgumentException("Capacity must be a power of two."); + buffer = new float[capacity]; + } + + public int Count => (writePos - readPos) & (buffer.Length - 1); + public int Space => (readPos - writePos - 1) & (buffer.Length - 1); + + public int Write(float[] data, int count) + { + int space = Space; + int toWrite = Math.Min(count, space); + int mask = buffer.Length - 1; + for (int i = 0; i < toWrite; i++) + buffer[(writePos + i) & mask] = data[i]; + writePos = (writePos + toWrite) & mask; + return toWrite; + } + + public int Read(float[] destination, int count) + { + int available = Count; + int toRead = Math.Min(count, available); + int mask = buffer.Length - 1; + for (int i = 0; i < toRead; i++) + destination[i] = buffer[(readPos + i) & mask]; + readPos = (readPos + toRead) & mask; + return toRead; + } + } +} diff --git a/Audio/RingBufferStream.cs b/Audio/RingBufferStream.cs new file mode 100644 index 0000000..58d3bab --- /dev/null +++ b/Audio/RingBufferStream.cs @@ -0,0 +1,43 @@ +using SFML.Audio; +using SFML.System; + +namespace FluidSim.Audio +{ + internal class RingBufferStream : SoundStream + { + private readonly RingBuffer ringBuffer; + + public RingBufferStream(RingBuffer buffer) + { + ringBuffer = buffer; + // 2 channels, 44.1 kHz, standard stereo mapping + Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight }); + } + + protected override bool OnGetData(out short[] samples) + { + const int monoBlockSize = 512; // number of mono samples we'll read + float[] temp = new float[monoBlockSize]; + int read = ringBuffer.Read(temp, monoBlockSize); + samples = new short[monoBlockSize * 2]; + + if (read > 0) + { + for (int i = 0; i < read; i++) + { + float clamped = Math.Clamp(temp[i], -1f, 1f); + short final = (short)(clamped * short.MaxValue); + samples[i * 2] = final; // left + samples[i * 2 + 1] = final; // right + } + } + for (int i = read * 2; i < samples.Length; i++) + samples[i] = 0; + + return true; + } + + protected override void OnSeek(Time timeOffset) => + throw new NotSupportedException(); + } +} diff --git a/Audio/SoundEngine.cs b/Audio/SoundEngine.cs new file mode 100644 index 0000000..c8cfb9d --- /dev/null +++ b/Audio/SoundEngine.cs @@ -0,0 +1,45 @@ +namespace FluidSim.Audio; + +public class SoundEngine : IDisposable +{ + private readonly RingBuffer ringBuffer; + private readonly RingBufferStream stream; + private bool isPlaying; + + public SoundEngine(int bufferCapacity = 16384) + { + ringBuffer = new RingBuffer(bufferCapacity); + stream = new RingBufferStream(ringBuffer); + } + + public void Start() + { + if (isPlaying) return; + stream.Play(); + isPlaying = true; + } + + public void Stop() + { + if (!isPlaying) return; + stream.Stop(); + isPlaying = false; + float[] drain = new float[ringBuffer.Count]; + ringBuffer.Read(drain, drain.Length); + } + + public int WriteSamples(float[] data, int count) => + ringBuffer.Write(data, count); + + public float Volume + { + get => stream.Volume; + set => stream.Volume = value; + } + + public void Dispose() + { + Stop(); + stream.Dispose(); + } +} \ No newline at end of file diff --git a/Audio/SoundProcessor.cs b/Audio/SoundProcessor.cs new file mode 100644 index 0000000..18b64d6 --- /dev/null +++ b/Audio/SoundProcessor.cs @@ -0,0 +1,29 @@ +using FluidSim.Components; +using FluidSim.Interfaces; +using System; + +namespace FluidSim.Audio +{ + public static class SoundProcessor + { + public static float MaxDeltaP { get; set; } = 100_000f; + public static float MaxArea { get; set; } = 1e-4f; + public static float MaxVelocity { get; set; } = 343f; + public static float ReferenceDensity { get; set; } = 1.225f; + public static float ReferenceSpeedOfSound { get; set; } = 343f; + public static float Gain { get; set; } = 1.0f; + + public static double ComputeSample(Connection conn) + { + Port portA = conn.PortA; + Port portB = conn.PortB; + + double pressureUp = portA.Pressure; + double pressureDown = portB.Pressure; + + // No flow or no pressure drop → silence + double deltaP = pressureUp - pressureDown; + return deltaP / 1; + } + } +} \ No newline at end of file diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 5b6e2aa..2ff7e6d 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -12,29 +12,52 @@ namespace FluidSim.Components private int _n; private double _dx, _dt, _gamma = 1.4, _area; private double[] _rho, _rhou, _E; + private double _hydraulicDiameter; - // Volume states at boundaries private double _rhoLeft, _pLeft, _rhoRight, _pRight; private bool _leftBCSet, _rightBCSet; - public double FrictionFactor { get; set; } = 0.02; + public double FrictionFactor { get; set; } public int GetCellCount() => _n; public double GetCellDensity(int i) => _rho[i]; public double GetCellPressure(int i) => Pressure(i); public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); - public Pipe1D(double length, double area, int nCells, int sampleRate) + /// + /// Create a pipe with CFL‑stable automatic cell count. + /// + /// Pipe length [m]. + /// Cross‑sectional area [m²]. + /// Simulation step rate [Hz]. + /// Speed of sound [m/s] (default 343). + /// Darcy friction factor (0 = inviscid). + /// CFL safety factor ≤ 1 (0.8 recommended). + public Pipe1D(double length, double area, int sampleRate, + double c0 = 343.0, double frictionFactor = 0.02, + double cflSafety = 0.8) { - _n = nCells; - _dx = length / nCells; - _dt = 1.0 / sampleRate; + if (area <= 0) throw new ArgumentException("Pipe area must be > 0"); _area = area; + _dt = 1.0 / sampleRate; + FrictionFactor = frictionFactor; + + // Nyquist‑based cell count (wave resolution) + double nNyquist = Math.Ceiling(length * sampleRate / c0); + // CFL‑stable cell count: dx ≥ maxSpeed·dt / cflSafety, maxSpeed = 2·c0 (supersonic safe) + double maxSpeed = 2.0 * c0; + double dxMinStable = maxSpeed * _dt / cflSafety; + double nStable = Math.Floor(length / dxMinStable); + + _n = Math.Max(2, (int)Math.Min(nNyquist, nStable)); + _dx = length / _n; _rho = new double[_n]; _rhou = new double[_n]; _E = new double[_n]; + _hydraulicDiameter = Math.Max(2.0 * Math.Sqrt(_area / Math.PI), 1e-9); + PortA = new Port(); PortB = new Port(); } @@ -56,7 +79,6 @@ namespace FluidSim.Components public double GetLeftDensity() => _rho[0]; public double GetRightDensity() => _rho[_n - 1]; - // ★ New: pass both density and pressure from the volume public void SetLeftVolumeState(double rhoVol, double pVol) { _rhoLeft = rhoVol; @@ -88,22 +110,13 @@ namespace FluidSim.Components // --- Left boundary (face 0) --- if (_leftBCSet) { - // Ghost = actual volume state (ρ_vol, u=0, p_vol) - double rhoL = _rhoLeft; - double uL = 0.0; - double pL = _pLeft; - - double rhoR = _rho[0]; - double uR = _rhou[0] / Math.Max(rhoR, 1e-12); - double pR = Pressure(0); - + double rhoL = _rhoLeft, uL = 0.0, pL = _pLeft; + double rhoR = _rho[0], uR = _rhou[0] / Math.Max(rhoR, 1e-12), pR = Pressure(0); HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]); } else { - Fm[0] = 0; - Fp[0] = Pressure(0); - Fe[0] = 0; + Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; } // --- Internal faces --- @@ -118,25 +131,16 @@ namespace FluidSim.Components // --- Right boundary (face n) --- if (_rightBCSet) { - double rhoL = _rho[n - 1]; - double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12); - double pL = Pressure(n - 1); - - // Ghost = actual volume state (ρ_vol, u=0, p_vol) - double rhoR = _rhoRight; - double uR = 0.0; - double pR = _pRight; - + double rhoL = _rho[n - 1], uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12), pL = Pressure(n - 1); + double rhoR = _rhoRight, uR = 0.0, pR = _pRight; HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]); } else { - Fm[n] = 0; - Fp[n] = Pressure(n - 1); - Fe[n] = 0; + Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; } - // --- Cell update --- + // --- Cell update (inviscid fluxes) --- for (int i = 0; i < n; i++) { double dM = (Fm[i + 1] - Fm[i]) / _dx; @@ -149,12 +153,41 @@ namespace FluidSim.Components if (_rho[i] < 1e-12) _rho[i] = 1e-12; double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i]; if (_E[i] < kinetic) _E[i] = kinetic; + + // Emergency reset if NaN + if (double.IsNaN(_rho[i]) || double.IsNaN(_rhou[i]) || double.IsNaN(_E[i])) + { + _rho[i] = 1.225; // reset to atmospheric air at 300 K + _rhou[i] = 0.0; + _E[i] = 101325.0 / (_gamma - 1.0); // internal energy at 1 atm + } } - // --- Friction disabled --- - // if (FrictionFactor > 0) { … } + // --- Friction (Darcy–Weisbach, energy‑conserving) --- + if (FrictionFactor > 0) + { + double D = _hydraulicDiameter; + double twoD = 2.0 * D; + for (int i = 0; i < n; i++) + { + double rho = _rho[i]; + double u = _rhou[i] / rho; + double absU = Math.Abs(u); + double src = FrictionFactor * rho * absU * u / twoD; + double kinOld = 0.5 * rho * u * u; + _rhou[i] -= _dt * src; + double uNew = _rhou[i] / rho; + double kinNew = 0.5 * rho * uNew * uNew; + _E[i] += (kinOld - kinNew); + } + } + + // --- Publish to ports --- + PortA.Pressure = Pressure(0); + PortA.Density = _rho[0]; + PortB.Pressure = Pressure(_n - 1); + PortB.Density = _rho[_n - 1]; - // --- Port flows --- PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0; PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0; @@ -170,14 +203,24 @@ namespace FluidSim.Components void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, out double fm, out double fp, out double fe) { - double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12)); - double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12)); + const double eps = 1e-12; + pL = Math.Max(pL, eps); + pR = Math.Max(pR, eps); + + double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, eps)); + double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, eps)); double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL; double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR; + double SL = Math.Min(uL - cL, uR - cR); double SR = Math.Max(uL + cL, uR + cR); - double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) - / (rL * (SL - uL) - rR * (SR - uR)); + + double denom = rL * (SL - uL) - rR * (SR - uR); + double Ss; + if (Math.Abs(denom) < eps) + Ss = 0.5 * (uL + uR); + else + Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom; double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL; double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR; @@ -186,16 +229,20 @@ namespace FluidSim.Components else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; } else if (Ss >= 0) { - double rsL = rL * (SL - uL) / (SL - Ss); - double ps = pL + rL * (SL - uL) * (Ss - uL); - double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL))); + double diffSL = SL - uL; + if (Math.Abs(diffSL) < eps) diffSL = eps; + double rsL = rL * diffSL / (SL - Ss); + double ps = pL + rL * diffSL * (Ss - uL); + double EsL = EL + (Ss - uL) * (Ss + pL / (rL * diffSL)); fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss; } else { - double rsR = rR * (SR - uR) / (SR - Ss); + double diffSR = SR - uR; + if (Math.Abs(diffSR) < eps) diffSR = eps; + double rsR = rR * diffSR / (SR - Ss); double ps = pL + rL * (SL - uL) * (Ss - uL); - double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR))); + double EsR = ER + (Ss - uR) * (Ss + pR / (rR * diffSR)); fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss; } } diff --git a/Components/SoundConnection.cs b/Components/SoundConnection.cs new file mode 100644 index 0000000..3ebe259 --- /dev/null +++ b/Components/SoundConnection.cs @@ -0,0 +1,9 @@ +using FluidSim.Interfaces; + +namespace FluidSim.Components +{ + public class SoundConnection : Connection + { + public SoundConnection(Port a, Port b) : base(a, b) { } + } +} \ No newline at end of file diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs index b5949a2..d380288 100644 --- a/Components/Volume0D.cs +++ b/Components/Volume0D.cs @@ -1,6 +1,4 @@ -using System; -using FluidSim.Interfaces; -using FluidSim.Utils; +using FluidSim.Interfaces; namespace FluidSim.Components { diff --git a/Core/OrificeBoundary.cs b/Core/OrificeBoundary.cs index 733eba2..87c1e85 100644 --- a/Core/OrificeBoundary.cs +++ b/Core/OrificeBoundary.cs @@ -1,5 +1,4 @@ -using System; -using FluidSim.Components; +using FluidSim.Components; namespace FluidSim.Core { diff --git a/Core/Simulation.cs b/Core/Simulation.cs deleted file mode 100644 index 2b310c5..0000000 --- a/Core/Simulation.cs +++ /dev/null @@ -1,67 +0,0 @@ -using System; -using FluidSim.Components; -using FluidSim.Utils; - -namespace FluidSim.Core -{ - public static class Simulation - { - private static Solver solver; - private static Volume0D volA, volB; - private static Pipe1D pipe; - private static Connection connA, connB; - private static int stepCount; - private static double time; - private static double dt; - - public static void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - - double V = 5.0 * Units.L; - volA = new Volume0D(V, 1.1 * Units.atm, Units.Celsius(20), sampleRate); - volB = new Volume0D(V, 1.0 * Units.atm, Units.Celsius(20), sampleRate); - - double length = 150 * Units.mm; - double diameter = 25 * Units.mm; - double area = Units.AreaFromDiameter(25, Units.mm); - int nCells = 10; - - pipe = new Pipe1D(length, area, nCells, sampleRate); - pipe.SetUniformState(volA.Density, 0.0, volA.Pressure); - pipe.FrictionFactor = 0.02; - - // Connections with orifice area equal to pipe area (flange joint) - connA = new Connection(volA.Port, pipe.PortA) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; - connB = new Connection(pipe.PortB, volB.Port) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; - - solver = new Solver(); - solver.AddVolume(volA); - solver.AddVolume(volB); - solver.AddPipe(pipe); - solver.AddConnection(connA); - solver.AddConnection(connB); - } - - public static float Process() - { - solver.Step(); - time += dt; - stepCount++; - Log(); - return 0f; - } - - public static void Log() - { - if (stepCount <= 50 || stepCount % 200 == 0) - { - Console.WriteLine( - $"t = {time * 1e3:F3} ms Step {stepCount:D4}: " + - $"PA = {volA.Pressure / 1e5:F6} bar, " + - $"PB = {volB.Pressure / 1e5:F6} bar, " + - $"FlowA = {pipe.PortA.MassFlowRate * 1e3:F2} g/s"); - } - } - } -} \ No newline at end of file diff --git a/Core/Solver.cs b/Core/Solver.cs index 0c7807f..6e74ab9 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -1,4 +1,4 @@ -using System.Collections.Generic; +using FluidSim.Audio; using FluidSim.Components; using FluidSim.Interfaces; @@ -10,17 +10,51 @@ namespace FluidSim.Core private readonly List _pipes = new(); private readonly List _connections = new(); + public float LastSample { get; private set; } + public void AddVolume(Volume0D v) => _volumes.Add(v); public void AddPipe(Pipe1D p) => _pipes.Add(p); public void AddConnection(Connection c) => _connections.Add(c); public void Step() { - // 1. Volumes publish state to their ports + // 1. Publish volume states to their own ports foreach (var v in _volumes) v.PushStateToPort(); - // 2. Set volume states as boundary conditions on pipes + // 2. Handle direct volume‑to‑volume connections + foreach (var conn in _connections) + { + if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB)) + { + Volume0D volA = _volumes.Find(v => v.Port == conn.PortA); + Volume0D volB = _volumes.Find(v => v.Port == conn.PortB); + + if (volA == null || volB == null) continue; + + double pA = volA.Pressure, rhoA = volA.Density; + double pB = volB.Pressure, rhoB = volB.Density; + + double mdot = OrificeBoundary.MassFlow(pA, rhoA, pB, rhoB, conn); + + if (mdot > 0) // A → B + { + volA.Port.MassFlowRate = -mdot; + volB.Port.MassFlowRate = mdot; + volB.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // fluid from A + volA.Port.SpecificEnthalpy = volA.SpecificEnthalpy; // outflow carries its own enthalpy + } + else // B → A (mdot negative) + { + volA.Port.MassFlowRate = -mdot; // positive + volB.Port.MassFlowRate = mdot; // negative + volA.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // fluid from B + volB.Port.SpecificEnthalpy = volB.SpecificEnthalpy; // outflow carries its own + } + } + } + + // 3. Pipe‑volume boundary conditions – unchanged foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) @@ -29,11 +63,11 @@ namespace FluidSim.Core SetVolumeBC(conn.PortB, conn.PortA); } - // 3. Run pipe simulations + // 4. Run pipe simulations foreach (var p in _pipes) p.Simulate(); - // 4. Transfer pipe‑port flows to volume ports + // 5. Transfer pipe‑to‑volume flows – unchanged foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) @@ -42,9 +76,21 @@ namespace FluidSim.Core TransferPipeToVolume(conn.PortB, conn.PortA); } - // 5. Integrate volumes + // 6. Integrate volumes foreach (var v in _volumes) v.Integrate(); + + // 7. COMPUTE AUDIO SAMPLE from all sound connections + double sample = 0f; + foreach (var conn in _connections) + { + if (conn is SoundConnection) + { + // Both ports have the same absolute mass flow; either works. + sample += SoundProcessor.ComputeSample(conn); + } + } + LastSample = (float)Math.Tanh(sample); } bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p); diff --git a/Core/SoundEngine.cs b/Core/SoundEngine.cs deleted file mode 100644 index bab240a..0000000 --- a/Core/SoundEngine.cs +++ /dev/null @@ -1,131 +0,0 @@ -using SFML.Audio; -using SFML.System; - -namespace FluidSim; - -#region Lock‑free ring buffer (unchanged) -internal class RingBuffer -{ - private readonly float[] buffer; - private volatile int readPos; - private volatile int writePos; - - public RingBuffer(int capacity) - { - if ((capacity & (capacity - 1)) != 0) - throw new ArgumentException("Capacity must be a power of two."); - buffer = new float[capacity]; - } - - public int Count => (writePos - readPos) & (buffer.Length - 1); - public int Space => (readPos - writePos - 1) & (buffer.Length - 1); - - public int Write(float[] data, int count) - { - int space = Space; - int toWrite = Math.Min(count, space); - int mask = buffer.Length - 1; - for (int i = 0; i < toWrite; i++) - buffer[(writePos + i) & mask] = data[i]; - writePos = (writePos + toWrite) & mask; - return toWrite; - } - - public int Read(float[] destination, int count) - { - int available = Count; - int toRead = Math.Min(count, available); - int mask = buffer.Length - 1; - for (int i = 0; i < toRead; i++) - destination[i] = buffer[(readPos + i) & mask]; - readPos = (readPos + toRead) & mask; - return toRead; - } -} -#endregion - -#region Stereo stream that consumes the ring buffer -internal class RingBufferStream : SoundStream -{ - private readonly RingBuffer ringBuffer; - - public RingBufferStream(RingBuffer buffer) - { - ringBuffer = buffer; - // 2 channels, 44.1 kHz, standard stereo mapping - Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight }); - } - - protected override bool OnGetData(out short[] samples) - { - const int monoBlockSize = 512; // number of mono samples we'll read - float[] temp = new float[monoBlockSize]; - int read = ringBuffer.Read(temp, monoBlockSize); - samples = new short[monoBlockSize * 2]; - - if (read > 0) - { - for (int i = 0; i < read; i++) - { - float clamped = Math.Clamp(temp[i], -1f, 1f); - short final = (short)(clamped * short.MaxValue); - samples[i * 2] = final; // left - samples[i * 2 + 1] = final; // right - } - } - for (int i = read * 2; i < samples.Length; i++) - samples[i] = 0; - - return true; - } - - protected override void OnSeek(Time timeOffset) => - throw new NotSupportedException(); -} -#endregion - -#region Public sound engine API (unchanged) -public class SoundEngine : IDisposable -{ - private readonly RingBuffer ringBuffer; - private readonly RingBufferStream stream; - private bool isPlaying; - - public SoundEngine(int bufferCapacity = 16384) - { - ringBuffer = new RingBuffer(bufferCapacity); - stream = new RingBufferStream(ringBuffer); - } - - public void Start() - { - if (isPlaying) return; - stream.Play(); - isPlaying = true; - } - - public void Stop() - { - if (!isPlaying) return; - stream.Stop(); - isPlaying = false; - float[] drain = new float[ringBuffer.Count]; - ringBuffer.Read(drain, drain.Length); - } - - public int WriteSamples(float[] data, int count) => - ringBuffer.Write(data, count); - - public float Volume - { - get => stream.Volume; - set => stream.Volume = value; - } - - public void Dispose() - { - Stop(); - stream.Dispose(); - } -} -#endregion \ No newline at end of file diff --git a/Program.cs b/Program.cs index 938ced0..a07b687 100644 --- a/Program.cs +++ b/Program.cs @@ -2,7 +2,8 @@ using SFML.Window; using SFML.System; using System.Diagnostics; -using FluidSim.Core; +using FluidSim.Scenarios; +using FluidSim.Audio; namespace FluidSim; @@ -10,6 +11,8 @@ public class Program { private const int SampleRate = 44100; private static volatile bool running = true; + // Global step counter – incremented every simulation step + private static long stepCount = 0; public static void Main() { @@ -22,47 +25,71 @@ public class Program soundEngine.Volume = 70; soundEngine.Start(); - double lastAudioTime = 0.0; var stopwatch = Stopwatch.StartNew(); - int warmupSamples = SampleRate / 2; + // --- Warmup: fill audio buffer with silence --- + int warmupSamples = SampleRate / 2; // 0.5 s float[] warmup = new float[warmupSamples]; for (int i = 0; i < warmupSamples; i++) warmup[i] = 0; - soundEngine.WriteSamples(warmup, warmupSamples); - lastAudioTime = stopwatch.Elapsed.TotalSeconds; + + // Reset timer after warmup – this is the “real‑time zero” + stopwatch.Restart(); + stepCount = 0; // simulation steps start now + + // --- Initialise the simulation scenario --- + Simulation.Initialize(SampleRate); const int chunkSize = 2048; float[] buffer = new float[chunkSize]; - Simulation.Initialize(SampleRate); + double lastLogTime = 0.0; // for periodic speed printout while (window.IsOpen) { window.DispatchEvents(); + // --- Compute how many audio samples are needed since last frame --- double currentTime = stopwatch.Elapsed.TotalSeconds; - double elapsed = currentTime - lastAudioTime; - int samplesNeeded = (int)(elapsed * SampleRate); + double elapsed = currentTime; // since stopwatch was reset + int samplesNeeded = (int)(elapsed * SampleRate) - (int)(stepCount); + // (stepCount is total generated samples, so we just need the remainder) + // --- Generate the required number of simulation steps --- while (samplesNeeded > 0 && running) { int toGenerate = Math.Min(samplesNeeded, chunkSize); for (int i = 0; i < toGenerate; i++) { buffer[i] = Simulation.Process(); + stepCount++; } soundEngine.WriteSamples(buffer, toGenerate); samplesNeeded -= toGenerate; } - lastAudioTime = currentTime; + // --- Display speed --- + double simTime = stepCount / (double)SampleRate; + double wallTime = stopwatch.Elapsed.TotalSeconds; + double speed = (wallTime > 0) ? simTime / wallTime : 0.0; + // Update window title with instant speed + window.SetTitle($"FluidSim | Speed: {speed:F3}× | Steps: {stepCount}"); + + // Console log once per second + if (wallTime - lastLogTime >= 1.0) + { + Console.WriteLine($"Speed: {speed:F3}× ({stepCount} steps, {wallTime:F2}s wall)"); + lastLogTime = wallTime; + } + + // --- Rendering (placeholder) --- window.Clear(Color.Black); window.Display(); } + // --- Cleanup --- soundEngine.Dispose(); window.Dispose(); } diff --git a/Scenarios/Simulation.cs b/Scenarios/Simulation.cs new file mode 100644 index 0000000..80dc89d --- /dev/null +++ b/Scenarios/Simulation.cs @@ -0,0 +1,93 @@ +using FluidSim.Audio; +using FluidSim.Components; +using FluidSim.Core; +using FluidSim.Utilities; + +namespace FluidSim.Scenarios +{ + public static class Simulation + { + private static Solver solver; + private static Volume0D cavity, ambient; + private static Pipe1D neck; + private static Connection connNeckCavity, connNeckAmbient; + private static int stepCount; + private static double time; + private static double dt; + + public static void Initialize(int sampleRate) + { + dt = 1.0 / sampleRate; + + // --- Cavity (the “bottle”) --- + double V_cav = 1.0 * Units.L; // 1 litre + cavity = new Volume0D(V_cav, 10 * Units.atm, Units.Celsius(20), sampleRate); + + // --- Ambient (a huge “room”) --- + double V_amb = 1000.0; // 1000 m³ ≈ constant pressure + ambient = new Volume0D(V_amb, 1.0 * Units.atm, Units.Celsius(20), sampleRate); + + // --- Neck (short pipe) --- + double length_neck = 80.0 * Units.mm; + double diam_neck = 10 * Units.mm; + double area_neck = Units.AreaFromRadius(diam_neck, Units.mm); // few cells enough for a short neck + neck = new Pipe1D(length_neck, area_neck, sampleRate); + neck.SetUniformState(ambient.Density, 0.0, ambient.Pressure); + neck.FrictionFactor = 0.02; // slight damping + + // --- Connections --- + // Cavity-to-neck (full area) + connNeckCavity = new Connection(cavity.Port, neck.PortA) + { + Area = area_neck, + DischargeCoefficient = 1.0, + Gamma = 1.4 + }; + // Neck-to-ambient (SoundConnection to capture the radiated tone) + connNeckAmbient = new SoundConnection(neck.PortB, ambient.Port) + { + Area = area_neck, + DischargeCoefficient = 1.0, + Gamma = 1.4 + }; + + // --- Solver --- + solver = new Solver(); + solver.AddVolume(cavity); + solver.AddVolume(ambient); + solver.AddPipe(neck); + solver.AddConnection(connNeckCavity); + solver.AddConnection(connNeckAmbient); + + // --- Sound tuning --- + SoundProcessor.MaxDeltaP = 0.1f * (float)Units.atm; // small Δp expected + SoundProcessor.MaxArea = (float)area_neck; + SoundProcessor.MaxVelocity = 343f; + SoundProcessor.ReferenceDensity = 1.2f; + SoundProcessor.ReferenceSpeedOfSound = 343f; + SoundProcessor.Gain = 10.0f; // amplify because Δp is small + } + + public static float Process() + { + solver.Step(); + time += dt; + stepCount++; + Log(); + return solver.LastSample; + } + + public static void Log() + { + if (stepCount <= 200 || stepCount % (int)(0.5 / dt) == 0) + { + Console.WriteLine( + $"t = {time * 1e3:F3} ms " + + $"Sample = {solver.LastSample:F6}, " + + $"P_cav = {cavity.Pressure / 1e5:F6} bar, " + + $"flow_cav = {cavity.Port.MassFlowRate / 1e3:F4} g/s, " + + $"flow_neck = {neck.PortB.MassFlowRate * 1e3:F4} g/s"); + } + } + } +} \ No newline at end of file diff --git a/Sources/EffortSource.cs b/Sources/EffortSource.cs deleted file mode 100644 index 78e186c..0000000 --- a/Sources/EffortSource.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Sources -{ - internal class EffortSource - { - } -} diff --git a/Sources/FlowSource.cs b/Sources/FlowSource.cs deleted file mode 100644 index 2004336..0000000 --- a/Sources/FlowSource.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Sources -{ - internal class FlowSource - { - } -} diff --git a/Utils/Units.cs b/Utilities/Units.cs similarity index 95% rename from Utils/Units.cs rename to Utilities/Units.cs index b151478..bade8e6 100644 --- a/Utils/Units.cs +++ b/Utilities/Units.cs @@ -1,6 +1,4 @@ -using System; - -namespace FluidSim.Utils +namespace FluidSim.Utilities { public static class Units {