grillkol/FluidSim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: grillkol:master
Branches Tags
grillkol:master grillkol:Testing grillkol:testing
..
compare: grillkol:testing
Branches Tags
grillkol:Testing grillkol:master grillkol:testing

1 Commits
master ... testing

Author SHA1 Message Date
max
a262410616 Attempting to generate sound, and set cells automatically 2026-05-02 23:26:52 +02:00
22 changed files with 558 additions and 1314 deletions
Expand all files Collapse all files

41
Audio/RingBuffer.cs Normal file
View File

@@ -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;
}
}
}

43
Audio/RingBufferStream.cs Normal file
View File

@@ -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();
}
}

45
Audio/SoundEngine.cs Normal file
View File

@@ -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();
}
}

29
Audio/SoundProcessor.cs Normal file
View File

@@ -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;
}
}
}

4
Interfaces/Connection.cs → Components/Connection.cs
View File

@@ -1,4 +1,6 @@
namespace FluidSim.Interfaces
using FluidSim.Interfaces;
namespace FluidSim.Components
{
/// <summary>Pure data link between two ports, with orifice parameters.</summary>
public class Connection

494
Components/Pipe1D.cs
View File

@@ -3,85 +3,65 @@ using FluidSim.Interfaces;
namespace FluidSim.Components
{
public enum BoundaryType
{
VolumeCoupling,
OpenEnd,
ClosedEnd
}
public class Pipe1D
{
public Port PortA { get; }
public Port PortB { get; }
public double Area => _area;
public double DampingMultiplier { get; set; } = 1.0;
private int _n;
private double _dx, _dt, _gamma, _area, _diameter;
private double _dx, _dt, _gamma = 1.4, _area;
private double[] _rho, _rhou, _E;
private double _hydraulicDiameter;
// Volume‑coupling ghost states for boundaries A and B
private double _rhoA, _pA;
private double _rhoB, _pB;
private bool _aBCSet, _bBCSet;
private double _rhoLeft, _pLeft, _rhoRight, _pRight;
private bool _leftBCSet, _rightBCSet;
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 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 BoundaryType ABCType => _aBCType;
public BoundaryType BBCType => _bBCType;
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
/// <summary>
/// Create a pipe with CFL‑stable automatic cell count.
/// </summary>
/// <param name="length">Pipe length [m].</param>
/// <param name="area">Cross‑sectional area [m²].</param>
/// <param name="sampleRate">Simulation step rate [Hz].</param>
/// <param name="c0">Speed of sound [m/s] (default 343).</param>
/// <param name="frictionFactor">Darcy friction factor (0 = inviscid).</param>
/// <param name="cflSafety">CFL safety factor ≤ 1 (0.8 recommended).</param>
public Pipe1D(double length, double area, int sampleRate,
double c0 = 343.0, double frictionFactor = 0.02,
double cflSafety = 0.8)
{
double dtGlobal = 1.0 / sampleRate;
int nCells;
if (forcedCellCount > 1)
{
nCells = forcedCellCount;
}
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;
_dx = length / _n;
_dt = dtGlobal;
if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
_area = area;
_gamma = 1.4;
_dt = 1.0 / sampleRate;
FrictionFactor = frictionFactor;
// Hydraulic diameter for a circular pipe
_diameter = 2.0 * Math.Sqrt(area / Math.PI);
// 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();
}
public void SetABoundaryType(BoundaryType type) => _aBCType = type;
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)
{
double e = p / ((_gamma - 1) * rho);
@@ -94,181 +74,23 @@ namespace FluidSim.Components
}
}
public void SetCellState(int i, double rho, double u, double p)
public double GetLeftPressure() => Pressure(0);
public double GetRightPressure() => Pressure(_n - 1);
public double GetLeftDensity() => _rho[0];
public double GetRightDensity() => _rho[_n - 1];
public void SetLeftVolumeState(double rhoVol, double pVol)
{
if (i < 0 || i >= _n) return;
_rho[i] = rho;
_rhou[i] = rho * u;
double e = p / ((_gamma - 1) * rho);
_E[i] = rho * e + 0.5 * rho * u * u;
_rhoLeft = rhoVol;
_pLeft = pVol;
_leftBCSet = true;
}
public void SetAVolumeState(double rhoVol, double pVol)
public void SetRightVolumeState(double rhoVol, double pVol)
{
_rhoA = rhoVol;
_pA = pVol;
_aBCSet = true;
}
public void SetBVolumeState(double rhoVol, double pVol)
{
_rhoB = rhoVol;
_pB = pVol;
_bBCSet = true;
}
public void ClearBC() => _aBCSet = _bBCSet = false;
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 u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
double local = u + c;
if (local > maxW) maxW = local;
}
maxW = Math.Max(maxW, 1e-8);
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
}
public void SimulateSingleStep(double 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) ----------
double rhoIntA = _rho[0];
double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
double pIntA = Pressure(0);
switch (_aBCType)
{
case BoundaryType.VolumeCoupling:
if (_aBCSet)
{
HLLCFlux(_rhoA, 0.0, _pA,
rhoIntA, uIntA, pIntA,
out Fm[0], out Fp[0], out Fe[0]);
}
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 ----------
for (int i = 0; i < n - 1; i++)
{
double rhoL = _rho[i];
double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
double pL = Pressure(i);
double rhoR = _rho[i + 1];
double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
double pR = Pressure(i + 1);
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 dP = (Fp[i + 1] - Fp[i]) / _dx;
double dE = (Fe[i + 1] - Fe[i]) / _dx;
_rho[i] -= dtSub * dM;
_rhou[i] -= dtSub * dP;
_E[i] -= dtSub * dE;
double rho = Math.Max(_rho[i], 1e-12);
double dampingRate = laminarCoeff / rho;
double dampingFactor = Math.Exp(-dampingRate * dtSub);
_rhou[i] *= dampingFactor;
if (_rho[i] < 1e-12) _rho[i] = 1e-12;
double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
double pMin = 100.0;
double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
if (_E[i] < eMin) _E[i] = eMin;
}
// ---------- Port quantities ----------
double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
PortA.MassFlowRate = mdotA_sub;
PortB.MassFlowRate = mdotB_sub;
PortA.Pressure = pIntA;
PortB.Pressure = pIntB;
PortA.Density = _rho[0];
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);
}
_rhoRight = rhoVol;
_pRight = pVol;
_rightBCSet = true;
}
private double GetCellTotalSpecificEnthalpy(int i)
@@ -280,114 +102,125 @@ namespace FluidSim.Components
return h + 0.5 * u * u;
}
private double Pressure(int i) =>
public void Simulate()
{
int n = _n;
double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];
// --- Left boundary (face 0) ---
if (_leftBCSet)
{
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;
}
// --- Internal faces ---
for (int i = 0; i < n - 1; i++)
{
double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
HLLCFlux(_rho[i], uL, Pressure(i), _rho[i + 1], uR, Pressure(i + 1),
out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
}
// --- Right boundary (face n) ---
if (_rightBCSet)
{
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;
}
// --- Cell update (inviscid fluxes) ---
for (int i = 0; i < n; i++)
{
double dM = (Fm[i + 1] - Fm[i]) / _dx;
double dP = (Fp[i + 1] - Fp[i]) / _dx;
double dE = (Fe[i + 1] - Fe[i]) / _dx;
_rho[i] -= _dt * dM;
_rhou[i] -= _dt * dP;
_E[i] -= _dt * dE;
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 (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];
PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0;
PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0;
PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0);
PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1);
_leftBCSet = _rightBCSet = false;
}
double Pressure(int i) =>
(_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
// ========== Characteristic‑based Open End ==========
private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
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));
const double eps = 1e-12;
pL = Math.Max(pL, eps);
pR = Math.Max(pR, eps);
// Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
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 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;
@@ -396,25 +229,22 @@ 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;
}
}
public double GetPressureAtFraction(double fraction)
{
int i = (int)(fraction * (_n - 1));
i = Math.Clamp(i, 0, _n - 1);
return Pressure(i);
}
}
}

9
Components/SoundConnection.cs Normal file
View File

@@ -0,0 +1,9 @@
using FluidSim.Interfaces;
namespace FluidSim.Components
{
public class SoundConnection : Connection
{
public SoundConnection(Port a, Port b) : base(a, b) { }
}
}

21
Components/Volume0D.cs
View File

@@ -1,6 +1,4 @@
using System;
using FluidSim.Interfaces;
using FluidSim.Utils;
using FluidSim.Interfaces;
namespace FluidSim.Components
{
@@ -50,26 +48,13 @@ namespace FluidSim.Components
Port.SpecificEnthalpy = SpecificEnthalpy;
}
// Original integrate (uses the constructor’s sample rate)
public void Integrate()
{
Integrate(_dt);
}
public void SetPressure(double newPressure)
{
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;
double dm = mdot * _dt;
double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
Mass += dm;
InternalEnergy += dE;

3
Core/OrificeBoundary.cs
View File

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

179
Core/Solver.cs
View File

@@ -1,5 +1,4 @@
using System;
using System.Collections.Generic;
using FluidSim.Audio;
using FluidSim.Components;
using FluidSim.Interfaces;
@@ -11,160 +10,116 @@ namespace FluidSim.Core
private readonly List<Pipe1D> _pipes = new();
private readonly List<Connection> _connections = new();
private double _dt;
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 SetTimeStep(double dt) => _dt = dt;
/// <summary>
/// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right).
/// </summary>
public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0)
public void Step()
{
if (isA)
{
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. Volumes publish state
// 1. Publish volume states to their own ports
foreach (var v in _volumes)
v.PushStateToPort();
// 2. Set volume BCs for volume‑coupled ends
// 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))
{
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;
// 4. Run pipe simulations
foreach (var p in _pipes)
nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt));
double dtSub = _dt / nSub;
for (int sub = 0; sub < nSub; sub++)
{
foreach (var p in _pipes)
p.SimulateSingleStep(dtSub);
p.Simulate();
// 5. Transfer pipe‑to‑volume flows – unchanged
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);
}
TransferPipeToVolume(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))
TransferAndIntegrate(conn.PortB, conn.PortA, dtSub);
}
TransferPipeToVolume(conn.PortB, conn.PortA);
}
if (sub < nSub - 1)
{
// 6. Integrate volumes
foreach (var v in _volumes)
v.PushStateToPort();
v.Integrate();
// 7. COMPUTE AUDIO SAMPLE from all sound connections
double sample = 0f;
foreach (var conn in _connections)
{
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
if (conn is SoundConnection)
{
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);
// Both ports have the same absolute mass flow; either works.
sample += SoundProcessor.ComputeSample(conn);
}
else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
}
LastSample = (float)Math.Tanh(sample);
}
bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p);
bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
void SetVolumeBC(Port pipePort, Port volPort)
{
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);
Pipe1D pipe = GetPipe(pipePort);
if (pipe == null) return;
bool isA = pipe.PortA == pipePort;
if (isA)
pipe.SetAVolumeState(volPort.Density, volPort.Pressure);
bool isLeft = pipe.PortA == pipePort;
if (isLeft)
pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
else
pipe.SetBVolumeState(volPort.Density, volPort.Pressure);
pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
}
private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub)
void TransferPipeToVolume(Port pipePort, Port volPort)
{
double mdot = pipePort.MassFlowRate;
volPort.MassFlowRate = -mdot;
if (mdot < 0) // pipe → volume
{
// pipePort.SpecificEnthalpy is already total (h + ½u²)
volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
}
// else volume’s own enthalpy (from PushStateToPort) is used
GetVolume(volPort)?.Integrate(dtSub);
// else: volume → pipe, volume’s own static enthalpy is used (already set)
}
}
}

131
Core/SoundEngine.cs
View File

@@ -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

23
Core/SoundProcessor.cs
View File

@@ -1,23 +0,0 @@
namespace FluidSim.Core
{
/// <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>
/// 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;
foreach (float s in samples)
sum += s;
sum *= MasterGain;
return sum;
}
}
}

25
Interfaces/SoundConnection.cs
View File

@@ -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;
}
}
}

182
Program.cs
View File

@@ -2,183 +2,95 @@
using SFML.Window;
using SFML.System;
using System.Diagnostics;
using FluidSim.Core;
using FluidSim.Scenarios;
using FluidSim.Audio;
namespace FluidSim;
public class Program
{
private const int SampleRate = 44100;
private const double DrawFrequency = 60.0;
private static Scenario scenario;
// Speed control
//private static double desiredSpeed = 1.0;
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; // remembers the last non‑1.0 scroll speed
private static bool isRealTime = true; // true when desiredSpeed == 1.0
private static volatile bool running = true;
// Global step counter – incremented every simulation step
private static long stepCount = 0;
public static void Main()
{
var mode = new VideoMode(new Vector2u(1280, 720));
var window = new RenderWindow(mode, "Pipe Resonator");
var window = new RenderWindow(mode, "Fluid Simulation");
window.SetVerticalSyncEnabled(true);
window.Closed += (_, _) => { running = false; window.Close(); };
window.MouseWheelScrolled += OnMouseWheel;
window.KeyPressed += OnKeyPressed;
var soundEngine = new SoundEngine(bufferCapacity: 16384);
var soundEngine = new SoundEngine(bufferCapacity: 2048);
soundEngine.Volume = 70;
soundEngine.Start();
//scenario = new PipeResonatorScenario();
//scenario = new HelmholtzResonatorScenario();
scenario = new SodShockTubeScenario();
scenario.Initialize(SampleRate);
var stopwatch = Stopwatch.StartNew();
double lastDrawTime = 0.0;
double drawInterval = 1.0 / DrawFrequency;
double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds;
// Resampling buffer
List<float> simBuffer = new List<float>(4096);
double readIndex = 0.0;
// --- 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);
for (int i = 0; i < 4; i++)
simBuffer.Add(scenario.Process());
// Reset timer after warmup – this is the “real‑time zero”
stopwatch.Restart();
stepCount = 0; // simulation steps start now
long totalSimSteps = simBuffer.Count;
long totalOutputSamples = 0;
// --- Initialise the simulation scenario ---
Simulation.Initialize(SampleRate);
double lastRealTime = stopwatch.Elapsed.TotalSeconds;
const int outputChunk = 256;
float[] outputBuf = new float[outputChunk];
const int chunkSize = 2048;
float[] buffer = new float[chunkSize];
double lastLogTime = 0.0; // for periodic speed printout
while (window.IsOpen)
{
window.DispatchEvents();
double currentRealTime = stopwatch.Elapsed.TotalSeconds;
double dtSpeed = currentRealTime - lastSpeedUpdateTime;
lastSpeedUpdateTime = currentRealTime;
// --- Compute how many audio samples are needed since last frame ---
double currentTime = stopwatch.Elapsed.TotalSeconds;
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)
// Smoothly transition currentSpeed → desiredSpeed
// When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp
double smoothingRate = 8.0; // higher = faster catch‑up
currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed));
// ---------- Generate audio ----------
double targetAudioClock = currentRealTime + 0.05;
while (totalOutputSamples < targetAudioClock * SampleRate && running)
// --- Generate the required number of simulation steps ---
while (samplesNeeded > 0 && running)
{
int toGenerate = (int)Math.Min(
(long)outputChunk,
(long)(targetAudioClock * SampleRate) - totalOutputSamples
);
if (toGenerate <= 0) break;
double maxIndex = readIndex + (toGenerate - 1) * currentSpeed + 2;
int requiredSimIndex = (int)Math.Ceiling(maxIndex);
while (simBuffer.Count - 1 < requiredSimIndex)
{
simBuffer.Add(scenario.Process());
totalSimSteps++;
}
int toGenerate = Math.Min(samplesNeeded, chunkSize);
for (int i = 0; i < toGenerate; i++)
{
int i0 = (int)readIndex;
int i1 = i0 + 1;
double frac = readIndex - i0;
buffer[i] = Simulation.Process();
stepCount++;
}
soundEngine.WriteSamples(buffer, toGenerate);
samplesNeeded -= toGenerate;
}
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);
// --- Display speed ---
double simTime = stepCount / (double)SampleRate;
double wallTime = stopwatch.Elapsed.TotalSeconds;
double speed = (wallTime > 0) ? simTime / wallTime : 0.0;
readIndex += currentSpeed;
// Update window title with instant speed
window.SetTitle($"FluidSim | Speed: {speed:F3}× | Steps: {stepCount}");
while (readIndex >= 1.0 && simBuffer.Count > 2)
// Console log once per second
if (wallTime - lastLogTime >= 1.0)
{
simBuffer.RemoveAt(0);
readIndex -= 1.0;
}
Console.WriteLine($"Speed: {speed:F3}× ({stepCount} steps, {wallTime:F2}s wall)");
lastLogTime = wallTime;
}
int accepted = soundEngine.WriteSamples(outputBuf, toGenerate);
totalOutputSamples += accepted;
if (accepted < toGenerate)
break;
}
// ---------- 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";
window.SetTitle(
$"{toggleHint} Sim: {simTime:F2}s | " +
$"Speed: {currentSpeed:F4}x → {desiredSpeed:F4}x | " +
$"Actual: {actualSpeed:F2}x"
);
// --- Rendering (placeholder) ---
window.Clear(Color.Black);
scenario.Draw(window);
window.Display();
lastDrawTime = currentRealTime;
}
}
// --- Cleanup ---
soundEngine.Dispose();
window.Dispose();
}
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;
}
private static void OnKeyPressed(object? sender, KeyEventArgs e)
{
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;
}
}
}

133
Scenarios/HelmholtzResonatorScenario.cs
View File

@@ -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);
}
}
}

184
Scenarios/PipeResonatorScenario.cs
View File

@@ -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);
}
}
}

23
Scenarios/Scenario.cs
View File

@@ -1,23 +0,0 @@
using SFML.Graphics;
namespace FluidSim.Core
{
public abstract class Scenario
{
/// <summary>
/// Initialize the scenario with a given audio sample rate.
/// </summary>
public abstract void Initialize(int sampleRate);
/// <summary>
/// Advance one simulation step and return an audio sample.
/// The step size is 1 / sampleRate seconds.
/// </summary>
public abstract float Process();
/// <summary>
/// Draw the current simulation state onto the given SFML render target.
/// </summary>
public abstract void Draw(RenderWindow target);
}
}

93
Scenarios/Simulation.cs Normal file
View File

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

158
Scenarios/SodShockTubeScenario.cs
View File

@@ -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);
}
}
}

10
Sources/EffortSource.cs
View File

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

10
Sources/FlowSource.cs
View File

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

12
Utils/Units.cs → Utilities/Units.cs
View File

@@ -1,6 +1,4 @@
using System;
namespace FluidSim.Utils
namespace FluidSim.Utilities
{
public static class Units
{
@@ -21,10 +19,10 @@ namespace FluidSim.Utils
public static double Celsius(double tC) => tC + 273.15;
public static double AreaFromRadius(double radius) =>
Math.PI * (radius) * (radius);
public static double AreaFromRadius(double radius, double unit = mm) =>
Math.PI * (radius * unit) * (radius * unit);
public static double AreaFromDiameter(double diameter) =>
Math.PI * 0.25 * (diameter) * (diameter);
public static double AreaFromDiameter(double diameter, double unit = mm) =>
Math.PI * 0.25 * (diameter * unit) * (diameter * unit);
}
}