FluidSim/Scenarios/HelmholtzScenario.cs

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.2f * 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);
            pipeSystem.DampingMultiplier = 500f;

            // --- Boundary system ---
            boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1);

            float ComputeResistance(float decayTimeSeconds, float rho, float L_eff, float A)
            {
                // R = 2 * rho * L_eff / (A * decayTimeSeconds)
                return 2f * rho * L_eff / (A * MathF.Max(decayTimeSeconds, 1e-6f));
            }

            // Use steady orifice – the pipe already provides the inertia
            boundaries.AddOrificeWithInertance(
                cavityPort, pipeIndex: 0, isLeftEnd: true,
                areaIndex: cavityOrificeIdx,
                dischargeCoeff: 0.9f,
                effectiveLength: neckLength,   // physical length (or L_eff)
                lossCoefficient: 7000        // start with this, adjust for decay time
            );

            // 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 ---
            // Slightly higher sub‑step count to ensure stability of the resonant oscillation
            solver = new Solver { SubStepCount = 6, 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()
        {
            stepCount++;
            if (stepCount <= 8192) return 0f; // let buffer pre‑fill

            solver.Step();

            float flow = boundaries.GetOpenEndMassFlow(openEndIdx);
            float sample = soundProcessor.Process(flow);

            if (stepCount % 10000 == 0)
            {
                float cavityP = cavity.Pressure;
                float cavityT = cavity.Temperature;
                float cavityRho = cavity.Density;
                float cCavity = MathF.Sqrt(1.4f * cavityP / MathF.Max(cavityRho, 1e-12f));

                // Temperature in the middle of the neck
                int midCell = 10;
                float pMid = pipeSystem.GetCellPressure(midCell);
                float rhoMid = pipeSystem.GetCellDensity(midCell);
                float tMid = pMid / MathF.Max(rhoMid * 287f, 1e-12f);

                // Neck effective length (physical + end correction)
                float neckLen = 0.05f;            // physical
                float neckDia = 0.02f;
                float neckArea = MathF.PI * 0.25f * neckDia * neckDia;
                float endCorr = 0.85f * neckDia; // unflanged end
                float L_eff = neckLen + endCorr;

                // Theoretical Helmholtz frequency from current cavity sound speed
                float fHelmholtz = cCavity / (2f * MathF.PI) *
                                   MathF.Sqrt(neckArea / (cavity.Volume * L_eff));

                Console.WriteLine(
                    $"Step {stepCount}: cav P={cavityP / 1e5f:F4} bar, T={cavityT:F1} K, " +
                    $"pipeMid T={tMid:F1} K, est f={fHelmholtz:F1} Hz");
            }

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