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

Compare commits

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

1 Commits
aba9b76530 ... testing

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

26
.vscode/launch.json vendored
View File

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

41
.vscode/tasks.json vendored
View File

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

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

24
Audio/AudioOutputStream.cs → Audio/RingBufferStream.cs
View File

@@ -3,31 +3,22 @@ using SFML.System;
namespace FluidSim.Audio
{
public class AudioOutputStream : SoundStream
internal class RingBufferStream : SoundStream
{
private readonly SimulationRingBuffer _sourceBuffer;
private double _speed = 1.0; // non‑volatile, accessed with Volatile.Read/Write
private readonly RingBuffer ringBuffer;
public AudioOutputStream(SimulationRingBuffer sourceBuffer)
public RingBufferStream(RingBuffer buffer)
{
_sourceBuffer = sourceBuffer;
// 2 channels, 44.1 kHz, stereo
ringBuffer = buffer;
// 2 channels, 44.1 kHz, standard stereo mapping
Initialize(2, 44100, new[] { SoundChannel.FrontLeft, SoundChannel.FrontRight });
}
/// <summary>Playback speed (0.01 … 1.0 or higher for catch‑up).</summary>
public double Speed
{
get => Volatile.Read(ref _speed);
set => Volatile.Write(ref _speed, value);
}
protected override bool OnGetData(out short[] samples)
{
const int monoBlockSize = 512;
const int monoBlockSize = 512; // number of mono samples we'll read
float[] temp = new float[monoBlockSize];
int read = _sourceBuffer.ReadInterpolated(temp, monoBlockSize, Speed);
int read = ringBuffer.Read(temp, monoBlockSize);
samples = new short[monoBlockSize * 2];
if (read > 0)
@@ -40,7 +31,6 @@ namespace FluidSim.Audio
samples[i * 2 + 1] = final; // right
}
}
// Fill rest with silence
for (int i = read * 2; i < samples.Length; i++)
samples[i] = 0;

98
Audio/SimulationRingBuffer.cs
View File

@@ -1,98 +0,0 @@
namespace FluidSim.Audio
{
public class SimulationRingBuffer
{
private readonly float[] _buffer;
private readonly int _capacity;
private int _writeHead; // monotonic, producer only
private int _readHead; // monotonic, consumer advances after consumption
// Consumer interpolation state
private double _readPosFrac;
private bool _consumerInit;
// Events for signalling
private readonly AutoResetEvent _spaceAvailable = new AutoResetEvent(false);
private readonly AutoResetEvent _dataAvailable = new AutoResetEvent(false);
public SimulationRingBuffer(int capacity)
{
if ((capacity & (capacity - 1)) != 0)
throw new ArgumentException("Capacity must be a power of two.");
_capacity = capacity;
_buffer = new float[capacity];
}
// ---------- Producer ----------
public int FreeSpace => _capacity - (_writeHead - Volatile.Read(ref _readHead));
/// <summary>Number of samples currently available for reading (integer count).</summary>
public int AvailableSamples => Volatile.Read(ref _writeHead) - Volatile.Read(ref _readHead);
public void Write(float sample)
{
while (FreeSpace == 0)
_spaceAvailable.WaitOne();
int w = _writeHead;
int mask = _capacity - 1;
_buffer[w & mask] = sample;
Volatile.Write(ref _writeHead, w + 1);
_dataAvailable.Set();
}
public int Write(float[] data, int count)
{
int free = FreeSpace;
int toWrite = Math.Min(count, free);
int w = _writeHead;
int mask = _capacity - 1;
for (int i = 0; i < toWrite; i++)
_buffer[(w + i) & mask] = data[i];
Volatile.Write(ref _writeHead, w + toWrite);
if (toWrite > 0)
_dataAvailable.Set();
return toWrite;
}
// ---------- Consumer ----------
public void ResetConsumer() => _consumerInit = false;
public int ReadInterpolated(float[] dest, int destCount, double speed)
{
if (!_consumerInit)
{
_readPosFrac = Volatile.Read(ref _readHead);
_consumerInit = true;
}
int mask = _capacity - 1;
int writeHead = Volatile.Read(ref _writeHead);
int produced = 0;
for (int i = 0; i < destCount; i++)
{
int idxFloor = (int)_readPosFrac;
int idxCeil = idxFloor + 1;
if (idxCeil >= writeHead)
break;
float y0 = _buffer[idxFloor & mask];
float y1 = _buffer[idxCeil & mask];
double frac = _readPosFrac - idxFloor;
dest[i] = (float)(y0 + (y1 - y0) * frac);
_readPosFrac += speed;
produced++;
}
int newReadHead = (int)_readPosFrac;
if (newReadHead > Volatile.Read(ref _readHead))
{
Volatile.Write(ref _readHead, newReadHead);
_spaceAvailable.Set();
}
return produced;
}
}
}

46
Audio/SoundEngine.cs
View File

@@ -1,45 +1,45 @@
namespace FluidSim.Audio
{
public class SoundEngine : IDisposable
{
private readonly AudioOutputStream _stream;
private bool _isPlaying;
namespace FluidSim.Audio;
public SoundEngine(SimulationRingBuffer sourceBuffer, int bufferCapacity = 16384)
public class SoundEngine : IDisposable
{
private readonly RingBuffer ringBuffer;
private readonly RingBufferStream stream;
private bool isPlaying;
public SoundEngine(int bufferCapacity = 16384)
{
_stream = new AudioOutputStream(sourceBuffer);
ringBuffer = new RingBuffer(bufferCapacity);
stream = new RingBufferStream(ringBuffer);
}
public void Start()
{
if (_isPlaying) return;
_stream.Play();
_isPlaying = true;
if (isPlaying) return;
stream.Play();
isPlaying = true;
}
public void Stop()
{
if (!_isPlaying) return;
_stream.Stop();
_isPlaying = false;
if (!isPlaying) return;
stream.Stop();
isPlaying = false;
float[] drain = new float[ringBuffer.Count];
ringBuffer.Read(drain, drain.Length);
}
public double Speed
{
get => _stream.Speed;
set => _stream.Speed = value;
}
public int WriteSamples(float[] data, int count) =>
ringBuffer.Write(data, count);
public float Volume
{
get => _stream.Volume;
set => _stream.Volume = value;
get => stream.Volume;
set => stream.Volume = value;
}
public void Dispose()
{
Stop();
_stream.Dispose();
}
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;
}
}
}

39
Components/Atmosphere.cs
View File

@@ -1,39 +0,0 @@
using FluidSim.Interfaces;
namespace FluidSim.Components
{
public class Atmosphere : IComponent
{
public float Pressure { get; set; } = 101325f;
public float Temperature { get; set; } = 300f;
public float GasConstant { get; set; } = 287f;
public float Gamma => 1.4f;
public float Density => Pressure / (GasConstant * Temperature);
public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / Density;
public Port Port { get; }
public Atmosphere()
{
Port = new Port { Owner = this };
UpdatePort();
}
public IReadOnlyList<Port> Ports => new[] { Port };
public void UpdateState(float dt)
{
UpdatePort();
}
private void UpdatePort()
{
Port.Pressure = Pressure;
Port.Density = Density;
Port.Temperature = Temperature;
Port.SpecificEnthalpy = SpecificEnthalpy;
Port.AirFraction = 1f;
}
}
}

17
Components/Connection.cs Normal file
View File

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

52
Components/Crankshaft.cs
View File

@@ -1,52 +0,0 @@
using System;
namespace FluidSim.Components
{
public class Crankshaft
{
public float AngularVelocity; // rad/s
public float CrankAngle; // rad, 0 … 4π
public float PreviousAngle;
public float Inertia = 0.2f;
public float FrictionConstant; // N·m
public float FrictionViscous; // N·m per rad/s
private float externalTorque;
public Crankshaft(float initialRPM = 400f)
{
AngularVelocity = initialRPM * 2f * MathF.PI / 60f;
CrankAngle = 0f;
PreviousAngle = 0f;
}
public void AddTorque(float torque) => externalTorque += torque;
public void Step(float dt)
{
if (float.IsNaN(AngularVelocity) || float.IsInfinity(AngularVelocity))
AngularVelocity = 0f;
if (float.IsNaN(externalTorque) || float.IsInfinity(externalTorque))
externalTorque = 0f;
PreviousAngle = CrankAngle;
float friction = FrictionConstant * MathF.Sign(AngularVelocity)
+ FrictionViscous * AngularVelocity;
float netTorque = externalTorque - friction;
float alpha = netTorque / Inertia;
AngularVelocity += alpha * dt;
if (AngularVelocity < 0f) AngularVelocity = 0f;
CrankAngle += AngularVelocity * dt;
if (CrankAngle >= 4f * MathF.PI)
CrankAngle -= 4f * MathF.PI;
else if (CrankAngle < 0f)
CrankAngle += 4f * MathF.PI;
externalTorque = 0f;
}
}
}

281
Components/Cylinder.cs
View File

@@ -1,281 +0,0 @@
using System;
using System.Collections.Generic;
using FluidSim.Interfaces;
namespace FluidSim.Components
{
public class Cylinder : IComponent
{
public Port IntakePort { get; }
public Port ExhaustPort { get; }
public Crankshaft Crankshaft { get; }
private readonly Port[] _ports;
IReadOnlyList<Port> IComponent.Ports => _ports;
public float Bore { get; }
public float Stroke { get; }
public float ConRodLength { get; }
public float CompressionRatio { get; }
public float IVO, IVC, EVO, EVC; // degrees
public float IntakeValveDiameter = 0.03f;
public float ExhaustValveDiameter = 0.028f;
public float IntakeValveLift = 0.005f;
public float ExhaustValveLift = 0.005f;
public float IntakeValveMaxArea => MathF.PI * IntakeValveDiameter * IntakeValveLift;
public float ExhaustValveMaxArea => MathF.PI * ExhaustValveDiameter * ExhaustValveLift;
public float SparkAdvance = 20f;
public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f;
public float StoichiometricAFR = 14.7f;
public float FuelLowerHeatingValue = 44e6f;
public float EnergyVariationFraction = 0.05f;
public float MisfireProbability = 0.0f;
public float CylinderWallArea = 0.02f;
public float HeatTransferCoefficient = 100f;
public float AmbientTemperature = 300f;
public float PhaseOffset; // rad
public float Volume => cylinderVolume;
public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f);
public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f);
public float Mass => _airMass + _exhaustMass;
public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume;
private float cylinderVolume, cylinderEnergy;
private float _airMass, _exhaustMass;
private float trappedAirMass, fuelMass, burnFraction;
private bool combustionActive, fuelInjected;
private float _energyFactor = 1f;
private readonly Random _random = new Random();
private const float Gamma = 1.4f;
private const float GasConstant = 287f;
private const float MaxPressurePa = 200e5f;
private const float MaxTemperatureK = 3500f;
public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio,
float ivo, float ivc, float evo, float evc, Crankshaft crankshaft)
{
Bore = bore; Stroke = stroke; ConRodLength = conRodLength;
CompressionRatio = compressionRatio;
IVO = ivo; IVC = ivc; EVO = evo; EVC = evc;
Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft));
cylinderVolume = clearanceVolume;
float initRho = 1.225f;
_airMass = initRho * clearanceVolume;
_exhaustMass = 0f;
cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f);
IntakePort = new Port { Owner = this };
ExhaustPort = new Port { Owner = this };
_ports = new[] { IntakePort, ExhaustPort };
}
private float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke;
private float clearanceVolume => SweptVolume / (CompressionRatio - 1f);
private float CrankRadius => Stroke * 0.5f;
private float Obliquity => CrankRadius / ConRodLength;
private float CrankDeg =>
((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f;
public float ComputeVolume(float thetaRad)
{
float r = CrankRadius, l = ConRodLength;
float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad);
float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
float x = r * (1f - cosTh) + l * (1f - term);
float area = MathF.PI * 0.25f * Bore * Bore;
return clearanceVolume + area * x;
}
private float ValveLift(float thetaDeg, float opens, float closes, float peakLift)
{
float deg = thetaDeg % 720f;
if (deg < 0f) deg += 720f;
float duration;
float effectiveOpen = opens;
float effectiveClose = closes;
if (closes < opens)
{
// Wrap‑around case (e.g., exhaust: opens near 480°, closes near 30°)
effectiveClose += 720f;
}
duration = effectiveClose - effectiveOpen;
if (duration <= 0f) return 0f;
// Map the angle into the [opens, opens+duration] window
float mapped = deg;
if (mapped < opens) mapped += 720f;
if (mapped < opens || mapped > effectiveClose) return 0f;
float rampDur = duration * 0.25f;
float holdDur = duration - 2f * rampDur;
if (mapped >= opens && mapped < opens + rampDur)
{
float t = (mapped - opens) / rampDur;
return peakLift * t * t * (3f - 2f * t);
}
else if (mapped >= opens + rampDur && mapped < opens + rampDur + holdDur)
{
return peakLift;
}
else if (mapped >= opens + rampDur + holdDur && mapped <= effectiveClose)
{
float t = (mapped - (opens + rampDur + holdDur)) / rampDur;
return peakLift * (1f - t) * (1f - t) * (1f + 2f * t);
}
return 0f;
}
public float IntakeValveArea =>
MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift);
public float ExhaustValveArea =>
MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift);
private float Wiebe(float angleSinceSpark)
{
if (angleSinceSpark < WiebeStart) return 0f;
float phi = (angleSinceSpark - WiebeStart) / WiebeDuration;
if (phi <= 0f) return 0f;
return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f));
}
public void PreStep(float dt)
{
float prevVolume = cylinderVolume;
float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset;
cylinderVolume = ComputeVolume(crankAngleRad);
float dV = cylinderVolume - prevVolume;
float pRel = Pressure - 101325f;
float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad);
float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh);
float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term);
float pistonArea = MathF.PI * 0.25f * Bore * Bore;
Crankshaft.AddTorque(pRel * pistonArea * dxdtheta);
cylinderEnergy -= Pressure * dV;
float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % 720f;
float currDeg = crankAngleRad * 180f / MathF.PI % 720f;
// Intake closing
if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC)
{
trappedAirMass = _airMass;
fuelMass = trappedAirMass / StoichiometricAFR;
fuelInjected = true;
}
// Spark
float sparkAngle = 0f - SparkAdvance;
if (sparkAngle < 0f) sparkAngle += 720f;
bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) ||
(prevDeg > sparkAngle + 360f && currDeg < sparkAngle);
if (crossedSpark && !combustionActive && fuelInjected)
{
if (_random.NextDouble() < MisfireProbability)
{
combustionActive = false;
}
else
{
combustionActive = true; burnFraction = 0f;
float range = EnergyVariationFraction;
_energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f);
}
}
// Combustion
if (combustionActive)
{
float angleSinceSpark = currDeg - sparkAngle;
if (angleSinceSpark < 0f) angleSinceSpark += 720f;
float newFraction = Wiebe(angleSinceSpark);
if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance))
{
newFraction = 1f; combustionActive = false;
float totalMass = _airMass + _exhaustMass;
_airMass = 0f; _exhaustMass = totalMass;
}
fuelInjected = false;
float dFraction = newFraction - burnFraction;
if (dFraction > 0f)
{
float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction;
cylinderEnergy += dQ;
_exhaustMass += fuelMass * dFraction;
burnFraction = newFraction;
}
}
// Heat loss
float dQ_loss = HeatTransferCoefficient * CylinderWallArea *
(Temperature - AmbientTemperature) * dt;
cylinderEnergy -= dQ_loss;
// Update port states
float p = Pressure, rho = Density, T = Temperature;
float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f);
float af = AirFraction;
IntakePort.Pressure = p; IntakePort.Density = rho;
IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af;
ExhaustPort.Pressure = p; ExhaustPort.Density = rho;
ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af;
}
public void UpdateState(float dt)
{
float dmAir = 0f, dmExhaust = 0f, dE = 0f;
foreach (var port in _ports)
{
float mdot = port.MassFlowRate;
float af = mdot >= 0f ? port.AirFraction : AirFraction;
dmAir += mdot * af * dt;
dmExhaust += mdot * (1f - af) * dt;
dE += mdot * port.SpecificEnthalpy * dt;
}
_airMass += dmAir; _exhaustMass += dmExhaust;
cylinderEnergy += dE;
float V = MathF.Max(cylinderVolume, 1e-12f);
float currentP = (Gamma - 1f) * cylinderEnergy / V;
if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f);
float currentRho = (_airMass + _exhaustMass) / V;
float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f);
if (currentT > MaxTemperatureK)
{
float pAtTlimit = currentRho * GasConstant * MaxTemperatureK;
cylinderEnergy = pAtTlimit * V / (Gamma - 1f);
}
float totalMass = _airMass + _exhaustMass;
if (totalMass < 1e-9f)
{
_airMass = 1e-9f; _exhaustMass = 0f;
cylinderEnergy = 101325f * V / (Gamma - 1f);
}
else if (cylinderEnergy < 0f)
{
cylinderEnergy = 101325f * V / (Gamma - 1f);
}
if (_airMass < 0f) _airMass = 0f;
if (_exhaustMass < 0f) _exhaustMass = 0f;
}
}
}

250
Components/Pipe1D.cs Normal file
View File

@@ -0,0 +1,250 @@
using System;
using FluidSim.Interfaces;
namespace FluidSim.Components
{
public class Pipe1D
{
public Port PortA { get; }
public Port PortB { get; }
public double Area => _area;
private int _n;
private double _dx, _dt, _gamma = 1.4, _area;
private double[] _rho, _rhou, _E;
private double _hydraulicDiameter;
private double _rhoLeft, _pLeft, _rhoRight, _pRight;
private bool _leftBCSet, _rightBCSet;
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);
/// <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)
{
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();
}
public void SetUniformState(double rho, double u, double p)
{
double e = p / ((_gamma - 1) * rho);
double Etot = rho * e + 0.5 * rho * u * u;
for (int i = 0; i < _n; i++)
{
_rho[i] = rho;
_rhou[i] = rho * u;
_E[i] = Etot;
}
}
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)
{
_rhoLeft = rhoVol;
_pLeft = pVol;
_leftBCSet = true;
}
public void SetRightVolumeState(double rhoVol, double pVol)
{
_rhoRight = rhoVol;
_pRight = pVol;
_rightBCSet = true;
}
private double GetCellTotalSpecificEnthalpy(int i)
{
double rho = Math.Max(_rho[i], 1e-12);
double u = _rhou[i] / rho;
double p = Pressure(i);
double h = _gamma / (_gamma - 1.0) * p / rho;
return h + 0.5 * u * u;
}
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));
void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
out double fm, out double fp, out double fe)
{
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 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;
if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
else if (Ss >= 0)
{
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 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 * diffSR));
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
}
}
}
}

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

148
Components/Volume0D.cs
View File

@@ -1,133 +1,69 @@
using System;
using System.Collections.Generic;
using FluidSim.Interfaces;
using FluidSim.Interfaces;
namespace FluidSim.Components
{
public class Volume0D : IComponent
public class Volume0D
{
public List<Port> Ports { get; } = new List<Port>();
public Port Port { get; private set; }
private float _airMass;
private float _exhaustMass;
public float InternalEnergy;
public float Volume;
public float Dvdt;
public float Gamma { get; set; } = 1.4f;
public float GasConstant { get; set; } = 287f;
public float AmbientPressure { get; set; } = 101325f;
public double Mass { get; private set; }
public double InternalEnergy { get; private set; }
// ---------- Thermal relaxation to environment ----------
/// <summary>Rate of heat transfer to the surroundings (1/s). 0 = adiabatic.</summary>
public float EnergyRelaxationRate { get; set; } = 10f;
/// <summary>Temperature to relax toward (K). Default is room temperature.</summary>
public float AmbientTemperature { get; set; } = 300f;
public double Gamma { get; set; } = 1.4;
public double GasConstant { get; set; } = 287.0;
public float Mass => _airMass + _exhaustMass;
public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f);
public float Density => Mass / MathF.Max(Volume, 1e-12f);
public float Pressure => (Gamma - 1f) * InternalEnergy / MathF.Max(Volume, 1e-12f);
public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f);
public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / MathF.Max(Density, 1e-12f);
public double Volume { get; set; }
public double dVdt { get; set; }
public Volume0D(float initialVolume, float initialPressure,
float initialTemperature, float gasConstant = 287f, float gamma = 1.4f)
private double _dt;
public double Density => Mass / Volume;
public double Pressure => (Gamma - 1.0) * InternalEnergy / Volume;
public double Temperature => Pressure / (Density * GasConstant);
public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
public Volume0D(double initialVolume, double initialPressure,
double initialTemperature, int sampleRate,
double gasConstant = 287.0, double gamma = 1.4)
{
GasConstant = gasConstant;
Gamma = gamma;
Volume = initialVolume;
Dvdt = 0f;
dVdt = 0.0;
_dt = 1.0 / sampleRate;
float rho0 = initialPressure / (GasConstant * initialTemperature);
_airMass = rho0 * Volume;
_exhaustMass = 0f;
InternalEnergy = (initialPressure * Volume) / (Gamma - 1f);
double rho0 = initialPressure / (GasConstant * initialTemperature);
Mass = rho0 * Volume;
InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
Port = new Port();
PushStateToPort();
}
public Port CreatePort()
public void PushStateToPort()
{
var port = new Port { Owner = this };
port.Pressure = Pressure;
port.Density = Density;
port.Temperature = Temperature;
port.SpecificEnthalpy = SpecificEnthalpy;
port.AirFraction = AirFraction;
Ports.Add(port);
return port;
Port.Pressure = Pressure;
Port.Density = Density;
Port.Temperature = Temperature;
Port.SpecificEnthalpy = SpecificEnthalpy;
}
public void SetPressure(float pressure, float? temperature = null)
public void Integrate()
{
float V = MathF.Max(Volume, 1e-12f);
float T = temperature ?? Temperature;
float rho = pressure / (GasConstant * T);
float totalMass = rho * V;
float af = AirFraction;
_airMass = totalMass * af;
_exhaustMass = totalMass * (1f - af);
InternalEnergy = pressure * V / (Gamma - 1f);
}
double mdot = Port.MassFlowRate;
double h_in = Port.SpecificEnthalpy;
public void UpdateState(float dt)
{
float totalMdotAir = 0f, totalMdotExhaust = 0f, totalEdot = 0f;
foreach (var port in Ports)
{
float mdot = port.MassFlowRate;
float af = mdot >= 0f ? port.AirFraction : AirFraction;
totalMdotAir += mdot * af;
totalMdotExhaust += mdot * (1f - af);
totalEdot += mdot * port.SpecificEnthalpy;
}
double dm = mdot * _dt;
double dE = (mdot * h_in) * _dt - Pressure * dVdt * _dt;
float dAir = totalMdotAir * dt;
float dExhaust = totalMdotExhaust * dt;
float dE = totalEdot * dt - Pressure * Dvdt * dt;
_airMass += dAir;
_exhaustMass += dExhaust;
Mass += dm;
InternalEnergy += dE;
// ---- Thermal relaxation ----
if (EnergyRelaxationRate > 0f)
{
float currentMass = Mass;
if (currentMass > 1e-12f)
{
// Target internal energy: current mass at ambient temperature
float targetE = currentMass * GasConstant * AmbientTemperature / (Gamma - 1f);
float relaxFactor = MathF.Exp(-EnergyRelaxationRate * dt);
InternalEnergy = targetE + (InternalEnergy - targetE) * relaxFactor;
}
}
// Hard physical bounds – prevent NaN and unphysical states
if (Mass < 1e-12) Mass = 1e-12;
if (InternalEnergy < 1e-12) InternalEnergy = 1e-12;
float V = MathF.Max(Volume, 1e-12f);
float totalMass = _airMass + _exhaustMass;
if (totalMass < 1e-9f)
{
_airMass = 1e-9f;
_exhaustMass = 0f;
InternalEnergy = AmbientPressure * V / (Gamma - 1f);
PushStateToPort();
}
else if (InternalEnergy < 0f)
{
InternalEnergy = AmbientPressure * V / (Gamma - 1f);
}
if (_airMass < 0f) _airMass = 0f;
if (_exhaustMass < 0f) _exhaustMass = 0f;
float p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afr = AirFraction;
foreach (var port in Ports)
{
port.Pressure = p;
port.Density = rho;
port.Temperature = T;
port.SpecificEnthalpy = h;
port.AirFraction = afr;
}
}
IReadOnlyList<Port> IComponent.Ports => Ports;
}
}

404
Core/BoundarySystem.cs
View File

@@ -1,404 +0,0 @@
using FluidSim.Components;
using FluidSim.Interfaces;
using System;
namespace FluidSim.Core
{
public class BoundarySystem
{
// ---------- Private constants ----------
private const float Gamma = 1.4f;
private const float Gm1 = Gamma - 1f; // 0.4
private const float Rgas = 287f; // J/(kg·K)
private const float GammaOverGm1 = Gamma / Gm1; // 3.5
public struct OrificeDesc
{
public Port VolumePort;
public int PipeIndex;
public bool IsLeftEnd;
public int AreaIndex;
public float DischargeCoeff;
// --- Inertance support ---
public bool UseInertance;
public float EffectiveLength;
public float CurrentMdot; // kg/s, positive = volume → pipe
// --- Loss coefficient (linear resistance) ---
public float LossCoefficient; // N·s/m⁵ or kg/(m⁴·s)
}
public struct OpenEndDesc
{
public int PipeIndex;
public bool IsLeftEnd;
public float AmbientPressure;
public float Gamma;
public float PipeArea;
public float LastMassFlowRate;
public float LastFacePressure;
}
private OrificeDesc[] _orifices;
private OpenEndDesc[] _openEnds;
private float[] _orificeAreas;
private PipeSystem _pipeSystem;
public BoundarySystem(PipeSystem pipeSystem, int maxOrifices, int maxOpenEnds)
{
_pipeSystem = pipeSystem;
_orifices = new OrificeDesc[maxOrifices];
_openEnds = new OpenEndDesc[maxOpenEnds];
_orificeAreas = new float[maxOrifices];
}
public int OrificeCount { get; private set; }
public int OpenEndCount { get; private set; }
// ---------- Add orifice (no inertance) ----------
public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd,
int areaIndex, float dischargeCoeff = 1f,
float lossCoefficient = 0f)
{
_orifices[OrificeCount] = new OrificeDesc
{
VolumePort = volumePort,
PipeIndex = pipeIndex,
IsLeftEnd = isLeftEnd,
AreaIndex = areaIndex,
DischargeCoeff = dischargeCoeff,
UseInertance = false,
EffectiveLength = 0f,
CurrentMdot = 0f,
LossCoefficient = lossCoefficient
};
OrificeCount++;
}
// ---------- Add orifice with inertance ----------
public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd,
int areaIndex, float dischargeCoeff,
float effectiveLength, float lossCoefficient = 0f)
{
// Reuse the base AddOrifice and then override fields
AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient);
ref var d = ref _orifices[OrificeCount - 1];
d.UseInertance = true;
d.EffectiveLength = effectiveLength;
d.LossCoefficient = lossCoefficient; // store the linear resistance
}
public void AddOpenEnd(int pipeIndex, bool isLeftEnd,
float ambientPressure, float pipeArea, float gamma = 1.4f)
{
int idx = OpenEndCount;
_openEnds[idx] = new OpenEndDesc
{
PipeIndex = pipeIndex,
IsLeftEnd = isLeftEnd,
AmbientPressure = ambientPressure,
Gamma = gamma,
PipeArea = pipeArea
};
OpenEndCount++;
}
public void SetOrificeAreas(float[] areas)
{
for (int i = 0; i < OrificeCount; i++)
_orificeAreas[i] = areas[i];
}
public float GetOpenEndMassFlow(int openEndIndex)
{
if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 0f;
return _openEnds[openEndIndex].LastMassFlowRate;
}
public float GetOpenEndPressure(int openEndIndex)
{
if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 101325f;
return _openEnds[openEndIndex].LastFacePressure;
}
// ---------- Resolve all orifices ----------
public void ResolveOrifices(float dt)
{
for (int i = 0; i < OrificeCount; i++)
{
ref var d = ref _orifices[i];
float area = _orificeAreas[d.AreaIndex];
// Gather volume state
float volP = d.VolumePort?.Pressure ?? 101325f;
float volRho = d.VolumePort?.Density ?? 1.2f;
float volT = d.VolumePort?.Temperature ?? 300f;
float volH = d.VolumePort?.SpecificEnthalpy ?? 0f;
float volAF = d.VolumePort?.AirFraction ?? 1f;
// Gather pipe interior state
var (pipeRho, pipeU, pipeP) = d.IsLeftEnd
? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
: _pipeSystem.GetInteriorStateRight(d.PipeIndex);
float pipeT = pipeP / MathF.Max(pipeRho * Rgas, 1e-12f);
float pipeAF = d.IsLeftEnd
? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
: _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
// ---- Handle closed orifice (area ≈ 0) as a wall ----
if (area < 1e-12f || d.VolumePort == null)
{
var (rInt, uInt, pInt) = d.IsLeftEnd
? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
: _pipeSystem.GetInteriorStateRight(d.PipeIndex);
float afInt = d.IsLeftEnd
? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
: _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
if (d.IsLeftEnd)
_pipeSystem.SetGhostLeft(d.PipeIndex, rInt, -uInt, pInt, afInt);
else
_pipeSystem.SetGhostRight(d.PipeIndex, rInt, -uInt, pInt, afInt);
if (d.VolumePort != null) d.VolumePort.MassFlowRate = 0f;
continue;
}
// ---- Preliminary isentropic solution ----
float mdotEst, rhoFaceEst, uFaceEst, pFaceEst;
if (volP >= pipeP)
{
IsentropicOrifice.Compute(volP, volRho, volT, pipeP, Gamma, Rgas, area, d.DischargeCoeff,
out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
}
else
{
IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, Gamma, Rgas, area, d.DischargeCoeff,
out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst);
mdotEst = -mdotEst;
}
// ---- Compute final mass flow with limiters ----
float mdotFinal, rhoFace, uFace, pFace, airFracGhost;
if (d.UseInertance)
{
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;
float dmdot_dt = (dp - Rlin * d.CurrentMdot) / MathF.Max(inertance, 1e-12f);
float mdotNew = d.CurrentMdot + dmdot_dt * dt;
// Limit outflow from volume (if volume owner is Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
if (mdotNew > maxOut) mdotNew = maxOut;
if (mdotNew < -maxOut) mdotNew = -maxOut;
}
// Limit inflow from pipe – pipe cell cannot be emptied in one step
{
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell); // true cell area, not orifice
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
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);
}
if (float.IsNaN(mdotNew)) mdotNew = 0f;
d.CurrentMdot = mdotNew;
mdotFinal = mdotNew;
rhoFace = mdotFinal >= 0 ? volRho : pipeRho;
pFace = pFaceEst;
uFace = MathF.Abs(mdotFinal) / MathF.Max(rhoFace * area, 1e-12f);
}
else
{
// Standard quasi‑steady orifice
mdotFinal = mdotEst;
rhoFace = rhoFaceEst;
uFace = uFaceEst;
pFace = pFaceEst;
// Limit outflow from volume (if Volume0D)
if (d.VolumePort.Owner is Volume0D vol0)
{
float maxOut = vol0.Mass / dt;
if (mdotFinal > maxOut) mdotFinal = maxOut;
}
// ***** CRITICAL: Limit inflow from pipe – pipe cell cannot be drained *****
if (mdotFinal < 0)
{
int adjCell = d.IsLeftEnd ? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float pipeCellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float maxFromPipe = pipeCellMass / dt;
if (mdotFinal < -maxFromPipe)
mdotFinal = -maxFromPipe;
}
d.CurrentMdot = mdotFinal;
// Limit outflow from cylinder into pipe (positive mdot = volume → pipe)
if (mdotFinal > 0f && d.VolumePort?.Owner is Cylinder cyl)
{
float maxOut = cyl.Mass / dt;
if (mdotFinal > maxOut)
mdotFinal = maxOut;
}
}
// ---- Air fraction for ghost ----
if (mdotFinal >= 0)
airFracGhost = volAF;
else
{
airFracGhost = pipeAF;
if (d.VolumePort != null) d.VolumePort.AirFraction = pipeAF;
}
// ---- Sign convention for velocity ----
if (mdotFinal >= 0 && d.IsLeftEnd) uFace = +uFace;
else if (mdotFinal >= 0 && !d.IsLeftEnd) uFace = -uFace;
else if (mdotFinal < 0 && d.IsLeftEnd) uFace = -uFace;
else if (mdotFinal < 0 && !d.IsLeftEnd) uFace = +uFace;
// ---- Set ghost cells ----
if (d.IsLeftEnd)
_pipeSystem.SetGhostLeft(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
else
_pipeSystem.SetGhostRight(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost);
// ---- Update volume port ----
if (d.VolumePort != null)
{
d.VolumePort.MassFlowRate = -mdotFinal;
if (-mdotFinal >= 0) // mass entering volume (out of pipe)
{
float pipeH = GammaOverGm1 * pipeP / MathF.Max(pipeRho, 1e-12f);
d.VolumePort.SpecificEnthalpy = pipeH;
}
else // mass leaving volume (into pipe)
{
d.VolumePort.SpecificEnthalpy = volH;
}
}
}
}
// ---------- Resolve open ends ----------
public void ResolveOpenEnds(float dt)
{
for (int i = 0; i < OpenEndCount; i++)
{
ref var d = ref _openEnds[i];
var (rhoInt, uInt, pInt) = d.IsLeftEnd
? _pipeSystem.GetInteriorStateLeft(d.PipeIndex)
: _pipeSystem.GetInteriorStateRight(d.PipeIndex);
float afInt = d.IsLeftEnd
? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex)
: _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex);
float gamma = d.Gamma;
float gm1 = gamma - 1f;
float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f));
float pAmb = d.AmbientPressure;
// Characteristic solution (isentropic expansion to ambient)
float Jplus = uInt + 2f * cInt / gm1;
float Jminus = uInt - 2f * cInt / gm1;
float s = pInt / MathF.Pow(rhoInt, gamma);
float rhoIso = MathF.Pow(pAmb / s, 1f / gamma);
float cIso = MathF.Sqrt(gamma * pAmb / MathF.Max(rhoIso, 1e-12f));
float uIso = d.IsLeftEnd
? (Jminus + 2f * cIso / gm1)
: (Jplus - 2f * cIso / gm1);
// Supersonic check
bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt);
if (!supersonic)
{
supersonic = d.IsLeftEnd ? (uIso <= -cIso) : (uIso >= cIso);
}
float rhoGhost, uGhost, pGhost, afGhost;
if (supersonic)
{
rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt;
}
else
{
rhoGhost = rhoIso; uGhost = uIso; pGhost = pAmb;
bool inflow = d.IsLeftEnd ? (uIso >= 0f) : (uIso <= 0f);
afGhost = inflow ? 1f : afInt;
}
// ------- Mass flow limiter -------
int adjCell = d.IsLeftEnd
? _pipeSystem.GetPipeStart(d.PipeIndex)
: _pipeSystem.GetPipeEnd(d.PipeIndex) - 1;
float pipeRhoAdj = _pipeSystem.GetCellDensity(adjCell);
float pipeAreaCell = _pipeSystem.GetCellArea(adjCell);
float pipeDxAdj = _pipeSystem.GetCellDx(adjCell);
float cellMass = pipeRhoAdj * pipeAreaCell * pipeDxAdj;
float area = d.PipeArea;
float mdotRaw = rhoGhost * uGhost * area; // positive out of pipe
if (d.IsLeftEnd) mdotRaw = -mdotRaw; // now positive = out of pipe
// Outflow limit
if (mdotRaw > 0 && mdotRaw * dt > cellMass)
{
mdotRaw = cellMass / dt;
}
// Inflow limit (allow up to 10× cell mass to avoid starving the pipe)
else if (mdotRaw < 0 && -mdotRaw * dt > 10f * cellMass)
{
mdotRaw = -10f * cellMass / dt;
}
// Recompute uGhost from the limited mdot, keeping rhoGhost, pGhost
float mdotMag = MathF.Abs(mdotRaw);
uGhost = mdotMag / MathF.Max(rhoGhost * area, 1e-12f);
if (d.IsLeftEnd)
uGhost = (mdotRaw >= 0f) ? -uGhost : uGhost;
else
uGhost = (mdotRaw >= 0f) ? uGhost : -uGhost;
// Apply ghost
if (d.IsLeftEnd)
_pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
else
_pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost);
d.LastMassFlowRate = mdotRaw;
d.LastFacePressure = pGhost;
}
}
}
}

11
Core/Constants.cs
View File

@@ -1,11 +0,0 @@
namespace FluidSim.Core
{
public static class Constants
{
public const float Gamma = 1.4f;
public const float R_gas = 287f;
public const float P_amb = 101325f;
public const float T_amb = 300f;
public static readonly float Rho_amb = P_amb / (R_gas * T_amb);
}
}

27
Core/GhostBuffer.cs
View File

@@ -1,27 +0,0 @@
namespace FluidSim.Core
{
public class GhostBuffer
{
public float[] Rho, U, P, Y;
public int PipeCount { get; }
public GhostBuffer(int pipeCount)
{
PipeCount = pipeCount;
int size = pipeCount * 2;
Rho = new float[size];
U = new float[size];
P = new float[size];
Y = new float[size];
}
public void Set(int pipeIndex, bool isLeftEnd, float rho, float u, float p, float y)
{
int idx = pipeIndex * 2 + (isLeftEnd ? 0 : 1);
Rho[idx] = rho;
U[idx] = u;
P[idx] = p;
Y[idx] = y;
}
}
}

31
Core/IsentropicOrifice.cs
View File

@@ -1,31 +0,0 @@
using System;
namespace FluidSim.Core
{
public static class IsentropicOrifice
{
public static void Compute(
float pUp, float rhoUp, float TUp,
float pDown, float gamma, float R, float area, float Cd,
out float mdot, out float rhoFace, out float uFace, out float pFace)
{
mdot = 0f; rhoFace = rhoUp; uFace = 0f; pFace = pUp;
if (area <= 0f || pUp <= 0f || rhoUp <= 0f || TUp <= 0f) return;
float pr = MathF.Min(pDown / pUp, 1f);
if (pr < 1e-6f) pr = 1e-6f;
float prCrit = MathF.Pow(2f / (gamma + 1f), gamma / (gamma - 1f));
if (pr < prCrit) pr = prCrit;
float exponent = (gamma - 1f) / gamma;
float M = MathF.Sqrt((2f / (gamma - 1f)) * (MathF.Pow(pr, -exponent) - 1f));
if (float.IsNaN(M)) M = 0f;
float aUp = MathF.Sqrt(gamma * R * TUp);
uFace = M * aUp;
rhoFace = rhoUp * MathF.Pow(pr, 1f / gamma);
pFace = pUp * pr;
mdot = rhoFace * uFace * area * Cd;
}
}
}

99
Core/OrificeBoundary.cs Normal file
View File

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

187
Core/OutdoorExhaustReverb.cs
View File

@@ -1,187 +0,0 @@
using System;
namespace FluidSim.Core
{
public class OutdoorExhaustReverb
{
// ========== Early reflection delays (stereo: left/right) ==========
private readonly DelayLine groundL, groundR;
private readonly DelayLine wall1L, wall1R;
private readonly DelayLine wall2L, wall2R;
// ========== Diffuse tail FDNs (left/right each with 8 channels) ==========
private const int FDN_CHANNELS = 8;
private readonly DelayLine[] fdnL, fdnR;
private readonly float[] stateL, stateR;
private readonly OrthonormalMixer mixerL, mixerR;
private readonly LowPassFilter[] filterL, filterR;
public float DryMix { get; set; } = 1.0f; // direct sound unchanged
public float EarlyMix { get; set; } = 0.12f; // very little early reflection (ground bounce)
public float TailMix { get; set; } = 0.18f; // subtle diffuse tail
public float Feedback { get; set; } = 0.35f; // lower feedback – outdoor doesn't ring
public float DampingFreq { get; set; } = 2500f; // air absorption – high frequencies die quickly
public OutdoorExhaustReverb(int sampleRate)
{
// Early reflections – left/right offset by ~1‑2 ms for stereo width
groundL = new DelayLine((int)(sampleRate * 0.008)); // 8 ms
groundR = new DelayLine((int)(sampleRate * 0.010)); // 10 ms
wall1L = new DelayLine((int)(sampleRate * 0.045));
wall1R = new DelayLine((int)(sampleRate * 0.047));
wall2L = new DelayLine((int)(sampleRate * 0.080));
wall2R = new DelayLine((int)(sampleRate * 0.082));
// FDN delay lengths – prime numbers, offset between L/R
int[] lengthsL = { 3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009 };
int[] lengthsR = { 3217, 4049, 5531, 7043, 8537, 10037, 11519, 13033 };
fdnL = new DelayLine[FDN_CHANNELS];
fdnR = new DelayLine[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
{
int lenL = Math.Min(lengthsL[i], (int)(sampleRate * 0.25));
int lenR = Math.Min(lengthsR[i], (int)(sampleRate * 0.25));
fdnL[i] = new DelayLine(lenL);
fdnR[i] = new DelayLine(lenR);
}
stateL = new float[FDN_CHANNELS];
stateR = new float[FDN_CHANNELS];
mixerL = new OrthonormalMixer(FDN_CHANNELS);
mixerR = new OrthonormalMixer(FDN_CHANNELS);
filterL = new LowPassFilter[FDN_CHANNELS];
filterR = new LowPassFilter[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
{
filterL[i] = new LowPassFilter(sampleRate, DampingFreq);
filterR[i] = new LowPassFilter(sampleRate, DampingFreq);
}
}
/// <summary>Stereo reverb – returns (left, right) sample pair.</summary>
public (float left, float right) ProcessStereo(float drySample)
{
// ---- Early reflections ----
float gL = groundL.ReadWrite(drySample * 0.8f);
float gR = groundR.ReadWrite(drySample * 0.8f);
float w1L = wall1L.ReadWrite(drySample * 0.5f);
float w1R = wall1R.ReadWrite(drySample * 0.5f);
float w2L = wall2L.ReadWrite(drySample * 0.4f);
float w2R = wall2R.ReadWrite(drySample * 0.4f);
float earlyL = (gL + w1L + w2L) * EarlyMix;
float earlyR = (gR + w1R + w2R) * EarlyMix;
// ---- Read diffuse tail ----
float[] delOutL = new float[FDN_CHANNELS];
float[] delOutR = new float[FDN_CHANNELS];
for (int i = 0; i < FDN_CHANNELS; i++)
{
delOutL[i] = fdnL[i].Read();
delOutR[i] = fdnR[i].Read();
}
// Mix via orthonormal matrix
float[] mixL = new float[FDN_CHANNELS];
float[] mixR = new float[FDN_CHANNELS];
mixerL.Process(delOutL, mixL);
mixerR.Process(delOutR, mixR);
// Feedback + air absorption
for (int i = 0; i < FDN_CHANNELS; i++)
{
stateL[i] = drySample * 0.15f + Feedback * mixL[i];
stateL[i] = filterL[i].Process(stateL[i]);
fdnL[i].Write(stateL[i]);
stateR[i] = drySample * 0.15f + Feedback * mixR[i];
stateR[i] = filterR[i].Process(stateR[i]);
fdnR[i].Write(stateR[i]);
}
float tailL = 0.0f, tailR = 0.0f;
for (int i = 0; i < FDN_CHANNELS; i++)
{
tailL += delOutL[i];
tailR += delOutR[i];
}
tailL *= TailMix;
tailR *= TailMix;
float left = drySample * DryMix + earlyL + tailL;
float right = drySample * DryMix + earlyR + tailR;
return (left, right);
}
/// <summary>Mono fallback – sums left+right / 2.</summary>
public float Process(float drySample)
{
var (l, r) = ProcessStereo(drySample);
return MathF.Tanh((l + r) * 0.5f);
}
// ========== Helper classes ==========
private class DelayLine
{
private float[] buffer;
private int writePos;
public DelayLine(int length)
{
buffer = new float[Math.Max(length, 1)];
}
public float Read() => buffer[writePos];
public void Write(float value)
{
buffer[writePos] = value;
writePos = (writePos + 1) % buffer.Length;
}
public float ReadWrite(float value)
{
float outVal = buffer[writePos];
buffer[writePos] = value;
writePos = (writePos + 1) % buffer.Length;
return outVal;
}
}
private class LowPassFilter
{
private float b0, a1, y1;
private float sampleRate;
public LowPassFilter(int sampleRate, float cutoff)
{
this.sampleRate = sampleRate;
SetCutoff(cutoff);
}
public void SetCutoff(float cutoff)
{
float w = 2 * (float)Math.PI * cutoff / sampleRate;
float a0 = 1 + w;
b0 = w / a0;
a1 = (1 - w) / a0;
}
public float Process(float x)
{
float y = b0 * x - a1 * y1;
y1 = y;
return y;
}
}
private class OrthonormalMixer
{
private int size;
public OrthonormalMixer(int size) => this.size = size;
public void Process(float[] input, float[] output)
{
float sum = 0;
for (int i = 0; i < size; i++) sum += input[i];
float factor = 2.0f / size;
for (int i = 0; i < size; i++)
output[i] = factor * sum - input[i];
}
}
}
}

693
Core/Pipesystem.cs
View File

@@ -1,693 +0,0 @@
using System;
using System.Diagnostics;
using System.Numerics;
namespace FluidSim.Core
{
public class PipeSystem
{
// ---------- Master arrays ----------
private float[] _rho, _rhou, _E, _Y;
private readonly float[] _area;
private readonly float[] _dx;
private readonly int[] _pipeStart;
private readonly int[] _pipeEnd;
private readonly int _totalCells; // original cell count (visible)
private readonly int _allCells; // total allocated (padded to Vector<float>.Count)
private readonly int _pipeCount;
// Derived state – _p is kept for visualization
private float[] _p;
// Flux arrays for faces INTERNAL to a single pipe (size = _allCells + 1)
// Only valid for faces that are NOT pipe boundaries.
private float[] _fluxM, _fluxP, _fluxE, _fluxY;
// Per‑pipe boundary flux buffers (size = _pipeCount)
private float[] _leftFluxM, _leftFluxP, _leftFluxE, _leftFluxY;
private float[] _rightFluxM, _rightFluxP, _rightFluxE, _rightFluxY;
// Damping and relaxation
private float[] _dampingFactors;
private float[] _relaxFactors;
private bool _applyDamping;
private bool _applyRelax;
// Ghost buffer (per‑pipe ghost states)
private readonly GhostBuffer _ghost;
// Precomputed flag: true if a face is a pipe boundary (start or end)
private readonly bool[] _isPipeBoundaryFace;
// ---------- Physical constants ----------
private const float Gamma = 1.4f;
private const float Gm1 = 0.4f;
private const float Gm1Inv = 1f / Gm1; // 2.5
private const float GammaOverGm1 = Gamma / Gm1; // 3.5
private float _coeffBase;
private float _relaxRate;
private float _ambientPressure = 101325f;
private float _ambientEnergyRef;
public float DampingMultiplier
{
set
{
_coeffBase = 0.1f * value;
_applyDamping = _coeffBase != 0f;
}
}
public float EnergyRelaxationRate
{
set
{
_relaxRate = value;
_applyRelax = _relaxRate != 0f;
}
}
public float AmbientPressure
{
set
{
_ambientPressure = value;
_ambientEnergyRef = value * Gm1Inv;
}
}
// ---------- Profiling ----------
public bool EnableProfiling { get; set; }
private long _profFluxTicks;
private long _profUpdateTicks;
private long _profCallCount;
// ---------- Construction ----------
public PipeSystem(int totalCells, int[] pipeStart, int[] pipeEnd,
float[] area, float[] dx,
float initialRho, float initialU, float initialP)
{
_pipeStart = pipeStart;
_pipeEnd = pipeEnd;
_pipeCount = pipeStart.Length;
_totalCells = totalCells;
_area = area;
_dx = dx;
// Pad to SIMD width so all vectorized loops cover the whole data
int vecSize = Vector<float>.Count;
_allCells = totalCells % vecSize == 0 ? totalCells : totalCells + vecSize - (totalCells % vecSize);
_rho = new float[_allCells];
_rhou = new float[_allCells];
_E = new float[_allCells];
_Y = new float[_allCells];
_p = new float[_allCells]; // pressure for drawing
int faceCount = _allCells + 1;
_fluxM = new float[faceCount];
_fluxP = new float[faceCount];
_fluxE = new float[faceCount];
_fluxY = new float[faceCount];
// Per‑pipe boundary flux buffers
_leftFluxM = new float[_pipeCount];
_leftFluxP = new float[_pipeCount];
_leftFluxE = new float[_pipeCount];
_leftFluxY = new float[_pipeCount];
_rightFluxM = new float[_pipeCount];
_rightFluxP = new float[_pipeCount];
_rightFluxE = new float[_pipeCount];
_rightFluxY = new float[_pipeCount];
_dampingFactors = new float[_allCells];
_relaxFactors = new float[_allCells];
_applyDamping = _coeffBase != 0f;
_applyRelax = _relaxRate != 0f;
_ghost = new GhostBuffer(_pipeCount);
_ambientEnergyRef = initialP * Gm1Inv;
// Mark faces that coincide with a pipe boundary (start or end)
_isPipeBoundaryFace = new bool[faceCount];
for (int p = 0; p < _pipeCount; p++)
{
_isPipeBoundaryFace[_pipeStart[p]] = true;
_isPipeBoundaryFace[_pipeEnd[p]] = true;
}
// Initialize uniform state
float initE = initialP / (Gm1 * initialRho);
float rhoE = initialRho * initE + 0.5f * initialRho * initialU * initialU;
for (int i = 0; i < totalCells; i++)
{
_rho[i] = initialRho;
_rhou[i] = initialRho * initialU;
_E[i] = rhoE;
_Y[i] = 1f;
}
}
// ---------- Ghost setters (for BoundarySystem) ----------
public void SetGhostLeft(int pipeIndex, float rho, float u, float p, float y)
=> _ghost.Set(pipeIndex, true, rho, u, p, y);
public void SetGhostRight(int pipeIndex, float rho, float u, float p, float y)
=> _ghost.Set(pipeIndex, false, rho, u, p, y);
// ---------- Public read methods ----------
public int TotalCells => _totalCells;
public int PipeCount => _pipeCount;
public int GetPipeStart(int pipeIdx) => _pipeStart[pipeIdx];
public int GetPipeEnd(int pipeIdx) => _pipeEnd[pipeIdx];
public float GetCellPressure(int i) => _p[i];
public float GetCellDensity(int i) => _rho[i];
public float GetCellDx(int i) => _dx[i];
public float GetCellArea(int i) => _area[i];
public float GetCellVelocity(int i)
{
float rho = _rho[i];
return rho > 1e-12f ? _rhou[i] / rho : 0f;
}
public float GetCellAirFraction(int i) => _Y[i];
public (float rho, float u, float p) GetInteriorStateLeft(int pipeIdx)
{
int i = _pipeStart[pipeIdx];
float rho = _rho[i];
float rhou = _rhou[i];
float u = rhou / MathF.Max(rho, 1e-12f);
float p = Gm1 * (_E[i] - 0.5f * rhou * u);
return (rho, u, p);
}
public (float rho, float u, float p) GetInteriorStateRight(int pipeIdx)
{
int i = _pipeEnd[pipeIdx] - 1;
float rho = _rho[i];
float rhou = _rhou[i];
float u = rhou / MathF.Max(rho, 1e-12f);
float p = Gm1 * (_E[i] - 0.5f * rhou * u);
return (rho, u, p);
}
public float GetInteriorAirFractionLeft(int pipeIdx) => _Y[_pipeStart[pipeIdx]];
public float GetInteriorAirFractionRight(int pipeIdx) => _Y[_pipeEnd[pipeIdx] - 1];
public void SetCellState(int i, float rho, float u, float p, float y = 1f)
{
if (i < 0 || i >= _totalCells) return;
_rho[i] = rho;
_rhou[i] = rho * u;
_E[i] = p * Gm1Inv + 0.5f * rho * u * u;
_Y[i] = y;
}
// ---------- Main step ----------
public void SimulateStep(float dt)
{
long t0 = 0, t1 = 0;
if (EnableProfiling)
{
_profCallCount++;
t0 = Stopwatch.GetTimestamp();
}
ComputeFluxes(dt);
if (EnableProfiling)
{
t1 = Stopwatch.GetTimestamp();
_profFluxTicks += (t1 - t0);
t0 = t1;
}
UpdateCells(dt);
if (EnableProfiling)
{
t1 = Stopwatch.GetTimestamp();
_profUpdateTicks += (t1 - t0);
}
}
// ---------- Flux computation ----------
private void ComputeFluxes(float dt)
{
float fm, fp, fe;
int vecSize = Vector<float>.Count;
// ---- 1. Left ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeStart[p];
int ghostIdx = p * 2;
float rL = _ghost.Rho[ghostIdx];
float uL = _ghost.U[ghostIdx];
float pL = _ghost.P[ghostIdx];
float YL = _ghost.Y[ghostIdx];
float cL = MathF.Sqrt(Gamma * pL / MathF.Max(rL, 1e-12f));
float rR = _rho[idx], rhouR = _rhou[idx];
float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
float uR = rhouR * invRhoR;
float pR = Gm1 * (_E[idx] - 0.5f * rhouR * uR);
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
float YR = _Y[idx];
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_leftFluxM[p] = fm; _leftFluxP[p] = fp; _leftFluxE[p] = fe;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_leftFluxY[p] = fy;
}
// ---- 2. Right ghost boundaries → per‑pipe buffers ----
for (int p = 0; p < _pipeCount; p++)
{
int idx = _pipeEnd[p] - 1;
int ghostIdx = p * 2 + 1;
float rR = _ghost.Rho[ghostIdx];
float uR = _ghost.U[ghostIdx];
float pR = _ghost.P[ghostIdx];
float YR = _ghost.Y[ghostIdx];
float cR = MathF.Sqrt(Gamma * pR / MathF.Max(rR, 1e-12f));
float rL = _rho[idx], rhouL = _rhou[idx];
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
float uL = rhouL * invRhoL;
float pL = Gm1 * (_E[idx] - 0.5f * rhouL * uL);
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
float YL = _Y[idx];
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_rightFluxM[p] = fm; _rightFluxP[p] = fp; _rightFluxE[p] = fe;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_rightFluxY[p] = fy;
}
// ---- 3. Interior faces (skip pipe boundaries) → global flux arrays ----
for (int face = 1; face < _totalCells; face++)
{
// Skip faces that belong to a pipe boundary (they are already handled)
if (_isPipeBoundaryFace[face])
continue;
// Try to vectorize a block of contiguous non‑boundary faces
if (face + vecSize - 1 < _totalCells)
{
bool canVectorize = true;
for (int f = face; f < face + vecSize; f++)
{
if (_isPipeBoundaryFace[f])
{
canVectorize = false;
break;
}
}
if (canVectorize)
{
// --- Vectorised block ---
var rhoL = new Vector<float>(_rho, face - 1);
var rhouL = new Vector<float>(_rhou, face - 1);
var EL = new Vector<float>(_E, face - 1);
var YL = new Vector<float>(_Y, face - 1);
var rhoR = new Vector<float>(_rho, face);
var rhouR = new Vector<float>(_rhou, face);
var ER = new Vector<float>(_E, face);
var YR = new Vector<float>(_Y, face);
var invRhoL = Vector<float>.One / Vector.Max(rhoL, new Vector<float>(1e-12f));
var invRhoR = Vector<float>.One / Vector.Max(rhoR, new Vector<float>(1e-12f));
var uL = rhouL * invRhoL;
var uR = rhouR * invRhoR;
var kinL = 0.5f * rhouL * uL;
var kinR = 0.5f * rhouR * uR;
var pL = Gm1 * (EL - kinL);
var pR = Gm1 * (ER - kinR);
var cL = Vector.SquareRoot(Gamma * pL * invRhoL);
var cR = Vector.SquareRoot(Gamma * pR * invRhoR);
var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL;
var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR;
var FmL = rhoL * uL;
var FpL = rhoL * uL * uL + pL;
var FeL = (rhoL * ELs + pL) * uL;
var FmR = rhoR * uR;
var FpR = rhoR * uR * uR + pR;
var FeR = (rhoR * ERs + pR) * uR;
var absUL = Vector.Abs(uL);
var absUR = Vector.Abs(uR);
var alpha = Vector.Max(absUL + cL, absUR + cR);
var fmVec = 0.5f * (FmL + FmR) - 0.5f * alpha * (rhoR - rhoL);
var fpVec = 0.5f * (FpL + FpR) - 0.5f * alpha * (rhouR - rhouL);
var feVec = 0.5f * (FeL + FeR) - 0.5f * alpha * (rhoR * ERs - rhoL * ELs);
var fyL = FmL * YL;
var fyR = FmR * YR;
var fyVec = 0.5f * (fyL + fyR) - 0.5f * alpha * (rhoR * YR - rhoL * YL);
fmVec.CopyTo(_fluxM, face);
fpVec.CopyTo(_fluxP, face);
feVec.CopyTo(_fluxE, face);
fyVec.CopyTo(_fluxY, face);
face += vecSize - 1; // loop increment will add 1
continue;
}
}
// --- Scalar fallback for a single interior face ---
{
int iL = face - 1, iR = face;
float rL = _rho[iL], rhouL = _rhou[iL];
float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f));
float uL = rhouL * invRhoL;
float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL);
float cL = MathF.Sqrt(Gamma * pL * invRhoL);
float YL = _Y[iL];
float rR = _rho[iR], rhouR = _rhou[iR];
float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f));
float uR = rhouR * invRhoR;
float pR = Gm1 * (_E[iR] - 0.5f * rhouR * uR);
float cR = MathF.Sqrt(Gamma * pR * invRhoR);
float YR = _Y[iR];
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe);
_fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy);
_fluxY[face] = fy;
}
}
}
// ---------- Cell update (per pipe, using correct boundary fluxes) ----------
private void UpdateCells(float dt)
{
int vecSize = Vector<float>.Count;
float dtRelax = -_relaxRate * dt;
// Precompute damping and relaxation factors globally
if (_applyDamping)
{
for (int i = 0; i < _totalCells; i++)
{
float rho = _rho[i];
_dampingFactors[i] = rho > 1e-12f
? MathF.Exp(-_coeffBase * dt / rho)
: 1f;
}
}
if (_applyRelax)
{
float relaxVal = MathF.Exp(dtRelax);
for (int i = 0; i < _totalCells; i++)
_relaxFactors[i] = relaxVal;
}
// Update each pipe separately
for (int p = 0; p < _pipeCount; p++)
{
int start = _pipeStart[p];
int end = _pipeEnd[p]; // exclusive
int len = end - start;
if (len == 0) continue;
// ------- Left boundary cell (i = start) ------
{
int i = start;
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
// left face: always the pipe's left boundary flux
float fluxM_L = _leftFluxM[p];
float fluxP_L = _leftFluxP[p];
float fluxE_L = _leftFluxE[p];
float fluxY_L = _leftFluxY[p];
// right face: depends on pipe length
float fluxM_R, fluxP_R, fluxE_R, fluxY_R;
if (len == 1)
{
// Only one cell: right face is the pipe's right boundary flux
fluxM_R = _rightFluxM[p];
fluxP_R = _rightFluxP[p];
fluxE_R = _rightFluxE[p];
fluxY_R = _rightFluxY[p];
}
else
{
// interior face (global flux at index i+1)
fluxM_R = _fluxM[i + 1];
fluxP_R = _fluxP[i + 1];
fluxE_R = _fluxE[i + 1];
fluxY_R = _fluxY[i + 1];
}
float dtdx = dt / _dx[i];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[i];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[i] = rhoNew;
_rhou[i] = rhouNew;
_E[i] = ENew;
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
// ------- Interior cells (i = start+1 to end-2) ------
if (len > 2)
{
int iCell = start + 1;
int iEnd = end - 1; // exclusive upper bound
// Vectorised path for interior cells (if available)
for (; iCell <= iEnd - vecSize; iCell += vecSize)
{
var rhoOld = new Vector<float>(_rho, iCell);
var rhouOld = new Vector<float>(_rhou, iCell);
var EOld = new Vector<float>(_E, iCell);
var YOld = new Vector<float>(_Y, iCell);
var fluxM_L = new Vector<float>(_fluxM, iCell);
var fluxP_L = new Vector<float>(_fluxP, iCell);
var fluxE_L = new Vector<float>(_fluxE, iCell);
var fluxY_L = new Vector<float>(_fluxY, iCell);
var fluxM_R = new Vector<float>(_fluxM, iCell + 1);
var fluxP_R = new Vector<float>(_fluxP, iCell + 1);
var fluxE_R = new Vector<float>(_fluxE, iCell + 1);
var fluxY_R = new Vector<float>(_fluxY, iCell + 1);
var dtdx = new Vector<float>(dt) / new Vector<float>(_dx, iCell);
var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
var ENew = EOld - dtdx * (fluxE_R - fluxE_L);
var rhoYOld = rhoOld * YOld;
var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping)
rhouNew *= new Vector<float>(_dampingFactors, iCell);
if (_applyRelax)
{
var ambRef = new Vector<float>(_ambientEnergyRef);
var relax = new Vector<float>(_relaxFactors, iCell);
ENew = ambRef + (ENew - ambRef) * relax;
}
rhoNew = Vector.Max(rhoNew, new Vector<float>(1e-12f));
var kinNew = 0.5f * rhouNew * rhouNew / rhoNew;
var eMin = new Vector<float>(100f * Gm1Inv) + kinNew;
ENew = Vector.Max(ENew, eMin);
rhoNew.CopyTo(_rho, iCell);
rhouNew.CopyTo(_rhou, iCell);
ENew.CopyTo(_E, iCell);
var yNew = rhoYNew / rhoNew;
yNew = Vector.Min(Vector.Max(yNew, Vector<float>.Zero), Vector<float>.One);
yNew.CopyTo(_Y, iCell);
}
// Scalar remainder for interior cells
for (; iCell < iEnd; iCell++)
{
float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell];
float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell];
float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1];
float dtdx = dt / _dx[iCell];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[iCell];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[iCell] = rhoNew;
_rhou[iCell] = rhouNew;
_E[iCell] = ENew;
_Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
}
// ------- Right boundary cell (i = end-1, if len > 1) ------
if (len > 1)
{
int i = end - 1;
float rhoOld = _rho[i], rhouOld = _rhou[i], EOld = _E[i], YOld = _Y[i];
// left face
float fluxM_L, fluxP_L, fluxE_L, fluxY_L;
if (len == 2)
{
// Only two cells: left face is the pipe's left boundary flux
fluxM_L = _leftFluxM[p];
fluxP_L = _leftFluxP[p];
fluxE_L = _leftFluxE[p];
fluxY_L = _leftFluxY[p];
}
else
{
// interior face (global flux at i)
fluxM_L = _fluxM[i];
fluxP_L = _fluxP[i];
fluxE_L = _fluxE[i];
fluxY_L = _fluxY[i];
}
// right face: always the pipe's right boundary flux
float fluxM_R = _rightFluxM[p];
float fluxP_R = _rightFluxP[p];
float fluxE_R = _rightFluxE[p];
float fluxY_R = _rightFluxY[p];
float dtdx = dt / _dx[i];
float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L);
float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L);
float ENew = EOld - dtdx * (fluxE_R - fluxE_L);
float rhoYOld = rhoOld * YOld;
float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L);
if (_applyDamping) rhouNew *= _dampingFactors[i];
if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[i];
rhoNew = MathF.Max(rhoNew, 1e-12f);
float kin = 0.5f * rhouNew * rhouNew / rhoNew;
float eMin = 100f * Gm1Inv + kin;
ENew = MathF.Max(ENew, eMin);
_rho[i] = rhoNew;
_rhou[i] = rhouNew;
_E[i] = ENew;
_Y[i] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f);
}
}
// Recompute pressure for all cells (for visualization)
for (int i = 0; i < _totalCells; i++)
{
float rho = _rho[i];
float rhou = _rhou[i];
float u = rhou / MathF.Max(rho, 1e-12f);
_p[i] = Gm1 * (_E[i] - 0.5f * rhou * u);
}
}
// ---------- Scalar flux helpers ----------
private static void LaxFlux(float rL, float uL, float pL, float cL,
float rR, float uR, float pR, float cR,
out float fm, out float fp, out float fe)
{
float EL = pL * Gm1Inv / rL + 0.5f * uL * uL;
float ER = pR * Gm1Inv / rR + 0.5f * uR * uR;
float FmL = rL * uL;
float FpL = rL * uL * uL + pL;
float FeL = (rL * EL + pL) * uL;
float FmR = rR * uR;
float FpR = rR * uR * uR + pR;
float FeR = (rR * ER + pR) * uR;
float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR);
fm = 0.5f * (FmL + FmR) - 0.5f * alpha * (rR - rL);
fp = 0.5f * (FpL + FpR) - 0.5f * alpha * (rR * uR - rL * uL);
fe = 0.5f * (FeL + FeR) - 0.5f * alpha * (rR * ER - rL * EL);
}
private static void ScalarFlux(float rL, float uL, float YL,
float rR, float uR, float YR,
float alpha, out float fy)
{
float FyL = rL * uL * YL;
float FyR = rR * uR * YR;
fy = 0.5f * (FyL + FyR) - 0.5f * alpha * (rR * YR - rL * YL);
}
public int GetRequiredSubSteps(float dtGlobal, float cflTarget = 0.8f)
{
float maxW = 0f;
for (int i = 0; i < _totalCells; i++)
{
float rho = MathF.Max(_rho[i], 1e-12f);
float u = MathF.Abs(_rhou[i] / rho);
float p = Gm1 * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
float c = MathF.Sqrt(Gamma * p / rho);
float w = u + c;
if (w > maxW) maxW = w;
}
maxW = MathF.Max(maxW, 1e-8f);
float minDx = _dx.Min(); // need using System.Linq;
return Math.Max(1, (int)MathF.Ceiling(dtGlobal * maxW / (cflTarget * minDx)));
}
// ---------- Profiling report ----------
public string GetProfileReport()
{
if (!EnableProfiling || _profCallCount == 0)
return "Pipe profiling disabled or no data.";
double freq = Stopwatch.Frequency;
long totalTicks = _profFluxTicks + _profUpdateTicks;
if (totalTicks == 0) return "No pipe profile data collected.";
double totalMs = totalTicks * 1000.0 / freq;
double avgCallUs = totalMs * 1000.0 / _profCallCount;
double fluxMs = _profFluxTicks * 1000.0 / freq;
double updateMs = _profUpdateTicks * 1000.0 / freq;
double fluxAvgUs = fluxMs * 1000.0 / _profCallCount;
double updateAvgUs = updateMs * 1000.0 / _profCallCount;
string report = $" Pipe kernel (over {_profCallCount} calls, total {totalMs:F2} ms, avg {avgCallUs:F2} µs/call):\n";
report += $" Fluxes (incl. primitives): {fluxMs:F2} ms ({_profFluxTicks * 100.0 / totalTicks:F1}%), avg {fluxAvgUs:F2} µs/call\n";
report += $" Update cells: {updateMs:F2} ms ({_profUpdateTicks * 100.0 / totalTicks:F1}%), avg {updateAvgUs:F2} µs/call\n";
_profFluxTicks = 0;
_profUpdateTicks = 0;
_profCallCount = 0;
return report;
}
}
}

175
Core/Solver.cs
View File

@@ -1,7 +1,4 @@
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using FluidSim.Audio;
using FluidSim.Components;
using FluidSim.Interfaces;
@@ -9,94 +6,120 @@ namespace FluidSim.Core
{
public class Solver
{
private readonly List<IComponent> _components = new();
private PipeSystem _pipeSystem;
private BoundarySystem _boundarySystem;
private double _dt;
private readonly List<Volume0D> _volumes = new();
private readonly List<Pipe1D> _pipes = new();
private readonly List<Connection> _connections = new();
public int SubStepCount { get; set; } = 4;
public bool EnableProfiling { get; set; } = false;
public float LastSample { get; private set; }
private long _stepCount;
private long _ticksOrifice, _ticksOpenEnd, _ticksPipe, _ticksUpdate;
public void SetTimeStep(double dt) => _dt = dt;
public void AddComponent(IComponent component) => _components.Add(component);
public void SetPipeSystem(PipeSystem pipeSystem)
{
_pipeSystem = pipeSystem;
}
public void SetBoundarySystem(BoundarySystem boundarySystem)
{
_boundarySystem = boundarySystem;
}
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()
{
if (_pipeSystem == null || _boundarySystem == null) return;
// 1. Publish volume states to their own ports
foreach (var v in _volumes)
v.PushStateToPort();
int nSub = _pipeSystem.GetRequiredSubSteps((float)_dt, 0.8f);
nSub = Math.Max(nSub, SubStepCount); // never go below fixed minimum
float dtSub = (float)(_dt / nSub);
for (int sub = 0; sub < nSub; sub++)
// 2. Handle direct volume‑to‑volume connections
foreach (var conn in _connections)
{
long t0;
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);
t0 = Stopwatch.GetTimestamp();
_boundarySystem.ResolveOrifices(dtSub);
_ticksOrifice += Stopwatch.GetTimestamp() - t0;
if (volA == null || volB == null) continue;
t0 = Stopwatch.GetTimestamp();
_boundarySystem.ResolveOpenEnds(dtSub);
_ticksOpenEnd += Stopwatch.GetTimestamp() - t0;
double pA = volA.Pressure, rhoA = volA.Density;
double pB = volB.Pressure, rhoB = volB.Density;
t0 = Stopwatch.GetTimestamp();
_pipeSystem.SimulateStep(dtSub);
_ticksPipe += Stopwatch.GetTimestamp() - t0;
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
}
}
}
long tUS = Stopwatch.GetTimestamp();
foreach (var comp in _components)
comp.UpdateState((float)_dt);
_ticksUpdate += Stopwatch.GetTimestamp() - tUS;
_stepCount++;
if (_stepCount % 5000 == 0 && EnableProfiling)
// 3. Pipe‑volume boundary conditions – unchanged
foreach (var conn in _connections)
{
double freq = Stopwatch.Frequency;
double total = _ticksOrifice + _ticksOpenEnd + _ticksPipe + _ticksUpdate;
double avgStepUs = (total / freq) * 1e6 / 5000.0;
int orificeCalls = 5000 * nSub;
int updateCalls = 5000;
double orificeMs = _ticksOrifice * 1000.0 / freq;
double openEndMs = _ticksOpenEnd * 1000.0 / freq;
double pipeMs = _ticksPipe * 1000.0 / freq;
double updateMs = _ticksUpdate * 1000.0 / freq;
double orificeAvgUs = orificeMs * 1000.0 / orificeCalls;
double openEndAvgUs = openEndMs * 1000.0 / orificeCalls;
double pipeAvgUs = pipeMs * 1000.0 / orificeCalls;
double updateAvgUs = updateMs * 1000.0 / updateCalls;
Console.WriteLine($"--- Solver ({5000} steps, nSub={nSub}) ---");
Console.WriteLine($" Average step: {avgStepUs:F2} µs");
Console.WriteLine($" Orifice: {orificeMs:F2} ms ({(double)_ticksOrifice / total * 100:F1}%), avg {orificeAvgUs:F2} µs/call");
Console.WriteLine($" OpenEnd: {openEndMs:F2} ms ({(double)_ticksOpenEnd / total * 100:F1}%), avg {openEndAvgUs:F2} µs/call");
Console.WriteLine($" Pipe: {pipeMs:F2} ms ({(double)_ticksPipe / total * 100:F1}%), avg {pipeAvgUs:F2} µs/call");
Console.WriteLine($" Update: {updateMs:F2} ms ({(double)_ticksUpdate / total * 100:F1}%), avg {updateAvgUs:F2} µs/call");
// Pipe internal breakdown (with per-phase averages)
if (_pipeSystem.EnableProfiling)
{
Console.WriteLine(_pipeSystem.GetProfileReport());
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
SetVolumeBC(conn.PortA, conn.PortB);
else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
SetVolumeBC(conn.PortB, conn.PortA);
}
_ticksOrifice = _ticksOpenEnd = _ticksPipe = _ticksUpdate = 0;
// 4. Run pipe simulations
foreach (var p in _pipes)
p.Simulate();
// 5. Transfer pipe‑to‑volume flows – unchanged
foreach (var conn in _connections)
{
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
TransferPipeToVolume(conn.PortA, conn.PortB);
else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
TransferPipeToVolume(conn.PortB, conn.PortA);
}
// 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);
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)
{
Pipe1D pipe = GetPipe(pipePort);
if (pipe == null) return;
bool isLeft = pipe.PortA == pipePort;
if (isLeft)
pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
else
pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
}
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 → pipe, volume’s own static enthalpy is used (already set)
}
}
}

34
Core/SoundProcessor.cs
View File

@@ -1,34 +0,0 @@
using System;
namespace FluidSim.Core
{
public class SoundProcessor
{
private readonly float dt;
private readonly float scaleFactor; // 1 / (4π r)
private float flowLP, prevMassFlowOut, smoothDMdt;
private readonly float lpAlpha, alpha;
public float Gain = 1f;
public SoundProcessor(int sampleRate, float listenerDistance = 1f)
{
dt = 1f / sampleRate;
scaleFactor = 1f / (4f * MathF.PI * listenerDistance);
float tau = 0.02f;
alpha = MathF.Exp(-dt / tau);
float tauLP = 0.005f;
lpAlpha = MathF.Exp(-dt / tauLP);
}
public float Process(float massFlowOut)
{
flowLP = lpAlpha * flowLP + (1f - lpAlpha) * massFlowOut;
float rawDerivative = (flowLP - prevMassFlowOut) / dt;
prevMassFlowOut = flowLP;
smoothDMdt = alpha * smoothDMdt + (1f - alpha) * rawDerivative;
float pressure = smoothDMdt * scaleFactor * Gain;
return MathF.Tanh(pressure);
}
}
}

34
Core/ThreadLoadTracker.cs
View File

@@ -1,34 +0,0 @@
using System;
using System.Threading;
namespace FluidSim
{
/// <summary>
/// Tracks the duty cycle of a worker thread using an exponential moving average.
/// Thread‑safe: one writer (the sim thread), any reader (UI thread).
/// </summary>
public class ThreadLoadTracker
{
private double _loadPercent; // 0 .. 100, accessed with Volatile.Read/Write
private const double Alpha = 0.1; // smoothing factor (higher = faster response)
/// <summary>
/// Update the load percentage with a new observation.
/// </summary>
/// <param name="busyMs">Time spent on real work in the last cycle.</param>
/// <param name="totalMs">Total time of the last cycle (work + idle). If zero, ignored.</param>
public void Record(double busyMs, double totalMs)
{
if (totalMs <= 0) return;
double instantLoad = busyMs / totalMs * 100.0;
// Exponential moving average
double old = Volatile.Read(ref _loadPercent);
double newLoad = old + Alpha * (instantLoad - old);
Volatile.Write(ref _loadPercent, newLoad);
}
/// <summary>Current smoothed load percentage (0‑100).</summary>
public double LoadPercent => Volatile.Read(ref _loadPercent);
}
}

11
FluidSim.csproj
View File

@@ -5,7 +5,7 @@
<TargetFramework>net10.0</TargetFramework>
<ImplicitUsings>enable</ImplicitUsings>
<Nullable>enable</Nullable>
<PublishAot>false</PublishAot>
<PublishAot>true</PublishAot>
<InvariantGlobalization>true</InvariantGlobalization>
</PropertyGroup>
@@ -13,13 +13,4 @@
<PackageReference Include="SFML.Net" Version="3.0.0" />
</ItemGroup>
<ItemGroup>
<None Update="fonts\FiraCodeNerdFont-Medium.ttf">
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Update="fonts\LiberationMono-Regular.ttf">
<CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
</None>
</ItemGroup>
</Project>

10
Interfaces/IComponent.cs
View File

@@ -1,10 +0,0 @@
using System.Collections.Generic;
namespace FluidSim.Interfaces
{
public interface IComponent
{
IReadOnlyList<Port> Ports { get; }
void UpdateState(float dt);
}
}

24
Interfaces/Port.cs
View File

@@ -2,23 +2,19 @@
{
public class Port
{
public float MassFlowRate; // kg/s, positive INTO owning component
public float SpecificEnthalpy; // J/kg
public float Pressure; // Pa
public float Density; // kg/m³
public float Temperature; // K
public float AirFraction; // mass fraction (0 = exhaust, 1 = air)
public object? Owner { get; set; }
public double Pressure; // Pa
public double MassFlowRate; // kg/s, positive INTO the component
public double SpecificEnthalpy; // J/kg, enthalpy of fluid entering this port
public double Density; // kg/m³
public double Temperature; // K
public Port()
{
MassFlowRate = 0f;
SpecificEnthalpy = 0f;
Pressure = 101325f;
Density = 1.225f;
Temperature = 300f;
AirFraction = 1f;
Pressure = 101325.0;
MassFlowRate = 0.0;
SpecificEnthalpy = 0.0;
Density = 1.225;
Temperature = 300.0;
}
}
}

242
Program.cs
View File

@@ -1,216 +1,96 @@
using FluidSim.Audio;
using FluidSim.Core;
using FluidSim.Tests;
using SFML.Graphics;
using SFML.System;
using SFML.Graphics;
using SFML.Window;
using System;
using SFML.System;
using System.Diagnostics;
using System.Threading;
using System.Threading.Tasks;
using FluidSim.Scenarios;
using FluidSim.Audio;
namespace FluidSim;
public class Program
{
private const int SampleRate = 44100;
private const double DrawFrequency = 60.0;
// Playback speed
private static double _desiredSpeed = 0.001;
private static double _currentDisplaySpeed = _desiredSpeed;
private const double MinSpeed = 0.0001;
private const double MaxSpeed = 1.0;
private const double ScrollFactor = 1.1;
private static double _lastNormalSpeed = 0.1;
private static bool _isRealTime = false;
private static volatile bool _timeWarpActive;
// Thread load tracking
private static ThreadLoadTracker _loadTracker = new ThreadLoadTracker();
// Audio & simulation
private static SimulationRingBuffer _simRingBuffer = null!;
private static SoundEngine _soundEngine = null!;
private static Scenario _scenario = null!; // cast to access ThrottleArea
private static Font? _overlayFont;
private static Text? _overlayText;
// Throttle control
private static float _throttleTarget = 1.0f; // 0‑1, set by arrow keys
private static float _throttleCurrent = 0.0f; // actual current fraction (lerped)
private const float ThrottleLerpRate = 10.0f; // times per second (speed of movement)
private static bool _wKeyHeld = false;
private static float _lastThrottleUpdateTime;
private const int TargetMaxFill = (int)(SampleRate * 0.2);
private static volatile bool running = true;
// Global step counter – incremented every simulation step
private static long stepCount = 0;
public static void Main()
{
var window = CreateWindow();
LoadFont();
_scenario = new SingleCylScenario();
_scenario.Initialize(SampleRate);
_lastThrottleUpdateTime = 0.0f;
var mode = new VideoMode(new Vector2u(1280, 720));
var window = new RenderWindow(mode, "Fluid Simulation");
window.SetVerticalSyncEnabled(true);
window.Closed += (_, _) => { running = false; window.Close(); };
_simRingBuffer = new SimulationRingBuffer(8192);
_soundEngine = new SoundEngine(_simRingBuffer) { Volume = 100 };
_soundEngine.Start();
var cts = new CancellationTokenSource();
Task.Run(() => SimulationLoop(cts.Token), cts.Token);
var soundEngine = new SoundEngine(bufferCapacity: 2048);
soundEngine.Volume = 70;
soundEngine.Start();
var stopwatch = Stopwatch.StartNew();
double lastDrawTime = 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);
// 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];
double lastLogTime = 0.0; // for periodic speed printout
while (window.IsOpen)
{
window.DispatchEvents();
double now = stopwatch.Elapsed.TotalSeconds;
// --- 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)
// ---- Playback speed smoothing ----
double targetSpeed = _timeWarpActive ? 1.0 : _desiredSpeed;
_currentDisplaySpeed += (targetSpeed - _currentDisplaySpeed) *
(1.0 - Math.Exp(-8.0 * (now - lastDrawTime)));
_soundEngine.Speed = _currentDisplaySpeed;
// ---- Throttle update ----
float dtThrottle = (float)now - _lastThrottleUpdateTime;
_lastThrottleUpdateTime = (float)now;
float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f;
// Snap to zero instantly when target is zero (key released)
if (throttleDesiredFraction == 0.0)
// --- Generate the required number of simulation steps ---
while (samplesNeeded > 0 && running)
{
_throttleCurrent = 0.0f;
int toGenerate = Math.Min(samplesNeeded, chunkSize);
for (int i = 0; i < toGenerate; i++)
{
buffer[i] = Simulation.Process();
stepCount++;
}
else
{
float smoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle);
_throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing;
soundEngine.WriteSamples(buffer, toGenerate);
samplesNeeded -= toGenerate;
}
_scenario.Throttle = _throttleCurrent;
// --- Display speed ---
double simTime = stepCount / (double)SampleRate;
double wallTime = stopwatch.Elapsed.TotalSeconds;
double speed = (wallTime > 0) ? simTime / wallTime : 0.0;
// ---- Drawing ----
if (now - lastDrawTime >= 1.0 / DrawFrequency)
// Update window title with instant speed
window.SetTitle($"FluidSim | Speed: {speed:F3}× | Steps: {stepCount}");
// Console log once per second
if (wallTime - lastLogTime >= 1.0)
{
if (_overlayText != null)
{
string toggleHint = _isRealTime ? "[Space] slow mo" : "[Space] real time";
_overlayText.DisplayedString =
$"{toggleHint} Speed: {_currentDisplaySpeed:F3}x RT: {(_currentDisplaySpeed * 100.0):F1}% Sim load: {_loadTracker.LoadPercent:F0}%\n" +
$"Throttle: {_throttleCurrent * 100:F0}% Target: {_throttleTarget * 100:F0}% [W] {(_wKeyHeld ? "BLIP" : "---")}";
Console.WriteLine($"Speed: {speed:F3}× ({stepCount} steps, {wallTime:F2}s wall)");
lastLogTime = wallTime;
}
// --- Rendering (placeholder) ---
window.Clear(Color.Black);
_scenario.Draw(window);
if (_overlayText != null) window.Draw(_overlayText);
window.Display();
lastDrawTime = now;
}
}
cts.Cancel();
_soundEngine.Dispose();
// --- Cleanup ---
soundEngine.Dispose();
window.Dispose();
}
private static void SimulationLoop(CancellationToken token)
{
while (!token.IsCancellationRequested)
{
long cycleStart = Stopwatch.GetTimestamp();
long workStart = Stopwatch.GetTimestamp();
float sample = _scenario.Process();
_simRingBuffer.Write(sample);
long workEnd = Stopwatch.GetTimestamp();
while (_simRingBuffer.AvailableSamples > TargetMaxFill &&
!token.IsCancellationRequested)
{
Thread.Sleep(1);
}
long cycleEnd = Stopwatch.GetTimestamp();
double busyMs = (workEnd - workStart) / (double)Stopwatch.Frequency * 1000.0;
double totalMs = (cycleEnd - cycleStart) / (double)Stopwatch.Frequency * 1000.0;
_loadTracker.Record(busyMs, totalMs);
}
}
// ---------- Window & input ----------
private static RenderWindow CreateWindow()
{
var mode = new VideoMode(new Vector2u(1280, 720));
var window = new RenderWindow(mode, "FluidSim");
window.SetVerticalSyncEnabled(false);
window.SetFramerateLimit(60);
window.Closed += (_, _) => window.Close();
window.MouseWheelScrolled += OnMouseWheel;
window.KeyPressed += OnKeyPressed;
window.KeyReleased += OnKeyReleased;
return window;
}
private static void LoadFont()
{
try { _overlayFont = new Font("fonts/FiraCodeNerdFont-Medium.ttf"); }
catch { _overlayFont = null; }
if (_overlayFont != null)
_overlayText = new Text(_overlayFont)
{
FillColor = Color.White,
Position = new Vector2f(10, 10)
};
}
private static void OnMouseWheel(object? sender, MouseWheelScrollEventArgs e)
{
if (_timeWarpActive) return;
if (e.Delta > 0) _desiredSpeed *= ScrollFactor;
else if (e.Delta < 0) _desiredSpeed /= ScrollFactor;
_desiredSpeed = Math.Clamp(_desiredSpeed, MinSpeed, MaxSpeed);
_lastNormalSpeed = _desiredSpeed;
_isRealTime = Math.Abs(_desiredSpeed - 1.0) < 1e-6;
}
private static void OnKeyPressed(object? sender, KeyEventArgs e)
{
switch (e.Code)
{
case Keyboard.Key.Space:
_timeWarpActive = !_timeWarpActive;
if (!_timeWarpActive)
{
_desiredSpeed = _lastNormalSpeed;
_isRealTime = false;
}
break;
case Keyboard.Key.W:
_wKeyHeld = true;
break;
case Keyboard.Key.Up:
_throttleTarget = MathF.Min(1.0f, _throttleTarget + 0.05f);
break;
case Keyboard.Key.Down:
_throttleTarget = MathF.Max(0.0f, _throttleTarget - 0.05f);
break;
}
}
private static void OnKeyReleased(object? sender, KeyEventArgs e)
{
if (e.Code == Keyboard.Key.W)
_wKeyHeld = false;
}
}

BIN
Report20260507-1444 (2).diagsession
View File

Binary file not shown.

118
Scenarios/HelmholtzScenario.cs
View File

@@ -1,118 +0,0 @@
using FluidSim.Components;
using FluidSim.Core;
using FluidSim.Interfaces;
using SFML.Graphics;
using SFML.System;
using System;
namespace FluidSim.Tests
{
public class HelmholtzScenario : Scenario
{
private Volume0D cavity;
private Port cavityPort;
private PipeSystem pipeSystem;
private int[] pipeStart = { 0 };
private int[] pipeEnd;
private BoundarySystem boundaries;
private int cavityOrificeIdx = 0;
private int openEndIdx = 0;
private Solver solver;
private double dt;
private int stepCount;
private SoundProcessor soundProcessor;
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
// --- Realistic Helmholtz resonator dimensions ---
float cavityVolume = 1e-3f; // 1 liter
float neckLength = 0.05f; // 5 cm
float neckDiameter = 0.02f; // 2 cm diameter
float neckArea = MathF.PI * 0.25f * neckDiameter * neckDiameter;
int neckCells = 20;
// --- Volume (cavity) ---
float initialPressure = 1.1f * 101325f; // slight overpressure
float initialTemperature = 300f;
cavity = new Volume0D(cavityVolume, initialPressure, initialTemperature);
cavityPort = cavity.CreatePort();
// --- Pipe (neck) ---
float[] areas = new float[neckCells];
float[] dxs = new float[neckCells];
float dx = neckLength / neckCells;
for (int i = 0; i < neckCells; i++)
{
areas[i] = neckArea;
dxs[i] = dx;
}
pipeEnd = new[] { neckCells };
float rho0 = 101325f / (287f * 300f);
pipeSystem = new PipeSystem(neckCells, pipeStart, pipeEnd, areas, dxs, rho0, 0f, 101325f);
// --- Boundary system ---
boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
// Standard orifice with built‑in minor loss (K = 0.5) – no inertance needed
boundaries.AddOrificeWithInertance(
cavityPort, pipeIndex: 0, isLeftEnd: true,
areaIndex: cavityOrificeIdx,
dischargeCoeff: 0.9f,
effectiveLength: neckLength // physical neck length
);
// Open end at right side of pipe
boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, neckArea);
float[] orificeAreas = new float[1] { neckArea };
boundaries.SetOrificeAreas(orificeAreas);
// --- Solver ---
solver = new Solver { SubStepCount = 8, EnableProfiling = false };
solver.SetTimeStep(dt);
solver.SetPipeSystem(pipeSystem);
solver.SetBoundarySystem(boundaries);
solver.AddComponent(cavity);
// --- Sound ---
soundProcessor = new SoundProcessor(sampleRate, 1f) { Gain = 2f };
Console.WriteLine("Helmholtz resonator ready.");
stepCount = 0;
}
public override float Process()
{
solver.Step();
stepCount++;
float flow = boundaries.GetOpenEndMassFlow(openEndIdx);
float sample = soundProcessor.Process(flow);
return sample;
}
public override void Draw(RenderWindow target)
{
float winW = target.GetView().Size.X;
float winH = target.GetView().Size.Y;
float cavityCenterX = 100f;
float cavityWidth = 80f, cavityHeight = 100f;
float cavityTopY = winH / 2f - cavityHeight / 2f;
DrawVolume(target, cavity, cavityCenterX, cavityTopY - 40f, cavityWidth, cavityHeight);
float pipeStartX = cavityCenterX + cavityWidth / 2f + 10f;
float pipeEndX = winW - 50f;
float pipeCenterY = winH / 2f;
DrawPipe(target, pipeSystem, 0, pipeCenterY, pipeStartX, pipeEndX);
}
}
}

161
Scenarios/Scenario.cs
View File

@@ -1,161 +0,0 @@
using SFML.Graphics;
using SFML.System;
using FluidSim.Core;
using FluidSim.Components;
namespace FluidSim.Tests
{
public abstract class Scenario
{
protected const float AmbientPressure = 101325f;
protected const float AmbientTemperature = 300f;
public float Throttle { get; set; }
public abstract void Initialize(int sampleRate);
public abstract float Process();
public abstract void Draw(RenderWindow target);
protected Color PressureColor(float pressurePa)
{
float bar = pressurePa / 1e5f;
byte r, g, b;
if (bar < 1f)
{
float f = Math.Clamp(bar, 0f, 1f);
r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
}
else
{
float f = Math.Min((bar - 1f) / 9f, 1f);
r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
}
return new Color(r, g, b);
}
protected Color TemperatureColor(float t)
{
t = Math.Clamp(t, 0f, 2000f);
byte r, g, b;
if (t < AmbientTemperature)
{
float f = t / AmbientTemperature;
r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f));
}
else
{
float f = (t - AmbientTemperature) / (2000f - AmbientTemperature);
r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0;
}
return new Color(r, g, b);
}
protected void DrawVolume(RenderWindow target, Volume0D volume,
float centerX, float topY, float width, float height)
{
var rect = new RectangleShape(new Vector2f(width, height))
{
FillColor = PressureColor(volume.Pressure),
Position = new Vector2f(centerX - width / 2f, topY)
};
target.Draw(rect);
var border = new RectangleShape(new Vector2f(width, height))
{
FillColor = Color.Transparent,
OutlineColor = Color.White,
OutlineThickness = 1f,
Position = rect.Position
};
target.Draw(border);
}
protected void DrawCylinder(RenderWindow target, Cylinder cylinder,
float centerX, float topY, float width, float maxHeight)
{
float fraction = cylinder.PistonFraction;
float currentHeight = maxHeight * fraction;
var wall = new RectangleShape(new Vector2f(width, maxHeight))
{
FillColor = new Color(60, 60, 60),
Position = new Vector2f(centerX - width / 2f, topY)
};
target.Draw(wall);
var gas = new RectangleShape(new Vector2f(width, currentHeight))
{
FillColor = PressureColor(cylinder.Pressure),
Position = new Vector2f(centerX - width / 2f, topY)
};
target.Draw(gas);
var piston = new RectangleShape(new Vector2f(width, 4f))
{
FillColor = Color.White,
Position = new Vector2f(centerX - width / 2f, topY + currentHeight)
};
target.Draw(piston);
float valveW = 6f, valveH = 10f, valveY = topY + 4f;
var iv = new RectangleShape(new Vector2f(valveW, valveH))
{
FillColor = cylinder.IntakeValveArea > 0f ? Color.Green : Color.Red,
Position = new Vector2f(centerX - width / 2f - valveW - 2f, valveY)
};
target.Draw(iv);
var ev = new RectangleShape(new Vector2f(valveW, valveH))
{
FillColor = cylinder.ExhaustValveArea > 0f ? Color.Green : Color.Red,
Position = new Vector2f(centerX + width / 2f + 2f, valveY)
};
target.Draw(ev);
}
protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex,
float pipeCenterY, float pipeStartX, float pipeEndX)
{
int start = pipeSystem.GetPipeStart(pipeIndex);
int end = pipeSystem.GetPipeEnd(pipeIndex);
int n = end - start;
if (n < 2) return;
float pipeLen = pipeEndX - pipeStartX;
float dx = pipeLen / (n - 1);
float baseRadius = 25f;
var centers = new float[n];
var radii = new float[n];
var temps = new float[n];
for (int i = 0; i < n; i++)
{
int cell = start + i;
float p = pipeSystem.GetCellPressure(cell);
float rho = pipeSystem.GetCellDensity(cell);
temps[i] = p / MathF.Max(rho * 287f, 1e-12f);
float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f);
radii[i] = baseRadius * (1f + dev * 2f);
if (radii[i] < 2f) radii[i] = 2f;
centers[i] = pipeStartX + i * dx;
}
int segments = 8;
var va = new VertexArray(PrimitiveType.TriangleStrip);
for (int i = 0; i < n; i++)
{
float x = centers[i], r = radii[i];
Color col = TemperatureColor(temps[i]);
va.Append(new Vertex(new Vector2f(x, pipeCenterY - r), col));
va.Append(new Vertex(new Vector2f(x, pipeCenterY + r), col));
if (i < n - 1)
{
for (int s = 1; s <= segments; s++)
{
float t = s / (float)segments;
float xi = centers[i] + (centers[i + 1] - centers[i]) * t;
float ri = radii[i] + (radii[i + 1] - radii[i]) * t;
float Ti = temps[i] + (temps[i + 1] - temps[i]) * t;
Color colS = TemperatureColor(Ti);
va.Append(new Vertex(new Vector2f(xi, pipeCenterY - ri), colS));
va.Append(new Vertex(new Vector2f(xi, pipeCenterY + ri), colS));
}
}
}
target.Draw(va);
}
}
}

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

246
Scenarios/SingleCylScenario.cs
View File

@@ -1,246 +0,0 @@
using FluidSim.Components;
using FluidSim.Core;
using FluidSim.Interfaces;
using FluidSim.Utils;
using SFML.Graphics;
using SFML.System;
using System;
namespace FluidSim.Tests
{
public class SingleCylScenario : Scenario
{
private Crankshaft crankshaft;
private Cylinder cylinder;
private PipeSystem pipeSystem;
private BoundarySystem boundaries;
private Solver solver;
private Volume0D intakePlenum;
private Port plenumInlet, plenumOutlet;
private Volume0D exhaustCollector;
private Port colIn, colOut;
private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx;
private float[] orificeAreas;
private int intakeOpenIdx, exhaustOpenIdx;
private SoundProcessor exhaustSound, intakeSound;
private OutdoorExhaustReverb reverb;
private double dt;
private int stepCount;
// Use a private field for the maximum throttle area, avoiding any base‑class conflicts
private float _maxThrottleArea;
// pipe area for open end calculations
private float pipeArea;
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
// Maximum throttle area – independent of base class
_maxThrottleArea = (float)Units.AreaFromDiameter(3 * Units.cm); // 1 cm²
// ---- Crankshaft ----
crankshaft = new Crankshaft(2000);
crankshaft.Inertia = 0.01f;
crankshaft.FrictionConstant = 2f;
crankshaft.FrictionViscous = 0.0f;
// ---- Cylinder ----
float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 11f;
float ivo = 350f, ivc = 580f, evo = 120f, evc = 370f;
cylinder = new Cylinder(bore, stroke, conRod, compRatio,
ivo, ivc, evo, evc, crankshaft)
{
IntakeValveDiameter = 0.03f,
IntakeValveLift = 0.005f,
ExhaustValveDiameter = 0.028f,
ExhaustValveLift = 0.005f
};
// ---- Pipe system ----
int[] pipeStart = { 0, 10, 20 };
int[] pipeEnd = { 10, 20, 70 };
int totalCells = pipeEnd[^1]; // automatically 70, stays in sync
float[] area = new float[totalCells];
float[] dx = new float[totalCells];
float pipeDiameter = 0.02f; // 2 cm
pipeArea = MathF.PI * 0.25f * pipeDiameter * pipeDiameter;
float areaVal = pipeArea;
float intakeLenBefore = 0.2f, intakeLenRunner = 0.2f, exhaustLen = 0.4f;
for (int i = 0; i < totalCells; i++)
{
area[i] = areaVal;
if (i < 10) dx[i] = intakeLenBefore / 10f;
else if (i < 20) dx[i] = intakeLenRunner / 10f;
else dx[i] = exhaustLen / 50f;
}
pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx,
1.225f, 0f, 101325f);
pipeSystem.DampingMultiplier = 1.0f;
pipeSystem.EnergyRelaxationRate = 0.5f;
pipeSystem.AmbientPressure = 101325f;
// ---- Volumes ----
intakePlenum = new Volume0D(100e-6f, 101325f, 300f); // 100 mL
plenumInlet = intakePlenum.CreatePort();
plenumOutlet = intakePlenum.CreatePort();
exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // 10 mL (unused but present)
colIn = exhaustCollector.CreatePort();
colOut = exhaustCollector.CreatePort();
// ---- Boundary system ----
boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2);
throttleAreaIdx = 0;
plenumRunnerAreaIdx = 1;
intakeValveIdx = 2;
exhaustValveIdx = 3;
// Intake open end (pipe0 left)
boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, pipeArea);
intakeOpenIdx = 0;
// Throttle orifice (plenum inlet to pipe0 right)
boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.2f);
// Plenum to runner (plenum outlet to pipe1 left)
boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1f);
// Intake valve (cylinder intake to pipe1 right)
boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1f);
// Exhaust valve (cylinder exhaust to pipe2 left)
boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1f);
// Exhaust open end (pipe2 right)
boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, pipeArea);
exhaustOpenIdx = 1;
orificeAreas = new float[4];
orificeAreas[plenumRunnerAreaIdx] = areaVal; // fixed plenum->runner area
// ---- Solver ----
solver = new Solver { SubStepCount = 4, EnableProfiling = false };
solver.SetTimeStep(dt);
solver.SetPipeSystem(pipeSystem);
solver.SetBoundarySystem(boundaries);
solver.AddComponent(cylinder);
solver.AddComponent(intakePlenum);
solver.AddComponent(exhaustCollector);
// ---- Sound ----
exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 20f };
reverb = new OutdoorExhaustReverb(sampleRate);
stepCount = 0;
Console.WriteLine("TestScenario ready.");
}
public override float Process()
{
crankshaft.Step((float)dt);
cylinder.PreStep((float)dt);
// Update variable orifice areas – use the private _maxThrottleArea
float throttledArea = _maxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1f);
orificeAreas[throttleAreaIdx] = throttledArea;
orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea;
orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea;
boundaries.SetOrificeAreas(orificeAreas);
solver.Step();
stepCount++;
// Retrieve open‑end mass flows for sound synthesis
float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx);
float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx);
float exhaustDry = exhaustSound.Process(exhaustFlow);
float intakeDry = intakeSound.Process(intakeFlow);
if (stepCount % 1000 == 0)
{
float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI);
float crankDeg = crankshaft.CrankAngle; // degrees (0–720)
Console.WriteLine($"Step {stepCount}, CA={crankDeg:F1} deg, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar");
Console.WriteLine($" intake flow: {intakeFlow:F6}, exhaust flow: {exhaustFlow:F6}");
// Pipe 0 (intake before throttle)
var (r0L, u0L, p0L) = pipeSystem.GetInteriorStateLeft(0);
var (r0R, u0R, p0R) = pipeSystem.GetInteriorStateRight(0);
Console.WriteLine($" Pipe0 L: rho={r0L:F4} u={u0L:F3} p={p0L/1e5:F3}bar | R: rho={r0R:F4} u={u0R:F3} p={p0R/1e5:F3}bar");
// Pipe 1 (runner)
var (r1L, u1L, p1L) = pipeSystem.GetInteriorStateLeft(1);
var (r1R, u1R, p1R) = pipeSystem.GetInteriorStateRight(1);
Console.WriteLine($" Pipe1 L: rho={r1L:F4} u={u1L:F3} p={p1L/1e5:F3}bar | R: rho={r1R:F4} u={u1R:F3} p={p1R/1e5:F3}bar");
// Pipe 2 (exhaust)
var (r2L, u2L, p2L) = pipeSystem.GetInteriorStateLeft(2);
var (r2R, u2R, p2R) = pipeSystem.GetInteriorStateRight(2);
Console.WriteLine($" Pipe2 L: rho={r2L:F4} u={u2L:F3} p={p2L/1e5:F3}bar | R: rho={r2R:F4} u={u2R:F3} p={p2R/1e5:F3}bar");
// Plenum and cylinder mass
Console.WriteLine($" Plenum P={intakePlenum.Pressure/1e5:F3}bar, mass={intakePlenum.Mass:E4} kg");
Console.WriteLine($" Cyl mass={cylinder.Mass:E4} kg");
}
return reverb.Process(intakeDry + exhaustDry);
}
public override void Draw(RenderWindow target)
{
float winW = target.GetView().Size.X;
float winH = target.GetView().Size.Y;
float intakeY = winH / 2f - 40f;
float exhaustY = winH / 2f + 80f;
float openEndX = 40f;
// Intake pipe before throttle (pipe 0)
float pipe1StartX = openEndX;
float pipe1EndX = pipe1StartX + 120f;
DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX);
// Throttle symbol
float throttleX = pipe1EndX + 5f;
var throttleRect = new RectangleShape(new Vector2f(8f, 30f))
{
FillColor = Color.Yellow,
Position = new Vector2f(throttleX, intakeY - 15f)
};
target.Draw(throttleRect);
// Plenum
float plenW = 60f, plenH = 80f;
float plenLeftX = throttleX + 10f;
float plenCenterX = plenLeftX + plenW / 2f;
float plenTopY = intakeY - plenH / 2f;
DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH);
// Runner pipe (pipe 1)
float runnerStartX = plenLeftX + plenW + 5f;
float runnerEndX = runnerStartX + 100f;
DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX);
// Cylinder
float cylCX = runnerEndX + 50f;
float cylTopY = intakeY - 120f;
float cylW = 80f, cylMaxH = 240f;
DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH);
// Exhaust pipe (pipe 2)
float exhStartX = cylCX + cylW / 2f + 20f;
float exhEndX = winW - 60f;
DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX);
}
}
}

102
Scenarios/TestScenario.cs
View File

@@ -1,102 +0,0 @@
using System;
using SFML.Graphics;
using SFML.System;
using FluidSim.Core;
namespace FluidSim.Tests
{
public class TestScenario : Scenario
{
private PipeSystem pipeSystem;
private BoundarySystem boundaries;
private Solver solver;
private int[] pipeStart = { 0 };
private int[] pipeEnd;
private double dt;
private int stepCount;
// Sound output: use pressure at open end
private SoundProcessor openEndSound;
private int openEndIdx = 0; // index of the open end in BoundarySystem (we added only one)
public override void Initialize(int sampleRate)
{
dt = 1.0 / sampleRate;
const int cellCount = 200;
float length = 2f;
float dia = 0.02f;
float area = MathF.PI * 0.25f * dia * dia;
float[] areas = new float[cellCount];
float[] dxs = new float[cellCount];
float dx = length / cellCount;
for (int i = 0; i < cellCount; i++)
{
areas[i] = area;
dxs[i] = dx;
}
pipeEnd = new[] { cellCount };
float rho0 = 101325f / (287f * 300f);
pipeSystem = new PipeSystem(cellCount, pipeStart, pipeEnd, areas, dxs,
rho0, 0f, 101325f);
pipeSystem.DampingMultiplier = 0f;
pipeSystem.EnergyRelaxationRate = 0f;
pipeSystem.AmbientPressure = 101325f;
// Pressure bubble near right end
float pBubble = 10f * 101325f;
float TBubble = 2000f;
float rhoBubble = pBubble / (287f * TBubble);
for (int i = 0; i <= 10; i++)
pipeSystem.SetCellState(i, rhoBubble, 0f, pBubble);
// Boundaries: left closed, right open
boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);
boundaries.AddOrifice(null, pipeIndex: 0, isLeftEnd: true, areaIndex: 0, 1f);
boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, area);
float[] orificeAreas = new float[1] { 0f };
boundaries.SetOrificeAreas(orificeAreas);
solver = new Solver { SubStepCount = 3};
solver.SetTimeStep(dt);
solver.SetPipeSystem(pipeSystem);
solver.SetBoundarySystem(boundaries);
solver.EnableProfiling = true;
pipeSystem.EnableProfiling = true;
// Simple sound processor: convert mass flow rate to audio
openEndSound = new SoundProcessor(sampleRate, 1f) { Gain = 2f };
Console.WriteLine("Pulse test ready.");
stepCount = 0;
}
public override float Process()
{
solver.Step();
stepCount++;
float flow = boundaries.GetOpenEndMassFlow(openEndIdx);
float sample = openEndSound.Process(flow);
return sample;
}
public override void Draw(RenderWindow target)
{
float winW = target.GetView().Size.X;
float winH = target.GetView().Size.Y;
float startX = 50f;
float endX = winW - 50f;
float y = winH / 2f;
DrawPipe(target, pipeSystem, 0, y, startX, endX);
}
}
}

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

BIN
fonts/FiraCodeNerdFont-Medium.ttf
View File

Binary file not shown.

1
trace.speedscope.speedscope.json
View File

File diff suppressed because one or more lines are too long