450 lines
18 KiB
C#
450 lines
18 KiB
C#
using System;
|
||
using System.Diagnostics;
|
||
using FluidSim.Interfaces;
|
||
|
||
namespace FluidSim.Components
|
||
{
|
||
/// <summary>
|
||
/// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme.
|
||
/// Ghost states are set externally via SetGhostLeft/Right; they are always required.
|
||
/// Now includes a passive scalar for air mass fraction.
|
||
/// </summary>
|
||
public class Pipe1D : IComponent
|
||
{
|
||
public const bool EnableDetailedProfiling = false; // set to false in release builds
|
||
|
||
public Port PortA { get; }
|
||
public Port PortB { get; }
|
||
public double Area { get; }
|
||
public double DampingMultiplier { get; set; } = 10.0;
|
||
public double EnergyRelaxationRate { get; set; } = 5.0; // 1/s
|
||
|
||
private double _ambientPressure = 101325.0;
|
||
public double AmbientPressure
|
||
{
|
||
get => _ambientPressure;
|
||
set
|
||
{
|
||
_ambientPressure = value;
|
||
_ambientEnergyReference = value / (_gamma - 1.0);
|
||
}
|
||
}
|
||
|
||
private readonly int _n;
|
||
private readonly double _dx;
|
||
private readonly double _diameter;
|
||
private readonly double _gamma = 1.4;
|
||
|
||
private double[] _rho, _rhou, _E;
|
||
private double[] _Y; // air mass fraction
|
||
private double[] _fluxM, _fluxP, _fluxE;
|
||
|
||
private double _rhoGhostL, _uGhostL, _pGhostL, _YGhostL;
|
||
private double _rhoGhostR, _uGhostR, _pGhostR, _YGhostR;
|
||
private bool _ghostLValid, _ghostRValid;
|
||
|
||
private double _laminarCoeff;
|
||
private double _ambientEnergyReference;
|
||
|
||
// ---------- Profiling accumulators ----------
|
||
private long _profPrecomputeTicks;
|
||
private long _profLeftFluxTicks;
|
||
private long _profInteriorLoopTicks;
|
||
private long _profRightFluxTicks;
|
||
private long _profPortUpdateTicks;
|
||
private long _profCallCount;
|
||
|
||
public Pipe1D(double length, double area, int cellCount)
|
||
{
|
||
if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4");
|
||
|
||
_n = cellCount;
|
||
_dx = length / _n;
|
||
Area = area;
|
||
_diameter = 2.0 * Math.Sqrt(area / Math.PI);
|
||
|
||
_rho = new double[_n];
|
||
_rhou = new double[_n];
|
||
_E = new double[_n];
|
||
_Y = new double[_n];
|
||
_fluxM = new double[_n + 1];
|
||
_fluxP = new double[_n + 1];
|
||
_fluxE = new double[_n + 1];
|
||
|
||
double mu_air = 1.8e-5;
|
||
double radius = _diameter * 0.5;
|
||
_laminarCoeff = 8.0 * mu_air / (radius * radius);
|
||
|
||
_ambientEnergyReference = 101325.0 / (_gamma - 1.0);
|
||
|
||
PortA = new Port { Owner = this };
|
||
PortB = new Port { Owner = this };
|
||
|
||
SetUniformState(1.225, 0.0, 101325.0);
|
||
}
|
||
|
||
IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };
|
||
|
||
public void UpdateState(double dt) { }
|
||
|
||
// ---------- Ghost interface ----------
|
||
public void SetGhostLeft(double rho, double u, double p, double airFraction)
|
||
{
|
||
_rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _YGhostL = airFraction; _ghostLValid = true;
|
||
}
|
||
public void SetGhostRight(double rho, double u, double p, double airFraction)
|
||
{
|
||
_rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _YGhostR = airFraction; _ghostRValid = true;
|
||
}
|
||
public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; }
|
||
|
||
public double GetInteriorAirFractionLeft() => _Y[0];
|
||
public double GetInteriorAirFractionRight() => _Y[_n - 1];
|
||
|
||
public (double rho, double u, double p) GetInteriorStateLeft()
|
||
{
|
||
double rho = Math.Max(_rho[0], 1e-12);
|
||
double u = _rhou[0] / rho;
|
||
double p = PressureScalar(0);
|
||
return (rho, u, p);
|
||
}
|
||
public (double rho, double u, double p) GetInteriorStateRight()
|
||
{
|
||
double rho = Math.Max(_rho[_n - 1], 1e-12);
|
||
double u = _rhou[_n - 1] / rho;
|
||
double p = PressureScalar(_n - 1);
|
||
return (rho, u, p);
|
||
}
|
||
public int CellCount => _n;
|
||
public double GetCellDensity(int i) => _rho[i];
|
||
public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
|
||
public double GetCellPressure(int i) => PressureScalar(i);
|
||
|
||
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
|
||
{
|
||
double maxW = 0.0;
|
||
for (int i = 0; i < _n; i++)
|
||
{
|
||
double rho = Math.Max(_rho[i], 1e-12);
|
||
double u = Math.Abs(_rhou[i] / rho);
|
||
double p = PressureScalar(i);
|
||
double c = Math.Sqrt(_gamma * p / rho);
|
||
double local = u + c;
|
||
if (local > maxW) maxW = local;
|
||
}
|
||
maxW = Math.Max(maxW, 1e-8);
|
||
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
|
||
}
|
||
|
||
// ---------- Main step (per sub‑step) ----------
|
||
public void SimulateSingleStep(double dtSub)
|
||
{
|
||
if (!_ghostLValid || !_ghostRValid)
|
||
throw new InvalidOperationException("Ghost cells not set before SimulateSingleStep.");
|
||
|
||
double dt = dtSub;
|
||
int n = _n;
|
||
double dt_dx = dt / _dx;
|
||
double coeff = _laminarCoeff * DampingMultiplier;
|
||
double relaxRate = EnergyRelaxationRate;
|
||
double gamma = _gamma;
|
||
double gm1 = gamma - 1.0;
|
||
|
||
// ---------- Profiling start ----------
|
||
long t0 = 0, t1 = 0;
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t0 = Stopwatch.GetTimestamp();
|
||
_profCallCount++;
|
||
}
|
||
|
||
// ---------- Phase 1: Pre‑compute pressure and speed of sound ----------
|
||
double[] p = new double[n];
|
||
double[] c = new double[n];
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
double rho = Math.Max(_rho[i], 1e-12);
|
||
double u = _rhou[i] / rho;
|
||
p[i] = gm1 * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
|
||
c[i] = Math.Sqrt(gamma * p[i] / rho);
|
||
}
|
||
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profPrecomputeTicks += (t1 - t0);
|
||
t0 = t1;
|
||
}
|
||
|
||
// ---------- Local flux functions ----------
|
||
void LaxFlux(double rL, double uL, double pL, double cL,
|
||
double rR, double uR, double pR, double cR,
|
||
out double fm, out double fp, out double fe)
|
||
{
|
||
double EL = pL / (gm1 * rL) + 0.5 * uL * uL;
|
||
double ER = pR / (gm1 * rR) + 0.5 * uR * uR;
|
||
double Fm_L = rL * uL;
|
||
double Fp_L = rL * uL * uL + pL;
|
||
double Fe_L = (rL * EL + pL) * uL;
|
||
double Fm_R = rR * uR;
|
||
double Fp_R = rR * uR * uR + pR;
|
||
double Fe_R = (rR * ER + pR) * uR;
|
||
double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
|
||
fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rR - rL);
|
||
fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rR * uR - rL * uL);
|
||
fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (rR * ER - rL * EL);
|
||
}
|
||
|
||
void ScalarFlux(double rL, double uL, double cL, double YL,
|
||
double rR, double uR, double cR, double YR,
|
||
double alpha, out double fy)
|
||
{
|
||
double Fm_L = rL * uL;
|
||
double Fm_R = rR * uR;
|
||
fy = 0.5 * (Fm_L * YL + Fm_R * YR) - 0.5 * alpha * (rR * YR - rL * YL);
|
||
}
|
||
|
||
// ---------- Phase 2: Left face flux (ghostL – cell 0) ----------
|
||
double rL_ghost = Math.Max(_rhoGhostL, 1e-12);
|
||
double pL_ghost = _pGhostL;
|
||
double uL_ghost = _uGhostL;
|
||
double cL_ghost = Math.Sqrt(gamma * pL_ghost / rL_ghost);
|
||
|
||
LaxFlux(rL_ghost, uL_ghost, pL_ghost, cL_ghost,
|
||
_rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), p[0], c[0],
|
||
out double fluxM_left, out double fluxP_left, out double fluxE_left);
|
||
|
||
double alphaLeft = Math.Max(Math.Abs(uL_ghost) + cL_ghost,
|
||
Math.Abs(_rhou[0] / Math.Max(_rho[0], 1e-12)) + c[0]);
|
||
ScalarFlux(rL_ghost, uL_ghost, cL_ghost, _YGhostL,
|
||
_rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), c[0], _Y[0],
|
||
alphaLeft, out double fluxY_left);
|
||
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profLeftFluxTicks += (t1 - t0);
|
||
t0 = t1;
|
||
}
|
||
|
||
// ---------- Phase 3: Interior loop (fluxes + cell updates) ----------
|
||
double fluxM_prev = fluxM_left;
|
||
double fluxP_prev = fluxP_left;
|
||
double fluxE_prev = fluxE_left;
|
||
double fluxY_prev = fluxY_left;
|
||
|
||
for (int i = 0; i < n - 1; i++)
|
||
{
|
||
int iL = i;
|
||
int iR = i + 1;
|
||
|
||
double rL = Math.Max(_rho[iL], 1e-12);
|
||
double uL = _rhou[iL] / rL;
|
||
double pL = p[iL];
|
||
double cL = c[iL];
|
||
double YL = _Y[iL];
|
||
|
||
double rR = Math.Max(_rho[iR], 1e-12);
|
||
double uR = _rhou[iR] / rR;
|
||
double pR = p[iR];
|
||
double cR = c[iR];
|
||
double YR = _Y[iR];
|
||
|
||
LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR,
|
||
out double fluxM_right, out double fluxP_right, out double fluxE_right);
|
||
|
||
double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
|
||
ScalarFlux(rL, uL, cL, YL, rR, uR, cR, YR, alpha, out double fluxY_right);
|
||
|
||
// Update cell i
|
||
double r = _rho[i];
|
||
double ru = _rhou[i];
|
||
double E = _E[i];
|
||
double Y = _Y[i];
|
||
|
||
double newR = r - dt_dx * (fluxM_right - fluxM_prev);
|
||
double newRu = ru - dt_dx * (fluxP_right - fluxP_prev);
|
||
double newE = E - dt_dx * (fluxE_right - fluxE_prev);
|
||
double oldRhoY = r * Y;
|
||
double newRhoY = oldRhoY - dt_dx * (fluxY_right - fluxY_prev);
|
||
|
||
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
|
||
newRu *= dampingFactor;
|
||
double relaxFactor = Math.Exp(-relaxRate * dt);
|
||
newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
|
||
|
||
newR = Math.Max(newR, 1e-12);
|
||
double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
|
||
double eMin = 100.0 / gm1 + kin;
|
||
newE = Math.Max(newE, eMin);
|
||
|
||
_rho[i] = newR;
|
||
_rhou[i] = newRu;
|
||
_E[i] = newE;
|
||
_Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0);
|
||
|
||
fluxM_prev = fluxM_right;
|
||
fluxP_prev = fluxP_right;
|
||
fluxE_prev = fluxE_right;
|
||
fluxY_prev = fluxY_right;
|
||
}
|
||
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profInteriorLoopTicks += (t1 - t0);
|
||
t0 = t1;
|
||
}
|
||
|
||
// ---------- Phase 4: Right face flux (cell n‑1 – ghostR) ----------
|
||
double rR_ghost = Math.Max(_rhoGhostR, 1e-12);
|
||
double pR_ghost = _pGhostR;
|
||
double uR_ghost = _uGhostR;
|
||
double cR_ghost = Math.Sqrt(gamma * pR_ghost / rR_ghost);
|
||
|
||
double rInt = _rho[n - 1];
|
||
double uInt = _rhou[n - 1] / Math.Max(rInt, 1e-12);
|
||
|
||
LaxFlux(rInt, uInt, p[n - 1], c[n - 1],
|
||
rR_ghost, uR_ghost, pR_ghost, cR_ghost,
|
||
out double fluxM_right_final, out double fluxP_right_final, out double fluxE_right_final);
|
||
|
||
double alphaRight = Math.Max(Math.Abs(uInt) + c[n - 1], Math.Abs(uR_ghost) + cR_ghost);
|
||
ScalarFlux(rInt, uInt, c[n - 1], _Y[n - 1],
|
||
rR_ghost, uR_ghost, cR_ghost, _YGhostR,
|
||
alphaRight, out double fluxY_right_final);
|
||
|
||
// Update last cell
|
||
{
|
||
int i = n - 1;
|
||
double r = _rho[i];
|
||
double ru = _rhou[i];
|
||
double E = _E[i];
|
||
double Y = _Y[i];
|
||
|
||
double newR = r - dt_dx * (fluxM_right_final - fluxM_prev);
|
||
double newRu = ru - dt_dx * (fluxP_right_final - fluxP_prev);
|
||
double newE = E - dt_dx * (fluxE_right_final - fluxE_prev);
|
||
double oldRhoY = r * Y;
|
||
double newRhoY = oldRhoY - dt_dx * (fluxY_right_final - fluxY_prev);
|
||
|
||
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
|
||
newRu *= dampingFactor;
|
||
double relaxFactor = Math.Exp(-relaxRate * dt);
|
||
newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
|
||
|
||
newR = Math.Max(newR, 1e-12);
|
||
double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
|
||
double eMin = 100.0 / gm1 + kin;
|
||
newE = Math.Max(newE, eMin);
|
||
|
||
_rho[i] = newR;
|
||
_rhou[i] = newRu;
|
||
_E[i] = newE;
|
||
_Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0);
|
||
}
|
||
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profRightFluxTicks += (t1 - t0);
|
||
t0 = t1;
|
||
}
|
||
|
||
// ---------- Phase 5: Update port states ----------
|
||
(double rhoA, double uA, double pA) = GetInteriorStateLeft();
|
||
PortA.Pressure = pA; PortA.Density = rhoA;
|
||
PortA.Temperature = pA / (rhoA * 287.0);
|
||
PortA.SpecificEnthalpy = gm1 / (gamma - 1.0) * pA / rhoA;
|
||
PortA.AirFraction = _Y[0];
|
||
|
||
(double rhoB, double uB, double pB) = GetInteriorStateRight();
|
||
PortB.Pressure = pB; PortB.Density = rhoB;
|
||
PortB.Temperature = pB / (rhoB * 287.0);
|
||
PortB.SpecificEnthalpy = gm1 / (gamma - 1.0) * pB / rhoB;
|
||
PortB.AirFraction = _Y[_n - 1];
|
||
|
||
if (EnableDetailedProfiling)
|
||
{
|
||
t1 = Stopwatch.GetTimestamp();
|
||
_profPortUpdateTicks += (t1 - t0);
|
||
}
|
||
}
|
||
|
||
private double PressureScalar(int i)
|
||
{
|
||
double rho = Math.Max(_rho[i], 1e-12);
|
||
return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
|
||
}
|
||
|
||
public void SetUniformState(double rho, double u, double p)
|
||
{
|
||
double e = p / ((_gamma - 1.0) * rho);
|
||
double E = rho * e + 0.5 * rho * u * u;
|
||
for (int i = 0; i < _n; i++)
|
||
{
|
||
_rho[i] = rho;
|
||
_rhou[i] = rho * u;
|
||
_E[i] = E;
|
||
_Y[i] = 1.0; // initially pure air
|
||
}
|
||
}
|
||
|
||
public void SetCellState(int i, double rho, double u, double p)
|
||
{
|
||
if (i < 0 || i >= _n) return;
|
||
double e = p / ((_gamma - 1.0) * rho);
|
||
double E = rho * e + 0.5 * rho * u * u;
|
||
_rho[i] = rho;
|
||
_rhou[i] = rho * u;
|
||
_E[i] = E;
|
||
_Y[i] = 1.0;
|
||
}
|
||
|
||
public void SetCellPressure(int i, double p)
|
||
{
|
||
if (i < 0 || i >= _n) return;
|
||
double rho = _rho[i];
|
||
double u = _rhou[i] / rho;
|
||
double e = p / ((_gamma - 1.0) * rho);
|
||
_E[i] = rho * e + 0.5 * rho * u * u;
|
||
}
|
||
|
||
// ---------- Public profiling interface ----------
|
||
public void ResetDetailCounters()
|
||
{
|
||
_profPrecomputeTicks = 0;
|
||
_profLeftFluxTicks = 0;
|
||
_profInteriorLoopTicks = 0;
|
||
_profRightFluxTicks = 0;
|
||
_profPortUpdateTicks = 0;
|
||
_profCallCount = 0;
|
||
}
|
||
|
||
public string GetDetailProfileReport()
|
||
{
|
||
if (!EnableDetailedProfiling)
|
||
return "Detailed profiling disabled.";
|
||
|
||
double freq = Stopwatch.Frequency;
|
||
long totalTicks = _profPrecomputeTicks + _profLeftFluxTicks +
|
||
_profInteriorLoopTicks + _profRightFluxTicks +
|
||
_profPortUpdateTicks;
|
||
|
||
if (totalTicks == 0) return "No profiling data.";
|
||
|
||
double totalSec = totalTicks / freq;
|
||
double avgCallSec = totalSec / _profCallCount;
|
||
double avgCallUs = avgCallSec * 1e6;
|
||
|
||
string report = $" Pipe detailed (over {_profCallCount} calls, total {totalSec * 1000:F2} ms):\n";
|
||
report += $" Avg per call: {avgCallUs:F2} µs\n";
|
||
report += $" Precompute p,c: {_profPrecomputeTicks * 100.0 / totalTicks:F1} % ({_profPrecomputeTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n";
|
||
report += $" Left face flux: {_profLeftFluxTicks * 100.0 / totalTicks:F1} % ({_profLeftFluxTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n";
|
||
report += $" Interior loop: {_profInteriorLoopTicks * 100.0 / totalTicks:F1} % ({_profInteriorLoopTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n";
|
||
report += $" Right face flux: {_profRightFluxTicks * 100.0 / totalTicks:F1} % ({_profRightFluxTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n";
|
||
report += $" Port update: {_profPortUpdateTicks * 100.0 / totalTicks:F1} % ({_profPortUpdateTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n";
|
||
return report;
|
||
}
|
||
}
|
||
} |