commit a859740a95cde3b4dce302ed55777973fcd04101 Author: maxwes08 Date: Thu Oct 22 23:52:20 2026 +0200 Port FluidSim core and scenarios to C++20 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..11ed569 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +/build/ +/build-*/ +/out/ +/.vs/ +CMakeUserPresets.json +/verification/bin/ +/verification/obj/ diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..16af7a7 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,75 @@ +cmake_minimum_required(VERSION 3.28) +project(FluidSimCpp VERSION 0.1.0 LANGUAGES CXX) + +option(FLUIDSIM_BUILD_DESKTOP "Build the SFML windowed simulator with keyboard controls and audio" ON) + +add_library(fluidsim_core + src/engine.cpp + src/crankcase.cpp + src/two_stroke_cylinder.cpp + src/two_stroke_scenario.cpp + src/helmholtz_scenario.cpp + src/pulse_test_scenario.cpp + src/audio_processing.cpp + src/outdoor_exhaust_reverb.cpp + src/mechanical_sound_processor.cpp + src/scenario.cpp + src/solver.cpp + src/vehicle.cpp + src/simulation_ring_buffer.cpp +) +target_include_directories(fluidsim_core PUBLIC include) +target_compile_features(fluidsim_core PUBLIC cxx_std_20) + +add_executable(fluidsim_cpp src/main.cpp) +target_link_libraries(fluidsim_cpp PRIVATE fluidsim_core) + +add_executable(fluidsim_parity_trace tools/parity_trace.cpp) +target_link_libraries(fluidsim_parity_trace PRIVATE fluidsim_core) + +if(FLUIDSIM_BUILD_DESKTOP) + find_package(SFML 3.0.2 COMPONENTS Graphics Audio QUIET) + if(NOT SFML_FOUND) + include(FetchContent) + set(SFML_BUILD_AUDIO ON CACHE BOOL "" FORCE) + set(SFML_BUILD_GRAPHICS ON CACHE BOOL "" FORCE) + set(SFML_BUILD_WINDOW ON CACHE BOOL "" FORCE) + set(SFML_BUILD_NETWORK OFF CACHE BOOL "" FORCE) + set(SFML_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE) + set(SFML_BUILD_TEST_SUITE OFF CACHE BOOL "" FORCE) + set(BUILD_SHARED_LIBS OFF CACHE BOOL "" FORCE) + FetchContent_Declare(SFML + GIT_REPOSITORY https://github.com/SFML/SFML.git + GIT_TAG 3.0.2 + GIT_SHALLOW TRUE + EXCLUDE_FROM_ALL + SYSTEM + ) + FetchContent_MakeAvailable(SFML) + endif() + add_executable(fluidsim_desktop src/desktop_main.cpp src/audio_output_stream.cpp src/sound_engine.cpp) + target_link_libraries(fluidsim_desktop PRIVATE fluidsim_core SFML::Graphics SFML::Audio) + if(MSVC) + target_compile_options(fluidsim_desktop PRIVATE /W4 /permissive-) + else() + target_compile_options(fluidsim_desktop PRIVATE -Wall -Wextra -Wpedantic) + endif() + file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/../fonts" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") +endif() + +if(MSVC) + target_compile_options(fluidsim_core PRIVATE /W4 /permissive-) + target_compile_options(fluidsim_cpp PRIVATE /W4 /permissive-) + target_compile_options(fluidsim_parity_trace PRIVATE /W4 /permissive-) +else() + target_compile_options(fluidsim_core PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fluidsim_cpp PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fluidsim_parity_trace PRIVATE -Wall -Wextra -Wpedantic) +endif() + +include(CTest) +if(BUILD_TESTING) + add_executable(fluidsim_physics_tests tests/physics_tests.cpp) + target_link_libraries(fluidsim_physics_tests PRIVATE fluidsim_core) + add_test(NAME fluidsim_physics_tests COMMAND fluidsim_physics_tests) +endif() diff --git a/README.md b/README.md new file mode 100644 index 0000000..dd11872 --- /dev/null +++ b/README.md @@ -0,0 +1,83 @@ +# FluidSimCpp + +Native C++20 port of FluidSim. The C# application remains the reference while +the component models and behavior are being ported and checked against it. + +## Build (Windows / Visual Studio) + +Open a Visual Studio Developer PowerShell (or initialize the MSVC environment +with `VsDevCmd.bat`), then run from this directory: + +```powershell +cmake -S . -B build -G "NMake Makefiles" -DCMAKE_BUILD_TYPE=Release +cmake --build build +ctest --test-dir build --output-on-failure +``` + +The first configure fetches SFML 3.0.2 with CMake FetchContent for the desktop +window and audio stream. `-DFLUIDSIM_BUILD_DESKTOP=OFF` builds only the headless +runner and physics tests. + +Run the interactive simulator with: + +```powershell +.\build\fluidsim_desktop.exe +``` + +Controls: hold `W` for throttle, `C` for clutch, `A` / `D` to shift, `Z` for +neutral, up/down to adjust throttle target, left/right to adjust engine load, +`M` toggles mechanical engine noise, mouse wheel changes simulation playback +speed, and Space toggles real time. + +The command-line runner is useful for repeatable checks: + +```powershell +.\build\fluidsim_cpp.exe --seconds 1 --throttle 0.35 +.\build\fluidsim_cpp.exe --scenario two-stroke --seconds 1 --throttle 0.3 +.\build\fluidsim_cpp.exe --scenario helmholtz --seconds 1 +.\build\fluidsim_cpp.exe --scenario pulse --seconds 1 +``` + +To compare deterministic checkpoints from the active C# twin scenario against +the C++ port (energy variation is disabled for repeatability), build the C++ +trace runner and execute the PowerShell comparison script: + +```powershell +cmake --build build-msvc --target fluidsim_parity_trace +.\verification\Compare-ScenarioTrace.ps1 +``` + +## Current port status + +The C++ core now includes the shared `EngineCylinder` model with four-stroke and +two-stroke cycle implementations, crankshaft, crankcase, atmosphere and +configurable 0-D volume components (including thermal relaxation and moving +volume work), profiled compressible 1-D pipes and boundaries, adaptive solver +stepping, the impulse-clutch motorcycle/road-load model, and the active 600 cc +twin topology. Audio includes the flow-derivative source, spatialized +intake/exhaust radiation, outdoor feedback-delay reverb, optional mechanical +noise, a lock-free simulation audio buffer, a `SoundEngine` lifecycle wrapper, +and an SFML desktop frontend. The shared 50-RPM-bin dyno history and graph are +also ported. A crankcase-scavenged 125 cc two-stroke scenario is available in +the core. + +This is still an in-progress port, not a claimed complete replacement. All four +current C# scenarios (600 cc twin, 125 cc two-stroke, Helmholtz resonator, and +pressure pulse) have headless C++ counterparts and can be selected by the CLI. +Remaining parity work includes completing the C# drawing/diagnostic controls in +the desktop UI and resolving numerical differences reported by the comparison +harness. The latest MSVC build completed, CTest passed, and the one-second +headless twin run completed with finite output. The desktop executable opened, +but its client area rendered blank in the captured run, so desktop rendering is +not yet verified. + +The one-second twin trace remains substantially outside its 1% tolerance. The +latest comparison's largest relative difference is 40.6745 for cylinder 1 +pressure at checkpoint 33075 (69.8 kPa in C# versus 2.91 MPa in C++). A dense +intake trace matches closely at first, then the C# plenum pressure begins +rising at checkpoint 63 while C++ remains near ambient. This divergence needs +to be resolved before claiming numerical parity. CTest covers cylinder +geometry and valve timing, dyno-bin interpolation, configurable 0-D volume +thermodynamics, crankcase and atmosphere state, profiled-pipe geometry, +inertive and closed boundaries, gearing and clutch momentum transfer, finite +runs for every scenario, and finite audio/reverb output. diff --git a/include/FluidSim/Atmosphere.hpp b/include/FluidSim/Atmosphere.hpp new file mode 100644 index 0000000..c9bc006 --- /dev/null +++ b/include/FluidSim/Atmosphere.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include "FluidSim/Engine.hpp" + +namespace fluidsim::port { + +// Fixed thermodynamic boundary matching the C# Atmosphere component. +class Atmosphere final : public Component { +public: + float pressure{101325.0f}; + float temperature{300.0f}; + float gas_constant_value{gas_constant}; + float gamma_value{gamma}; + Port port{}; + + Atmosphere() { update_port(); } + + float density() const { return pressure / (gas_constant_value * temperature); } + float specific_enthalpy() const { + return gamma_value / (gamma_value - 1.0f) * pressure / density(); + } + std::vector ports() override { return {&port}; } + void update_state(float /*dt*/) override { update_port(); } + +private: + void update_port() { + port.pressure = pressure; + port.density = density(); + port.temperature = temperature; + port.enthalpy = specific_enthalpy(); + port.air_fraction = 1.0f; + } +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/AudioOutputStream.hpp b/include/FluidSim/AudioOutputStream.hpp new file mode 100644 index 0000000..f3d4e41 --- /dev/null +++ b/include/FluidSim/AudioOutputStream.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include "FluidSim/SimulationRingBuffer.hpp" +#include +#include +#include +#include + +namespace fluidsim::port { + +class AudioOutputStream final : public sf::SoundStream { +public: + explicit AudioOutputStream(SimulationRingBuffer& ring, unsigned int sample_rate = 44100); + void set_speed(double speed) { speed_.store(speed, std::memory_order_relaxed); } + double speed() const { return speed_.load(std::memory_order_relaxed); } + +private: + static constexpr std::size_t block_size = 512; + SimulationRingBuffer& ring_; + std::array temp_left_{}, temp_right_{}; + std::array output_{}; + std::atomic speed_{1.0}; + float previous_input_left_{}, previous_output_left_{}; + float previous_input_right_{}, previous_output_right_{}; + + bool onGetData(Chunk& data) override; + void onSeek(sf::Time time_offset) override; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/AudioProcessing.hpp b/include/FluidSim/AudioProcessing.hpp new file mode 100644 index 0000000..a05a0e7 --- /dev/null +++ b/include/FluidSim/AudioProcessing.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include + +namespace fluidsim::port { + +struct AudioFrame { float left{}, right{}; }; +struct Vec3 { float x{}, y{}, z{}; }; + +// Flow-derivative sound source used by the C# two-stroke and pulse scenarios. +class SoundProcessor { +public: + float gain{1.0f}; + explicit SoundProcessor(int sample_rate = 44100, float listener_distance = 1.0f); + float process(float mass_flow_out); +private: + float dt_{}, scale_factor_{}, lowpass_alpha_{}, derivative_alpha_{}; + float flow_lowpass_{}, previous_flow_{}, smooth_derivative_{}; +}; + +class RadiationSource { +public: + float flow_tau{0.005f}, smooth_tau{0.020f}, gain{15.0f}; + Vec3 pipe_end{}, pipe_axis{1.0f, 0.0f, 0.0f}; + Vec3 listener{0.0f, 0.0f, 1.5f}, listener_right{1.0f, 0.0f, 0.0f}; + float ear_spacing{0.18f}, head_shadow_cutoff{4000.0f}, speed_of_sound{343.0f}; + float directivity_balance{0.30f}, air_absorption_exponent{0.002f}, min_air_cutoff{1500.0f}; + + explicit RadiationSource(int sample_rate = 44100); + AudioFrame process(float mass_flow); + float last_mono_pressure() const { return last_mono_pressure_; } + +private: + class DelayLine; + class OnePole; + int sample_rate_{}; + float dt_{}, flow_lp_{}, previous_flow_lp_{}, smooth_derivative_{}; + float last_mono_pressure_{}; + std::vector left_delay_, right_delay_; + std::size_t write_position_{}, left_offset_{}, right_offset_{}; + float air_filter_{}, shadow_left_{}, shadow_right_{}; + + static float dot(Vec3 a, Vec3 b); + static float length(Vec3 value); + static Vec3 subtract(Vec3 a, Vec3 b); + float one_pole(float input, float cutoff, float& state) const; + float delayed(float input, std::vector& buffer, std::size_t offset); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Component.hpp b/include/FluidSim/Component.hpp new file mode 100644 index 0000000..476cbd9 --- /dev/null +++ b/include/FluidSim/Component.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include + +namespace fluidsim::port { + +struct Port { + float mass_flow{}; // kg/s, positive into the owning component + float enthalpy{}; // J/kg + float pressure{101325.0f}; + float density{1.1768f}; + float temperature{300.0f}; + float air_fraction{1.0f}; + float* owner_air_mass{}; + float* owner_exhaust_mass{}; +}; + +// Small equivalent of the C# IComponent lifecycle. Components update their +// conserved state after the coupled 1-D pipe solve has accumulated port flow. +class Component { +public: + virtual ~Component() = default; + virtual std::vector ports() = 0; + virtual void update_state(float dt) = 0; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Constants.hpp b/include/FluidSim/Constants.hpp new file mode 100644 index 0000000..6e24ede --- /dev/null +++ b/include/FluidSim/Constants.hpp @@ -0,0 +1,16 @@ +#pragma once + +namespace fluidsim::port { + +inline constexpr float pi = 3.14159265358979323846f; +inline constexpr float gamma = 1.4f; +// PipeSystem.cs stores these as independent float constants rather than +// deriving them from Gamma at runtime; preserve that rounding in the port. +inline constexpr float pipe_gm1 = 0.4f; +inline constexpr float pipe_gm1_inverse = 1.0f / pipe_gm1; +inline constexpr float gas_constant = 287.0f; +inline constexpr float ambient_pressure = 101325.0f; +inline constexpr float ambient_temperature = 300.0f; +inline constexpr float ambient_density = ambient_pressure / (gas_constant * ambient_temperature); + +} // namespace fluidsim::port diff --git a/include/FluidSim/Crankcase.hpp b/include/FluidSim/Crankcase.hpp new file mode 100644 index 0000000..10d9589 --- /dev/null +++ b/include/FluidSim/Crankcase.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include "FluidSim/Component.hpp" +#include "FluidSim/Engine.hpp" +#include + +namespace fluidsim::port { + +class Crankcase final : public Component { +public: + Crankcase(Crankshaft& crankshaft, float crank_radius, float connecting_rod_length, + float bore, float clearance_volume, float initial_pressure = 101325.0f, + float initial_temperature = 300.0f); + + Port intake_port{}, transfer_port{}; + float pressure() const { return pressure_; } + float temperature() const { return temperature_; } + float density() const { return mass_ / (volume_ > 1e-12f ? volume_ : 1e-12f); } + float volume() const { return volume_; } + float mass() const { return mass_; } + float air_fraction() const { return air_fraction_; } + + std::vector ports() override { return {&intake_port, &transfer_port}; } + void pre_step(float dt); + void update_state(float dt) override; + +private: + Crankshaft& crankshaft_; + float crank_radius_{}, connecting_rod_length_{}, piston_area_{}; + float clearance_volume_{}, obliquity_{}; + float mass_{}, energy_{}, air_fraction_{1.0f}; + float pressure_{101325.0f}, temperature_{300.0f}, volume_{}, previous_volume_{}; + + void sync_ports(); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/DynoCurve.hpp b/include/FluidSim/DynoCurve.hpp new file mode 100644 index 0000000..dc91288 --- /dev/null +++ b/include/FluidSim/DynoCurve.hpp @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include +#include + +namespace fluidsim::port { + +struct DynoPoint { + float power_kw{}; + float torque_nm{}; +}; + +// Shared rolling peak curve from the C# Scenario base class. Each bin stores +// the greatest observed power and torque independently, with missing bins +// interpolated when the engine crosses an RPM range too quickly to sample it. +class DynoCurve { +public: + static constexpr float rpm_bin_size = 50.0f; + + void reset() { + bins_.clear(); + last_bin_ = invalid_bin; + } + + void update(float rpm, float power_kw, float torque_nm) { + if (!std::isfinite(rpm) || rpm <= 0.0f || !std::isfinite(power_kw) || !std::isfinite(torque_nm)) + return; + const auto bin = static_cast(rpm / rpm_bin_size); + if (bin >= bins_.max_size()) return; + while (bins_.size() <= bin) bins_.push_back({}); + + if (last_bin_ != invalid_bin && bin > last_bin_ + 1) { + const float last_power = bins_[last_bin_].power_kw > 0.0f ? bins_[last_bin_].power_kw : 0.0f; + const float last_torque = bins_[last_bin_].torque_nm > 0.0f ? bins_[last_bin_].torque_nm : 0.0f; + for (std::size_t missing = last_bin_ + 1; missing < bin; ++missing) { + const float t = static_cast(missing - last_bin_) / + static_cast(bin - last_bin_); + const float interpolated_power = last_power + (power_kw - last_power) * t; + const float interpolated_torque = last_torque + (torque_nm - last_torque) * t; + auto& point = bins_[missing]; + if (interpolated_power > point.power_kw || point.power_kw <= 0.0f) + point.power_kw = interpolated_power; + if (interpolated_torque > point.torque_nm || point.torque_nm <= 0.0f) + point.torque_nm = interpolated_torque; + } + } + + auto& current = bins_[bin]; + if (power_kw > current.power_kw || current.power_kw <= 0.0f) current.power_kw = power_kw; + if (torque_nm > current.torque_nm || current.torque_nm <= 0.0f) current.torque_nm = torque_nm; + last_bin_ = bin; + } + + const std::vector& bins() const { return bins_; } + +private: + static constexpr std::size_t invalid_bin = static_cast(-1); + std::vector bins_; + std::size_t last_bin_{invalid_bin}; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Engine.hpp b/include/FluidSim/Engine.hpp new file mode 100644 index 0000000..ce3ec1a --- /dev/null +++ b/include/FluidSim/Engine.hpp @@ -0,0 +1,75 @@ +#pragma once + +#include "FluidSim/Component.hpp" +#include "FluidSim/Constants.hpp" +#include +#include +#include + +namespace fluidsim::port { + +class Crankshaft { +public: + float omega{400.0f * 2.0f * pi / 60.0f}; + float angle{}, previous_angle{}, inertia{0.2f}; + float friction_constant{}, friction_viscous{}; + float load_torque{}, drivetrain_torque{}; + float average_power{}, average_torque{}; + + void add_torque(float torque); + void step(float dt); +private: + float external_torque_{}; + std::array power_buffer_{}, torque_buffer_{}; + std::size_t power_index_{}, power_count_{}, torque_index_{}, torque_count_{}; + float power_sum_{}, torque_sum_{}; +}; + +class EngineCylinder : public Component { +public: + float bore{0.075f}, stroke{0.068f}, rod{0.112f}, compression{11.5f}; + float intake_valve_diameter{0.030f}, exhaust_valve_diameter{0.028f}; + float intake_valve_lift{0.005f}, exhaust_valve_lift{0.005f}; + float phase{}; + float spark_advance{20.0f}, wiebe_a{5.0f}, wiebe_m{2.0f}, burn_duration{60.0f}, burn_start{5.0f}; + float afr{14.7f}, fuel_lhv{44.0e6f}, energy_variation{0.05f}, misfire_probability{}; + float wall_area{0.02f}, heat_transfer{100.0f}, ambient_temperature{300.0f}; + bool rev_limit_enabled{}; + float rev_limit_rpm{7500.0f}, rev_limit_cut_time{0.05f}, rev_limit_hysteresis{200.0f}; + float air_mass{}, exhaust_mass{}, energy{}, volume{}, pressure{101325.0f}; + float fuel_mass{}, burn_fraction{}; + bool combustion_active{}, fuel_injected{}, rev_limit_active{}; + Port intake{}, exhaust{}; + + std::vector ports() override { return {&intake, &exhaust}; } + float swept_volume() const; + float clearance_volume() const; + float total_mass() const; + float temperature() const; + float air_fraction() const; + float compute_volume(float theta) const; + void initialize(const Crankshaft& crank); + virtual float valve_area(float cycle_degrees, bool inlet) const = 0; + virtual void pre_step(Crankshaft& crank, float dt, float throttle) = 0; + void update_state(float dt) override; + +protected: + float energy_factor_{1.0f}; + float rev_limit_timer_{}; + std::mt19937 random_{0xF17D51u}; + std::uniform_real_distribution unit_{0.0f, 1.0f}; + + float random_unit(); + static float wrap_degrees(float value, float cycle_degrees = 720.0f); + static bool crossed_window_close(float previous, float current, float open, float close); + void update_rev_limiter(float rpm, float dt); + void sync_ports(); +}; + +class FourStrokeCylinder final : public EngineCylinder { +public: + float intake_open{340.0f}, intake_close{590.0f}, exhaust_open{130.0f}, exhaust_close{390.0f}; + float valve_area(float cycle_degrees, bool inlet) const override; + void pre_step(Crankshaft& crank, float dt, float throttle) override; +}; +} // namespace fluidsim::port diff --git a/include/FluidSim/Fluid.hpp b/include/FluidSim/Fluid.hpp new file mode 100644 index 0000000..3519766 --- /dev/null +++ b/include/FluidSim/Fluid.hpp @@ -0,0 +1,387 @@ +#pragma once + +#include "FluidSim/Engine.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace fluidsim::port { + +class Volume0D : public Component { +public: + explicit Volume0D(float volume_m3, float pressure = 101325.0f, float temperature = 300.0f, + float gas_constant_value = 287.0f, float gamma_value = 1.4f) + : gamma(gamma_value), gas_constant(gas_constant_value), volume_(volume_m3), + air_mass_(pressure * volume_m3 / (gas_constant_value * temperature)), + energy_(pressure * volume_m3 / (gamma_value - 1.0f)) {} + + float gamma{1.4f}; + float gas_constant{287.0f}; + float ambient_pressure{101325.0f}; + float ambient_temperature{300.0f}; + float energy_relaxation_rate{10.0f}; + float volume_change_rate{}; + + Port& create_port() { + ports_.emplace_back(); + auto& p = ports_.back(); + p.owner_air_mass = &air_mass_; + p.owner_exhaust_mass = &exhaust_mass_; + sync_port(p); + return p; + } + float pressure() const { return (gamma - 1.0f) * energy_ / std::max(volume_, 1e-12f); } + float density() const { return mass() / std::max(volume_, 1e-12f); } + float temperature() const { return pressure() / std::max(density() * gas_constant, 1e-9f); } + float volume() const { return volume_; } + float mass() const { return air_mass_ + exhaust_mass_; } + float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); } + void set_pressure(float value, float temp = -1.0f) { + const float target_t = temp > 0.0f ? temp : temperature(); + const float total = value * volume_ / (gas_constant * target_t); + const float fraction = air_fraction(); + air_mass_ = total * fraction; + exhaust_mass_ = total * (1.0f - fraction); + energy_ = value * volume_ / (gamma - 1.0f); + sync_ports(); + } + void update(float dt) { + if (dt <= 0.0f) return; + float d_air = 0.0f, d_exhaust = 0.0f, d_energy = 0.0f; + const float fraction = air_fraction(); + for (const Port& p : ports_) { + const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction; + d_air += p.mass_flow * flow_fraction; + d_exhaust += p.mass_flow * (1.0f - flow_fraction); + d_energy += p.mass_flow * p.enthalpy; + } + air_mass_ += d_air * dt; + exhaust_mass_ += d_exhaust * dt; + energy_ += (d_energy - pressure() * volume_change_rate) * dt; + if (energy_relaxation_rate > 0.0f && mass() > 1e-12f) { + const float target = mass() * gas_constant * ambient_temperature / (gamma - 1.0f); + energy_ = target + (energy_ - target) * std::exp(-energy_relaxation_rate * dt); + } + if (mass() < 1e-9f || energy_ < 0.0f) { + air_mass_ = 1e-9f; + exhaust_mass_ = 0.0f; + energy_ = ambient_pressure * volume_ / (gamma - 1.0f); + } + air_mass_ = std::max(air_mass_, 0.0f); + exhaust_mass_ = std::max(exhaust_mass_, 0.0f); + sync_ports(); + } + void update_state(float dt) override { update(dt); } + std::vector ports() override { + std::vector result; + result.reserve(ports_.size()); + for (Port& port : ports_) result.push_back(&port); + return result; + } + const std::deque& ports() const { return ports_; } +private: + float volume_, air_mass_{}, exhaust_mass_{}, energy_; + std::deque ports_; + void sync_port(Port& p) const { + p.pressure = pressure(); + p.density = density(); + p.temperature = temperature(); + p.enthalpy = gamma / (gamma - 1.0f) * p.pressure / std::max(p.density, 1e-12f); + p.air_fraction = air_fraction(); + } + void sync_ports() { for (Port& p : ports_) sync_port(p); } +}; + +struct GasCell { + float rho{1.225f}, momentum{}; + float energy{1.225f * (101325.0f / (pipe_gm1 * 1.225f))}; + float air{1.0f}; +}; +struct GasState { float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f}; }; + +class Pipe1D { +public: + float area{}, dx{}; + float damping_multiplier{1.0f}, energy_relaxation_rate{0.5f}, ambient_pressure{101325.0f}; + std::vector cell_areas, cell_lengths; + std::vector cells; + GasState left_ghost{}, right_ghost{}; + Pipe1D(float area_m2, float length_m, int cell_count) + : area(area_m2), dx(length_m / static_cast(cell_count)), + cell_areas(static_cast(cell_count), area_m2), + cell_lengths(static_cast(cell_count), dx), + cells(static_cast(cell_count)), next_cells_(static_cast(cell_count)) {} + Pipe1D(std::vector areas_m2, std::vector lengths_m) + : area(areas_m2.empty() ? 0.0f : areas_m2.front()), + dx(lengths_m.empty() ? 0.0f : lengths_m.front()), + cell_areas(std::move(areas_m2)), cell_lengths(std::move(lengths_m)), + cells(cell_areas.size()), next_cells_(cell_areas.size()) { + if (cell_areas.empty() || cell_lengths.size() != cell_areas.size()) + throw std::invalid_argument("pipe area and length profiles must be non-empty and have equal length"); + for (std::size_t i = 0; i < cell_areas.size(); ++i) { + if (cell_areas[i] <= 0.0f || cell_lengths[i] <= 0.0f) + throw std::invalid_argument("pipe cell area and length must be positive"); + } + } + GasState state(int idx) const { + const auto& c = cells.at(static_cast(idx)); + const float rho = c.rho; + const float u = c.momentum / std::max(rho, 1e-12f); + const float p = pipe_gm1 * (c.energy - 0.5f * c.momentum * u); + return {rho, u, p, std::clamp(c.air, 0.0f, 1.0f)}; + } + float cell_mass(int idx) const { + const auto cell = static_cast(idx); + return cells.at(cell).rho * cell_areas.at(cell) * cell_lengths.at(cell); + } + float cell_area(int idx) const { return cell_areas.at(static_cast(idx)); } + float cell_length(int idx) const { return cell_lengths.at(static_cast(idx)); } + float minimum_cell_length() const { + return *std::min_element(cell_lengths.begin(), cell_lengths.end()); + } + int required_substeps(float dt, float target_cfl = 0.8f) const { + float maximum_wave_speed = 0.0f; + for (int i = 0; i < static_cast(cells.size()); ++i) { + const GasState s = state(i); + maximum_wave_speed = std::max(maximum_wave_speed, + std::abs(s.velocity) + std::sqrt(gamma * s.pressure / std::max(s.rho, 1e-12f))); + } + return std::max(1, static_cast(std::ceil(dt * maximum_wave_speed / + (target_cfl * std::max(minimum_cell_length(), 1e-12f))))); + } + void advance(float dt) { + const int n = static_cast(cells.size()); + if (n == 0 || dt <= 0.0f) return; + const float ambient_energy = ambient_pressure * pipe_gm1_inverse; + for (int i = 0; i < n; ++i) { + const GasState l = i == 0 ? left_ghost : flux_state(i - 1), r = flux_state(i); + const GasState a = flux_state(i), b = i == n - 1 ? right_ghost : flux_state(i + 1); + const Flux fl = flux(l, r), fr = flux(a, b); + const float factor = dt / cell_lengths[static_cast(i)]; + float rho = cells[i].rho - factor * (fr.m - fl.m); + float mom = cells[i].momentum - factor * (fr.p - fl.p); + float en = cells[i].energy - factor * (fr.e - fl.e); + const float rho_air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y); + if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || !std::isfinite(rho_air)) { + next_cells_[i] = cells[i]; + continue; + } + rho = std::max(rho, 1e-12f); + const float damping = std::exp(-0.1f * damping_multiplier * dt / + std::max(cells[i].rho, 1e-12f)); + mom *= damping; + if (energy_relaxation_rate > 0.0f) { + const float relax = std::exp(-energy_relaxation_rate * dt); + en = ambient_energy + (en - ambient_energy) * relax; + } + en = std::max(en, 100.0f * pipe_gm1_inverse + 0.5f * mom * mom / rho); + next_cells_[i] = {rho, mom, en, std::clamp(rho_air / rho, 0.0f, 1.0f)}; + } + cells.swap(next_cells_); + } +private: + struct Flux { float m{}, p{}, e{}, y{}; }; + std::vector next_cells_; + GasState flux_state(int idx) const { + const auto& c = cells[static_cast(idx)]; + const float rho = c.rho; + const float inverse_rho = 1.0f / std::max(rho, 1e-12f); + const float u = c.momentum * inverse_rho; + const float p = pipe_gm1 * (c.energy - 0.5f * c.momentum * u); + return {rho, u, p, std::clamp(c.air, 0.0f, 1.0f)}; + } + static Flux flux(const GasState& a, const GasState& b) { + const float ua = a.velocity, ub = b.velocity; + const float inv_rho_a = 1.0f / std::max(a.rho, 1e-12f); + const float inv_rho_b = 1.0f / std::max(b.rho, 1e-12f); + const float ea = a.pressure * pipe_gm1_inverse * inv_rho_a + 0.5f * ua * ua; + const float eb = b.pressure * pipe_gm1_inverse * inv_rho_b + 0.5f * ub * ub; + const float c = std::max(std::abs(ua) + std::sqrt(gamma * a.pressure * inv_rho_a), + std::abs(ub) + std::sqrt(gamma * b.pressure * inv_rho_b)); + return { + 0.5f * (a.rho * ua + b.rho * ub) - 0.5f * c * (b.rho - a.rho), + 0.5f * (a.rho * ua * ua + a.pressure + b.rho * ub * ub + b.pressure) - 0.5f * c * (b.rho * ub - a.rho * ua), + 0.5f * ((a.rho * ea + a.pressure) * ua + (b.rho * eb + b.pressure) * ub) - 0.5f * c * (b.rho * eb - a.rho * ea), + 0.5f * (a.rho * ua * a.air + b.rho * ub * b.air) - 0.5f * c * (b.rho * b.air - a.rho * a.air) + }; + } +}; + +class BoundarySystem { +public: + explicit BoundarySystem(std::vector& pipes) : pipes_(pipes) {} + void add_orifice(Port& port, int pipe, bool left, float& area, float cd = 1.0f) { + Orifice desc; + desc.port = &port; + desc.pipe = pipe; + desc.left = left; + desc.area = &area; + desc.cd = cd; + orifices_.push_back(desc); + } + void add_orifice_with_inertance(Port& port, int pipe, bool left, float& area, + float cd, float effective_length, float loss_coefficient = 0.0f) { + add_orifice(port, pipe, left, area, cd); + auto& desc = orifices_.back(); + desc.use_inertance = true; + desc.effective_length = effective_length; + desc.loss_coefficient = loss_coefficient; + } + void add_closed_end(int pipe, bool left) { + Orifice desc; + desc.pipe = pipe; + desc.left = left; + desc.closed_end = true; + orifices_.push_back(desc); + } + void add_open_end(int pipe, bool left, float pressure, float area) { opens_.push_back({pipe, left, pressure, area}); } + float open_end_mass_flow(std::size_t index) const { return opens_.at(index).last_mass_flow; } + float open_end_pressure(std::size_t index) const { return opens_.at(index).last_face_pressure; } + void resolve(float dt) { + for (auto& o : orifices_) resolve_orifice(o, dt); + for (auto& o : opens_) resolve_open_end(o, dt); + } +private: + struct Orifice { + Port* port{}; + int pipe{}; + bool left{}; + float* area{}; + float cd{1.0f}; + bool closed_end{}, use_inertance{}; + float effective_length{}, loss_coefficient{}, current_mass_flow{}; + }; + struct Open { int pipe; bool left; float pressure, area; float last_mass_flow{}, last_face_pressure{}; }; + struct OrificeSolution { float mass_flow{}, rho{}, velocity{}, pressure{}; }; + std::vector& pipes_; + std::vector orifices_; + std::vector opens_; + + static OrificeSolution isentropic(float pup, float rho_up, float t_up, float pdown, float area, float cd) { + OrificeSolution out{0.0f, rho_up, 0.0f, pup}; + if (area <= 0.0f || pup <= 0.0f || rho_up <= 0.0f || t_up <= 0.0f) return out; + float ratio = std::clamp(pdown / pup, 1e-6f, 1.0f); + const float critical = std::pow(2.0f / (gamma + 1.0f), gamma / (gamma - 1.0f)); + ratio = std::max(ratio, critical); + const float mach = std::sqrt(std::max(0.0f, 2.0f / (gamma - 1.0f) * + (std::pow(ratio, -(gamma - 1.0f) / gamma) - 1.0f))); + const float a_up = std::sqrt(gamma * gas_constant * t_up); + out.velocity = mach * a_up; + out.rho = rho_up * std::pow(ratio, 1.0f / gamma); + out.pressure = pup * ratio; + out.mass_flow = out.rho * out.velocity * area * cd; + return out; + } + static float available_mass(const Port& port) { + return (port.owner_air_mass ? *port.owner_air_mass : 0.0f) + + (port.owner_exhaust_mass ? *port.owner_exhaust_mass : 0.0f); + } + void resolve_orifice(Orifice& o, float dt) { + Pipe1D& pipe = pipes_.at(static_cast(o.pipe)); + const int idx = o.left ? 0 : static_cast(pipe.cells.size()) - 1; + const GasState interior = pipe.state(idx); + const float area = o.area ? std::max(*o.area, 0.0f) : 0.0f; + if (o.closed_end || !o.port || area <= 1e-12f) { + GasState wall = interior; + wall.velocity = -interior.velocity; + if (o.left) pipe.left_ghost = wall; else pipe.right_ghost = wall; + if (o.port) o.port->mass_flow = 0.0f; + return; + } + + const float dp = o.port->pressure - interior.pressure; + OrificeSolution solution{0.0f, interior.rho, 0.0f, interior.pressure}; + // Preserve the C# tie-break: at equal pressure it treats the volume as + // upstream, retaining the volume-side density in the ghost state. That + // matters when equal-pressure volume and pipe states have different + // temperatures/densities, as on the initialized intake plenum. + if (o.port->pressure >= interior.pressure) { + solution = isentropic(o.port->pressure, o.port->density, o.port->temperature, interior.pressure, area, o.cd); + } else { + solution = isentropic(interior.pressure, interior.rho, + interior.pressure / std::max(interior.rho * gas_constant, 1e-12f), o.port->pressure, area, o.cd); + solution.mass_flow = -solution.mass_flow; + } + float mdot = solution.mass_flow; // positive means port component -> pipe + const float adjacent_mass = pipe.cell_mass(idx); + if (o.use_inertance) { + const float rho_up = o.current_mass_flow >= 0.0f ? o.port->density : interior.rho; + const float inertance = rho_up * o.effective_length / area; + const float derivative = (dp - o.loss_coefficient * o.current_mass_flow) / + std::max(inertance, 1e-12f); + mdot = o.current_mass_flow + derivative * dt; + } + if (mdot > 0.0f && (o.port->owner_air_mass || o.port->owner_exhaust_mass)) + mdot = std::min(mdot, available_mass(*o.port) / std::max(dt, 1e-12f)); + if (mdot < 0.0f) mdot = std::max(mdot, -adjacent_mass / std::max(dt, 1e-12f)); + + if (o.use_inertance) { + o.current_mass_flow = mdot; + const bool leaves_port = mdot >= 0.0f; + const float rho_up = leaves_port ? o.port->density : interior.rho; + const float temperature_up = leaves_port ? o.port->temperature : + interior.pressure / std::max(interior.rho * gas_constant, 1e-12f); + const float max_velocity = 0.9f * std::sqrt(gamma * gas_constant * temperature_up); + const float max_flow = rho_up * max_velocity * area; + mdot = std::clamp(mdot, -max_flow, max_flow); + o.current_mass_flow = mdot; + if (mdot >= 0.0f) { + solution = isentropic(o.port->pressure, o.port->density, o.port->temperature, + interior.pressure, area, o.cd); + } else { + solution = isentropic(interior.pressure, interior.rho, temperature_up, + o.port->pressure, area, o.cd); + solution.mass_flow = -solution.mass_flow; + } + solution.mass_flow = mdot; + solution.rho = rho_up; + } + + const bool from_port = mdot >= 0.0f; + const float rho_face = std::max(solution.rho, 1e-9f); + const float u = std::abs(mdot) / std::max(rho_face * area, 1e-12f); + GasState ghost{rho_face, 0.0f, std::max(solution.pressure, 1.0f), from_port ? o.port->air_fraction : interior.air}; + ghost.velocity = ((from_port == o.left) ? u : -u); + if (o.left) pipe.left_ghost = ghost; else pipe.right_ghost = ghost; + o.port->mass_flow = -mdot; + if (!from_port) o.port->air_fraction = interior.air; + o.port->enthalpy = from_port + ? o.port->enthalpy + : gamma / (gamma - 1.0f) * interior.pressure / std::max(interior.rho, 1e-12f); + (void)dp; + } + void resolve_open_end(Open& o, float dt) { + Pipe1D& pipe = pipes_.at(static_cast(o.pipe)); + const int idx = o.left ? 0 : static_cast(pipe.cells.size()) - 1; + const GasState interior = pipe.state(idx); + const float c = std::sqrt(gamma * interior.pressure / std::max(interior.rho, 1e-12f)); + const float gm1 = gamma - 1.0f; + const float j_plus = interior.velocity + 2.0f * c / gm1; + const float j_minus = interior.velocity - 2.0f * c / gm1; + const float entropy = interior.pressure / std::pow(std::max(interior.rho, 1e-12f), gamma); + const float rho_iso = std::pow(std::max(o.pressure / entropy, 1e-12f), 1.0f / gamma); + const float c_iso = std::sqrt(gamma * o.pressure / std::max(rho_iso, 1e-12f)); + const float u_iso = o.left ? j_minus + 2.0f * c_iso / gm1 : j_plus - 2.0f * c_iso / gm1; + const bool supersonic = o.left ? (interior.velocity <= -c || u_iso <= -c_iso) + : (interior.velocity >= c || u_iso >= c_iso); + GasState ghost = supersonic ? interior : GasState{rho_iso, u_iso, o.pressure, + ((o.left ? u_iso >= 0.0f : u_iso <= 0.0f) ? 1.0f : interior.air)}; + + const float cell_mass = pipe.cell_mass(idx); + float mdot_out = ghost.rho * ghost.velocity * o.area * (o.left ? -1.0f : 1.0f); + if (mdot_out > 0.0f) mdot_out = std::min(mdot_out, cell_mass / std::max(dt, 1e-12f)); + else mdot_out = std::max(mdot_out, -10.0f * cell_mass / std::max(dt, 1e-12f)); + const float speed = std::abs(mdot_out) / std::max(ghost.rho * o.area, 1e-12f); + ghost.velocity = o.left ? (mdot_out >= 0.0f ? -speed : speed) + : (mdot_out >= 0.0f ? speed : -speed); + if (o.left) pipe.left_ghost = ghost; else pipe.right_ghost = ghost; + o.last_mass_flow = mdot_out; + o.last_face_pressure = ghost.pressure; + } +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/HelmholtzScenario.hpp b/include/FluidSim/HelmholtzScenario.hpp new file mode 100644 index 0000000..e9eae17 --- /dev/null +++ b/include/FluidSim/HelmholtzScenario.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include "FluidSim\AudioProcessing.hpp" +#include "FluidSim\Fluid.hpp" +#include "FluidSim\Solver.hpp" +#include + +namespace fluidsim::port { + +class HelmholtzScenario { +public: + explicit HelmholtzScenario(int sample_rate = 44100); + float process(); + float cavity_pressure() const { return cavity_.pressure(); } + float open_end_flow() const { return boundaries_.open_end_mass_flow(0); } + const Pipe1D& neck() const { return pipes_.front(); } + +private: + Volume0D cavity_{1e-3f, 1.1f * 101325.0f, 300.0f}; + Port* cavity_port_{}; + std::vector pipes_; + float dt_{}, neck_area_{}, cavity_orifice_area_{}; + BoundarySystem boundaries_; + Solver solver_; + SoundProcessor sound_; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/MechanicalSoundProcessor.hpp b/include/FluidSim/MechanicalSoundProcessor.hpp new file mode 100644 index 0000000..4c48027 --- /dev/null +++ b/include/FluidSim/MechanicalSoundProcessor.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "FluidSim/Engine.hpp" +#include + +namespace fluidsim::port { + +// Synthetic valve, piston, combustion, chain, and gear sounds driven by live +// engine state; direct port of Core/MechanicalSoundProcessor.cs. +class MechanicalSoundProcessor { +public: + float master_gain{0.6f}, valve_tick_gain{0.9f}, piston_slap_gain{1.2f}; + float combustion_gain{0.5f}, gear_whine_gain{0.12f}, cam_chain_gain{0.25f}; + + explicit MechanicalSoundProcessor(int sample_rate = 44100); + float process(const Crankshaft& crank, const EngineCylinder& cylinder, float load); + +private: + class Resonator { + public: + Resonator(float frequency_hz, float q, float sample_rate); + float process(float input); + private: + float b1_{}, b2_{}, gain_{}, y1_{}, y2_{}; + }; + class OnePoleLowPass { + public: + OnePoleLowPass(float cutoff_hz, float sample_rate); + float process(float input); + private: + float alpha_{}, state_{}; + }; + + float dt_{}; + Resonator intake_tick_, exhaust_tick_; + OnePoleLowPass slap_filter_, combustion_filter_, chain_filter_; + float slap_envelope_{}, combustion_envelope_{}, chain_envelope_{}; + float slap_decay_{}, combustion_decay_{}, chain_decay_{}; + float gear_phase1_{}, gear_phase2_{}, gear_phase3_{}; + float previous_intake_area_{}, previous_exhaust_area_{}, previous_pressure_{}; + float previous_cycle_position_{-1.0f}, previous_crank_fraction_{-1.0f}; + std::mt19937 random_{0x5EEDu}; + std::uniform_real_distribution unit_{-1.0f, 1.0f}; + + float random_sample(); + static bool crossed_fraction(float previous, float current, float target); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/OutdoorExhaustReverb.hpp b/include/FluidSim/OutdoorExhaustReverb.hpp new file mode 100644 index 0000000..7289d40 --- /dev/null +++ b/include/FluidSim/OutdoorExhaustReverb.hpp @@ -0,0 +1,55 @@ +#pragma once + +#include "FluidSim/AudioProcessing.hpp" +#include + +namespace fluidsim::port { + +// Eight-line frequency-dependent feedback delay network with early reflections +// and transient bloom, ported from Core/OutdoorExhaustReverb.cs. +class OutdoorExhaustReverb { +public: + float dry_mix{1.00f}, early_mix{0.10f}, tail_mix{0.28f}; + float pre_delay_ms{5.0f}, ground_delay_ms{9.0f}, wall_delay_ms{45.0f}; + float low_decay_seconds{1.50f}, mid_decay_seconds{0.70f}, high_decay_seconds{0.20f}; + float band_low_hz{300.0f}, band_high_hz{3000.0f}; + float transient_feedback_boost{0.35f}, transient_hf_boost{3.0f}; + float transient_decay_ms{110.0f}, saturation_drive{1.2f}; + + explicit OutdoorExhaustReverb(int sample_rate = 44100); + AudioFrame process_stereo(float input); + float process(float input); + +private: + class DelayLine { + public: + explicit DelayLine(std::size_t size = 1); + float read() const; + float read_write(float value); + void write(float value); + private: + std::vector buffer_; + std::size_t position_{}; + }; + class OnePoleLowPass { + public: + OnePoleLowPass(float cutoff_hz, int sample_rate); + float process(float value); + private: + float alpha_{}, state_{}; + }; + + int sample_rate_{}; + float dt_{}, gain_low_{}, gain_mid_{}, gain_high_{}; + std::vector fdn_; + std::vector low_filters_, mid_filters_; + std::vector read_buffer_, mix_buffer_; + DelayLine pre_delay_, ground_left_, ground_right_, wall_left_, wall_right_; + OnePoleLowPass ground_filter_left_, ground_filter_right_, wall_filter_left_, wall_filter_right_; + float fast_envelope_{}, slow_envelope_{}, bloom_{}; + float fast_attack_{}, fast_release_{}, slow_attack_{}, slow_release_{}, bloom_decay_{}; + + static std::size_t delay_samples(int sample_rate, float milliseconds); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/PulseTestScenario.hpp b/include/FluidSim/PulseTestScenario.hpp new file mode 100644 index 0000000..9ea7bb2 --- /dev/null +++ b/include/FluidSim/PulseTestScenario.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include "FluidSim/AudioProcessing.hpp" +#include "FluidSim/Fluid.hpp" +#include "FluidSim/Solver.hpp" +#include + +namespace fluidsim::port { + +class PulseTestScenario { +public: + explicit PulseTestScenario(int sample_rate = 44100); + float process(); + const Pipe1D& pipe() const { return pipes_.front(); } + +private: + std::vector pipes_; + float dt_{}, pipe_area_{}; + BoundarySystem boundaries_; + Solver solver_; + SoundProcessor sound_; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Scenario.hpp b/include/FluidSim/Scenario.hpp new file mode 100644 index 0000000..46a4f6a --- /dev/null +++ b/include/FluidSim/Scenario.hpp @@ -0,0 +1,82 @@ +#pragma once + +#include "FluidSim/AudioProcessing.hpp" +#include "FluidSim/Engine.hpp" +#include "FluidSim/Fluid.hpp" +#include "FluidSim/MechanicalSoundProcessor.hpp" +#include "FluidSim/OutdoorExhaustReverb.hpp" +#include "FluidSim/Solver.hpp" +#include "FluidSim/Vehicle.hpp" +#include +#include +#include +#include + +namespace fluidsim::port { + +struct ScenarioSnapshot { + float rpm{}, speed_kmh{}, kinetic_energy{}, average_power{}, average_torque{}; + float cylinder_pressure1{}, cylinder_pressure2{}; + float cylinder_temperature1{}, cylinder_temperature2{}; + float cylinder_piston_fraction1{}, cylinder_piston_fraction2{}; + float valve1_intake_area{}, valve1_exhaust_area{}, valve2_intake_area{}, valve2_exhaust_area{}; + float intake_pressure{}, exhaust_pressure{}; + float throttle{}, clutch{}, intake_flow{}, exhaust_flow{}; + int gear{}; + bool rev_limit1{}, rev_limit2{}, mechanical_sound_enabled{}; + std::array, 6> pipe_states; +}; + +// Active C# SingleCylScenario topology: a 600 cc parallel twin, independent +// intake runners, separate short exhaust headers, shared collector and tailpipe. +class FourStrokeTwinScenario { +public: + explicit FourStrokeTwinScenario(int sample_rate = 44100); + AudioFrame process(); + void set_controls(float throttle_value, float load_value, float clutch_value); + void toggle_mechanical_sound(); + void shift_up(); + void shift_down(); + void set_neutral(); + void set_first_gear(); + ScenarioSnapshot snapshot() const; + + float throttle{}; + float load{}; + float clutch{1.0f}; + bool mechanical_sound_enabled{}; + float mechanical_mix{0.02f}, reverb_send_gain{1.0f}; + float exhaust_send{0.5f}, intake_send{0.5f}, mechanical_send{}; + Crankshaft crankshaft{}; + FourStrokeCylinder cylinder1{}, cylinder2{}; + Vehicle vehicle{}; + Volume0D intake_plenum{1.0e-3f, 101325.0f, 300.0f}; + Volume0D exhaust_collector{0.15e-3f, 101325.0f, 700.0f}; + std::vector pipes; + float intake_flow{}, exhaust_flow{}; + float last_mix_l{}, last_mix_r{}; + std::uint64_t sample_count{}; + +private: + Port* throttle_port_{}; + Port* runner_feed1_{}; + Port* runner_feed2_{}; + Port* collector_in1_{}; + Port* collector_in2_{}; + Port* collector_out_{}; + float throttle_area_{}, runner_area_{}, valve1_in_{}, valve2_in_{}; + float valve1_out_{}, valve2_out_{}, header_area_{}, tail_area_{}; + float dt_{}; + BoundarySystem boundaries_; + Solver solver_; + RadiationSource intake_radiation_; + RadiationSource exhaust_radiation_; + OutdoorExhaustReverb reverb_; + MechanicalSoundProcessor mechanical_sound_; + mutable std::mutex mutex_; + + static float area_for_diameter(float diameter); + static std::vector create_pipes(); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/SimulationRingBuffer.hpp b/include/FluidSim/SimulationRingBuffer.hpp new file mode 100644 index 0000000..2a68dc8 --- /dev/null +++ b/include/FluidSim/SimulationRingBuffer.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fluidsim::port { + +class SimulationRingBuffer { +public: + explicit SimulationRingBuffer(std::size_t capacity); + std::size_t available_samples() const; + std::size_t free_space() const; + bool write(float left, float right, std::stop_token stop = {}); + bool wait_for_space(std::size_t target_fill, std::stop_token stop); + void reset_consumer(); + std::size_t read_interpolated(std::span left, std::span right, double speed); + +private: + std::vector left_, right_; + std::size_t mask_{}; + std::atomic write_head_{}; + std::atomic read_head_{}; + double read_position_{}; + bool consumer_initialized_{}; + mutable std::mutex wait_mutex_; + std::condition_variable space_available_; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Solver.hpp b/include/FluidSim/Solver.hpp new file mode 100644 index 0000000..11248e4 --- /dev/null +++ b/include/FluidSim/Solver.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include "FluidSim/Component.hpp" +#include "FluidSim/Fluid.hpp" +#include +#include + +namespace fluidsim::port { + +class Solver { +public: + Solver(std::vector& pipes, BoundarySystem& boundaries); + int minimum_substeps{4}; + void add_component(Component& component); + void step(float dt); + std::uint64_t step_count() const { return step_count_; } +private: + std::vector& pipes_; + BoundarySystem& boundaries_; + std::vector components_; + std::uint64_t step_count_{}; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/SoundEngine.hpp b/include/FluidSim/SoundEngine.hpp new file mode 100644 index 0000000..95a6e28 --- /dev/null +++ b/include/FluidSim/SoundEngine.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include "FluidSim/AudioOutputStream.hpp" + +namespace fluidsim::port { + +// Lifecycle wrapper matching the C# Audio.SoundEngine. The output stream owns +// the device callback, while this class gives the application one explicit +// start/stop and playback-control surface. +class SoundEngine { +public: + explicit SoundEngine(SimulationRingBuffer& source_buffer); + ~SoundEngine(); + + SoundEngine(const SoundEngine&) = delete; + SoundEngine& operator=(const SoundEngine&) = delete; + + void start(); + void stop(); + void set_speed(double speed); + double speed() const; + void set_volume(float volume); + float volume() const; + bool is_playing() const { return is_playing_; } + +private: + AudioOutputStream stream_; + bool is_playing_{}; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/ThreadLoadTracker.hpp b/include/FluidSim/ThreadLoadTracker.hpp new file mode 100644 index 0000000..f379ab2 --- /dev/null +++ b/include/FluidSim/ThreadLoadTracker.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include + +namespace fluidsim::port { + +class ThreadLoadTracker { +public: + void record(double busy_milliseconds, double budget_milliseconds) { + if (budget_milliseconds <= 0.0) return; + constexpr double alpha = 0.05; + const double instant = busy_milliseconds / budget_milliseconds * 100.0; + const double previous = load_percent_.load(std::memory_order_relaxed); + load_percent_.store(previous + alpha * (instant - previous), std::memory_order_relaxed); + last_busy_ms_.store(busy_milliseconds, std::memory_order_relaxed); + } + double load_percent() const { return load_percent_.load(std::memory_order_relaxed); } + double last_step_microseconds() const { return last_busy_ms_.load(std::memory_order_relaxed) * 1000.0; } +private: + std::atomic load_percent_{}; + std::atomic last_busy_ms_{}; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/TwoStrokeCylinder.hpp b/include/FluidSim/TwoStrokeCylinder.hpp new file mode 100644 index 0000000..d511376 --- /dev/null +++ b/include/FluidSim/TwoStrokeCylinder.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "FluidSim/Crankcase.hpp" + +namespace fluidsim::port { + +// Two-stroke variant of the reciprocating cylinder model. The shared four- +// stroke implementation supplies geometry, state balance, and combustion +// parameters; cycle events and pumping torque are specialized here. +class TwoStrokeCylinder final : public EngineCylinder { +public: + TwoStrokeCylinder(float bore, float stroke, float connecting_rod_length, + float compression_ratio, float transfer_duration_degrees, + float exhaust_duration_degrees); + + float transfer_open() const { return 180.0f - transfer_duration_ * 0.5f; } + float transfer_close() const { return 180.0f + transfer_duration_ * 0.5f; } + float exhaust_open() const { return 180.0f - exhaust_duration_ * 0.5f; } + float exhaust_close() const { return 180.0f + exhaust_duration_ * 0.5f; } + void set_crankcase(Crankcase& crankcase) { crankcase_ = &crankcase; } + + float valve_area(float cycle_degrees, bool inlet) const override; + void pre_step(Crankshaft& crank, float dt, float throttle) override; + +private: + float transfer_duration_{}, exhaust_duration_{}; + Crankcase* crankcase_{}; + + static float lift(float angle, float opens, float closes, float peak_lift); + void handle_cycle_events(float previous_degrees, float current_degrees); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/TwoStrokeScenario.hpp b/include/FluidSim/TwoStrokeScenario.hpp new file mode 100644 index 0000000..b919dde --- /dev/null +++ b/include/FluidSim/TwoStrokeScenario.hpp @@ -0,0 +1,53 @@ +#pragma once + +#include "FluidSim/AudioProcessing.hpp" +#include "FluidSim/Crankcase.hpp" +#include "FluidSim/Fluid.hpp" +#include "FluidSim/OutdoorExhaustReverb.hpp" +#include "FluidSim/Solver.hpp" +#include "FluidSim/TwoStrokeCylinder.hpp" +#include +#include + +namespace fluidsim::port { + +struct TwoStrokeSnapshot { + float rpm{}, cylinder_pressure{}, crankcase_pressure{}, intake_pressure{}; + float throttle{}, transfer_area{}, exhaust_area{}; + bool reed_open{}; +}; + +// Port of the C# 125 cc crankcase-scavenged two-stroke scenario. It retains +// the plenum / reed / crankcase / transfer / cylinder / exhaust topology. +class TwoStrokeScenario { +public: + explicit TwoStrokeScenario(int sample_rate = 44100); + AudioFrame process(); + void set_throttle(float value); + TwoStrokeSnapshot snapshot() const; + + float throttle{}; + Crankshaft crankshaft{}; + TwoStrokeCylinder cylinder; + Crankcase crankcase; + Volume0D intake_plenum{0.5e-3f, 101325.0f, 300.0f}; + std::vector pipes; + +private: + float dt_{}, max_throttle_area_{}; + float reed_in_area_{}, reed_out_area_{}, transfer_pipe_area_{}, transfer_valve_area_{}, exhaust_valve_area_{}; + bool reed_open_{}; + std::uint64_t step_count_{}; + Port* plenum_inlet_{}; + Port* plenum_outlet_{}; + BoundarySystem boundaries_; + Solver solver_; + SoundProcessor intake_sound_; + SoundProcessor exhaust_sound_; + OutdoorExhaustReverb reverb_; + + static float area_for_diameter(float diameter); + static std::vector create_pipes(); +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Units.hpp b/include/FluidSim/Units.hpp new file mode 100644 index 0000000..d3b7d47 --- /dev/null +++ b/include/FluidSim/Units.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include "FluidSim/Constants.hpp" + +namespace fluidsim::port::units { + +inline constexpr double mm = 1e-3; +inline constexpr double cm = 1e-2; +inline constexpr double inch = 0.0254; +inline constexpr double mm2 = 1e-6; +inline constexpr double cm2 = 1e-4; +inline constexpr double milliliter = 1e-6; +inline constexpr double liter = 1e-3; +inline constexpr double pascal = 1.0; +inline constexpr double kilopascal = 1e3; +inline constexpr double bar = 1e5; +inline constexpr double atmosphere = 101325.0; +inline constexpr double psi = 6894.76; +inline constexpr double gram = 1e-3; +inline constexpr double pound = 0.453592; + +inline constexpr double celsius(double value) { return value + 273.15; } +inline constexpr double area_from_radius(double radius) { return static_cast(pi) * radius * radius; } +inline constexpr double area_from_diameter(double diameter) { return static_cast(pi) * 0.25 * diameter * diameter; } + +} // namespace fluidsim::port::units diff --git a/include/FluidSim/Vehicle.hpp b/include/FluidSim/Vehicle.hpp new file mode 100644 index 0000000..3633c80 --- /dev/null +++ b/include/FluidSim/Vehicle.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include + +namespace fluidsim::port { + +// Longitudinal motorcycle model and the friction clutch coupling it to the crank. +class Vehicle { +public: + int current_gear{}; + std::array gear_ratios{4.1f, 2.7f, 2.0f, 1.55f, 1.25f, 0.98f}; + float final_drive_ratio{2.7f}, primary_reduction{1.8f}; + float clutch_input{1.0f}, clutch_capacity{80.0f}, clutch_engagement_rate{8.0f}; + float mass{160.0f}, wheel_radius{0.32f}, rotational_inertia{1.2f}; + float rolling_resistance{0.015f}, drag_coefficient{0.7f}, frontal_area{0.65f}; + float air_density{1.225f}, gravity{9.80665f}; + float speed{}, engagement{}; + + float total_gear_ratio() const; + float kinetic_energy() const; + float update(float engine_omega, float engine_inertia, float dt); + void shift_up(); + void shift_down(); + void set_neutral(); + void set_first_gear(); + +private: + float shift_cut_timer_{}; + void apply_road_resistance(float dt); +}; + +} // namespace fluidsim::port diff --git a/src/audio_output_stream.cpp b/src/audio_output_stream.cpp new file mode 100644 index 0000000..34d5cd5 --- /dev/null +++ b/src/audio_output_stream.cpp @@ -0,0 +1,36 @@ +#include "FluidSim/AudioOutputStream.hpp" + +#include +#include + +namespace fluidsim::port { + +AudioOutputStream::AudioOutputStream(SimulationRingBuffer& ring, unsigned int sample_rate) + : ring_(ring) { + initialize(2, sample_rate, {sf::SoundChannel::FrontLeft, sf::SoundChannel::FrontRight}); +} + +bool AudioOutputStream::onGetData(Chunk& data) { + const std::size_t read = ring_.read_interpolated(temp_left_, temp_right_, speed()); + constexpr float dc_coefficient = 0.9995f; + for (std::size_t i = 0; i < read; ++i) { + const float left = temp_left_[i] - previous_input_left_ + dc_coefficient * previous_output_left_; + const float right = temp_right_[i] - previous_input_right_ + dc_coefficient * previous_output_right_; + previous_input_left_ = temp_left_[i]; + previous_input_right_ = temp_right_[i]; + previous_output_left_ = left; + previous_output_right_ = right; + const float finite_left = std::isfinite(left) ? std::clamp(left, -1.0f, 1.0f) : 0.0f; + const float finite_right = std::isfinite(right) ? std::clamp(right, -1.0f, 1.0f) : 0.0f; + output_[i * 2] = static_cast(finite_left * 32767.0f); + output_[i * 2 + 1] = static_cast(finite_right * 32767.0f); + } + std::fill(output_.begin() + static_cast(read * 2), output_.end(), 0); + data.samples = output_.data(); + data.sampleCount = output_.size(); + return true; +} + +void AudioOutputStream::onSeek(sf::Time /*time_offset*/) {} + +} // namespace fluidsim::port diff --git a/src/audio_processing.cpp b/src/audio_processing.cpp new file mode 100644 index 0000000..a132db9 --- /dev/null +++ b/src/audio_processing.cpp @@ -0,0 +1,88 @@ +#include "FluidSim/AudioProcessing.hpp" +#include "FluidSim/Engine.hpp" + +#include +#include + +namespace fluidsim::port { + +SoundProcessor::SoundProcessor(int sample_rate, float listener_distance) { + const int safe_rate = std::max(sample_rate, 1); + dt_ = 1.0f / static_cast(safe_rate); + scale_factor_ = 1.0f / (4.0f * pi * std::max(listener_distance, 1e-3f)); + derivative_alpha_ = std::exp(-dt_ / 0.020f); + lowpass_alpha_ = std::exp(-dt_ / 0.005f); +} + +float SoundProcessor::process(float mass_flow_out) { + flow_lowpass_ = lowpass_alpha_ * flow_lowpass_ + (1.0f - lowpass_alpha_) * mass_flow_out; + const float raw_derivative = (flow_lowpass_ - previous_flow_) / dt_; + previous_flow_ = flow_lowpass_; + smooth_derivative_ = derivative_alpha_ * smooth_derivative_ + + (1.0f - derivative_alpha_) * raw_derivative; + return std::tanh(smooth_derivative_ * scale_factor_ * gain); +} + +RadiationSource::RadiationSource(int sample_rate) + : sample_rate_(std::max(sample_rate, 1)), dt_(1.0f / static_cast(std::max(sample_rate, 1))) { + const std::size_t max_delay = static_cast(std::max(1.0f, + 0.5f / std::max(speed_of_sound, 1.0f) * static_cast(sample_rate_))); + left_delay_.resize(max_delay); + right_delay_.resize(max_delay); +} + +AudioFrame RadiationSource::process(float mass_flow) { + const float flow_alpha = 1.0f - std::exp(-dt_ / std::max(flow_tau, 1e-5f)); + const float derivative_decay = std::exp(-dt_ / std::max(smooth_tau, 1e-5f)); + flow_lp_ += flow_alpha * (mass_flow - flow_lp_); + const float raw_derivative = (flow_lp_ - previous_flow_lp_) / dt_; + previous_flow_lp_ = flow_lp_; + smooth_derivative_ = derivative_decay * smooth_derivative_ + (1.0f - derivative_decay) * raw_derivative; + + const Vec3 to_listener = subtract(listener, pipe_end); + const float distance = std::max(length(to_listener), 0.1f); + const Vec3 direction{to_listener.x / distance, to_listener.y / distance, to_listener.z / distance}; + const float cosine = dot(pipe_axis, direction); + const float directivity = (1.0f - directivity_balance) + directivity_balance * cosine; + last_mono_pressure_ = std::tanh(smooth_derivative_ * gain * directivity / (4.0f * 3.14159265358979323846f * distance)); + + const float air_cutoff = std::max(min_air_cutoff, + 20000.0f * std::exp(-distance * air_absorption_exponent)); + const float filtered = one_pole(last_mono_pressure_, air_cutoff, air_filter_); + + const Vec3 ear_offset{listener_right.x * ear_spacing * 0.5f, + listener_right.y * ear_spacing * 0.5f, + listener_right.z * ear_spacing * 0.5f}; + const Vec3 left_ear = subtract(listener, ear_offset); + const Vec3 right_ear{listener.x + ear_offset.x, listener.y + ear_offset.y, listener.z + ear_offset.z}; + const float left_distance = std::max(length(subtract(left_ear, pipe_end)), 0.1f); + const float right_distance = std::max(length(subtract(right_ear, pipe_end)), 0.1f); + left_offset_ = std::min(left_delay_.size() - 1, static_cast(left_distance / speed_of_sound / dt_)); + right_offset_ = std::min(right_delay_.size() - 1, static_cast(right_distance / speed_of_sound / dt_)); + + const auto delay_sample = [this, filtered](std::vector& buffer, std::size_t offset, float ear_distance) { + const std::size_t read = (write_position_ + buffer.size() - offset) % buffer.size(); + const float value = buffer[read]; + buffer[write_position_] = filtered * (length(subtract(listener, pipe_end)) / ear_distance); + return value; + }; + float left = delay_sample(left_delay_, left_offset_, left_distance); + float right = delay_sample(right_delay_, right_offset_, right_distance); + write_position_ = (write_position_ + 1) % left_delay_.size(); + if (left_distance > right_distance) left = one_pole(left, head_shadow_cutoff, shadow_left_); + else if (right_distance > left_distance) right = one_pole(right, head_shadow_cutoff, shadow_right_); + return {left, right}; +} + +float RadiationSource::dot(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; } +float RadiationSource::length(Vec3 v) { return std::sqrt(v.x * v.x + v.y * v.y + v.z * v.z); } +Vec3 RadiationSource::subtract(Vec3 a, Vec3 b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; } + +float RadiationSource::one_pole(float input, float cutoff, float& state) const { + const float angular = 2.0f * 3.14159265358979323846f * cutoff / static_cast(sample_rate_); + const float alpha = 1.0f - std::exp(-angular); + state += alpha * (input - state); + return state; +} + +} // namespace fluidsim::port diff --git a/src/crankcase.cpp b/src/crankcase.cpp new file mode 100644 index 0000000..152b9d6 --- /dev/null +++ b/src/crankcase.cpp @@ -0,0 +1,102 @@ +#include "FluidSim/Crankcase.hpp" + +#include +#include + +namespace fluidsim::port { +namespace { +constexpr float cv = gas_constant / (gamma - 1.0f); +} + +Crankcase::Crankcase(Crankshaft& crankshaft, float crank_radius, float connecting_rod_length, + float bore, float clearance_volume, float initial_pressure, + float initial_temperature) + : crankshaft_(crankshaft), crank_radius_(crank_radius), + connecting_rod_length_(connecting_rod_length), + piston_area_(pi * 0.25f * bore * bore), clearance_volume_(clearance_volume), + obliquity_(crank_radius / connecting_rod_length) { + const float theta = crankshaft_.angle; + const float sine = std::sin(theta), cosine = std::cos(theta); + const float term = std::sqrt(std::max(1.0f - obliquity_ * obliquity_ * sine * sine, 0.0f)); + const float displacement = crank_radius_ * (1.0f - cosine) + connecting_rod_length_ * (1.0f - term); + const float maximum_volume = clearance_volume_ + piston_area_ * 2.0f * crank_radius_; + volume_ = maximum_volume - piston_area_ * displacement; + previous_volume_ = volume_; + pressure_ = initial_pressure; + temperature_ = initial_temperature; + mass_ = initial_pressure * volume_ / (gas_constant * initial_temperature); + energy_ = mass_ * cv * initial_temperature; + intake_port.owner_air_mass = transfer_port.owner_air_mass = &mass_; + sync_ports(); +} + +void Crankcase::pre_step(float /*dt*/) { + previous_volume_ = volume_; + const float theta = std::fmod(crankshaft_.angle, 2.0f * pi); + const float sine = std::sin(theta), cosine = std::cos(theta); + const float term = std::sqrt(std::max(1.0f - obliquity_ * obliquity_ * sine * sine, 0.0f)); + const float displacement = crank_radius_ * (1.0f - cosine) + connecting_rod_length_ * (1.0f - term); + const float maximum_volume = clearance_volume_ + piston_area_ * 2.0f * crank_radius_; + volume_ = std::max(maximum_volume - piston_area_ * displacement, 1e-12f); + if (mass_ > 1e-12f) { + temperature_ = energy_ / (mass_ * cv); + pressure_ = mass_ * gas_constant * temperature_ / volume_; + } + sync_ports(); +} + +void Crankcase::update_state(float dt) { + float delta_mass = 0.0f, delta_air_mass = 0.0f, delta_energy = 0.0f; + for (Port* port : {&intake_port, &transfer_port}) { + const float flow = port->mass_flow; + const float fraction = flow >= 0.0f ? port->air_fraction : air_fraction_; + delta_mass += flow * dt; + delta_air_mass += flow * fraction * dt; + delta_energy += flow * port->enthalpy * dt; + } + const float old_air_mass = mass_ * air_fraction_; + mass_ += delta_mass; + energy_ += delta_energy; + if (mass_ > 1e-12f) air_fraction_ = std::clamp((old_air_mass + delta_air_mass) / mass_, 0.0f, 1.0f); + + const float volume_change = volume_ - previous_volume_; + energy_ -= pressure_ * volume_change; + if (mass_ < 1e-9f) { + mass_ = 1e-9f; + energy_ = mass_ * cv * 300.0f; + air_fraction_ = 1.0f; + } + if (energy_ < 0.0f) energy_ = mass_ * cv * 300.0f; + + temperature_ = energy_ / (mass_ * cv); + pressure_ = mass_ * gas_constant * temperature_ / std::max(volume_, 1e-12f); + if (pressure_ < 1.0f) { + temperature_ = volume_ / (mass_ * gas_constant); + energy_ = mass_ * cv * temperature_; + pressure_ = 1.0f; + } + constexpr float max_pressure = 5.0e5f; + if (pressure_ > max_pressure) { + const float target_mass = max_pressure * volume_ / (gas_constant * temperature_); + if (mass_ > target_mass) { + mass_ = target_mass; + energy_ = mass_ * cv * temperature_; + pressure_ = max_pressure; + } + } + sync_ports(); +} + +void Crankcase::sync_ports() { + const float rho = mass_ / std::max(volume_, 1e-12f); + const float enthalpy = gamma / (gamma - 1.0f) * pressure_ / std::max(rho, 1e-12f); + for (Port* port : {&intake_port, &transfer_port}) { + port->pressure = pressure_; + port->density = rho; + port->temperature = temperature_; + port->enthalpy = enthalpy; + port->air_fraction = air_fraction_; + } +} + +} // namespace fluidsim::port diff --git a/src/desktop_main.cpp b/src/desktop_main.cpp new file mode 100644 index 0000000..f837322 --- /dev/null +++ b/src/desktop_main.cpp @@ -0,0 +1,402 @@ +#include "FluidSim/Scenario.hpp" +#include "FluidSim/SoundEngine.hpp" +#include "FluidSim/DynoCurve.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { +using namespace fluidsim::port; + +constexpr unsigned int window_width = 1280; +constexpr unsigned int window_height = 720; +constexpr std::size_t ring_capacity = 8192; +constexpr std::size_t target_fill = 1024; +constexpr std::size_t preroll_samples = 1024; +constexpr float throttle_lerp_rate = 10.0f; + +sf::Color pressure_color(float pressure) { + const float bar = pressure / 1.0e5f; + if (bar < 1.0f) { + const float f = std::clamp(bar, 0.0f, 1.0f); + return {0, static_cast(255.0f * f), static_cast(255.0f * (1.0f - f))}; + } + const float f = std::min((bar - 1.0f) / 9.0f, 1.0f); + return {static_cast(255.0f * f), static_cast(255.0f * (1.0f - f)), 0}; +} + +sf::Color temperature_color(float temperature) { + const float t = std::clamp(temperature, 0.0f, 2000.0f); + if (t < 300.0f) { + const float f = t / 300.0f; + return {0, static_cast(255.0f * f), static_cast(255.0f * (1.0f - f))}; + } + const float f = (t - 300.0f) / 1700.0f; + return {static_cast(255.0f * f), static_cast(255.0f * (1.0f - f)), 0}; +} + +void draw_label(sf::RenderWindow& window, const sf::Font* font, const std::string& text, + sf::Vector2f position, sf::Color color = sf::Color::White, unsigned int size = 14) { + if (!font) return; + sf::Text label(*font); + label.setString(text); + label.setPosition(position); + label.setFillColor(color); + label.setCharacterSize(size); + window.draw(label); +} + +void draw_volume(sf::RenderWindow& window, float pressure, sf::FloatRect bounds) { + sf::RectangleShape body({bounds.size.x, bounds.size.y}); + body.setPosition(bounds.position); + body.setFillColor(pressure_color(pressure)); + body.setOutlineColor(sf::Color::White); + body.setOutlineThickness(1.0f); + window.draw(body); +} + +void draw_pipe(sf::RenderWindow& window, const ScenarioSnapshot& snapshot, std::size_t pipe_index, + float center_y, float x0, float x1) { + const auto& cells = snapshot.pipe_states[pipe_index]; + if (cells.empty()) return; + const float cell_width = (x1 - x0) / static_cast(cells.size()); + for (std::size_t i = 0; i < cells.size(); ++i) { + const auto& cell = cells[i]; + const float thickness = std::clamp(6.0f + cell.pressure / 101325.0f * 0.7f, 6.0f, 20.0f); + sf::RectangleShape segment({cell_width + 1.0f, thickness}); + segment.setPosition({x0 + static_cast(i) * cell_width, center_y - thickness * 0.5f}); + segment.setFillColor(temperature_color(cell.pressure / std::max(cell.rho * 287.0f, 1e-9f))); + window.draw(segment); + } +} + +void draw_cylinder(sf::RenderWindow& window, const sf::Font* font, float pressure, + float piston_fraction, float top, const char* name, bool intake_open, bool exhaust_open) { + constexpr float x = 470.0f, width = 74.0f, height = 142.0f; + const float fraction = std::clamp(piston_fraction, 0.0f, 1.0f); + const float gas_height = height * fraction; + sf::RectangleShape wall({width, height}); + wall.setPosition({x, top}); + wall.setFillColor({52, 59, 68}); + window.draw(wall); + sf::RectangleShape gas({width, gas_height}); + gas.setPosition({x, top}); + gas.setFillColor(pressure_color(pressure)); + window.draw(gas); + sf::RectangleShape piston({width + 8.0f, 5.0f}); + piston.setPosition({x - 4.0f, top + gas_height}); + piston.setFillColor(sf::Color::White); + window.draw(piston); + sf::RectangleShape intake_valve({7.0f, 12.0f}); + intake_valve.setPosition({x - 12.0f, top + 5.0f}); + intake_valve.setFillColor(intake_open ? sf::Color::Green : sf::Color::Red); + window.draw(intake_valve); + sf::RectangleShape exhaust_valve({7.0f, 12.0f}); + exhaust_valve.setPosition({x + width + 5.0f, top + 5.0f}); + exhaust_valve.setFillColor(exhaust_open ? sf::Color::Green : sf::Color::Red); + window.draw(exhaust_valve); + sf::RectangleShape border({width, height}); + border.setPosition({x, top}); + border.setFillColor(sf::Color::Transparent); + border.setOutlineColor(sf::Color(190, 200, 210)); + border.setOutlineThickness(1.0f); + window.draw(border); + draw_label(window, font, name, {x + 13.0f, top - 22.0f}, sf::Color::White, 14); +} + +void draw_graph_line(sf::RenderWindow& window, sf::Vector2f start, sf::Vector2f end, + sf::Color color, float thickness = 1.5f) { + const float dx = end.x - start.x; + const float dy = end.y - start.y; + const float length = std::sqrt(dx * dx + dy * dy); + if (length <= 0.0f) return; + sf::RectangleShape line({length, thickness}); + line.setPosition(start); + line.setRotation(sf::degrees(std::atan2(dy, dx) * 180.0f / pi)); + line.setFillColor(color); + window.draw(line); +} + +void draw_dyno_curve(sf::RenderWindow& window, const sf::Font* font, const DynoCurve& curve, + float graph_x, float graph_y, float graph_width, float graph_height, + float current_rpm, float current_power_kw) { + const auto& bins = curve.bins(); + if (bins.empty()) return; + + float max_power = 0.01f, max_torque = 0.01f, max_rpm = 1000.0f; + for (std::size_t index = 0; index < bins.size(); ++index) { + if (bins[index].power_kw > 0.0f || bins[index].torque_nm > 0.0f) { + max_power = std::max(max_power, bins[index].power_kw); + max_torque = std::max(max_torque, bins[index].torque_nm); + max_rpm = std::max(max_rpm, static_cast(index) * DynoCurve::rpm_bin_size + + DynoCurve::rpm_bin_size * 0.5f); + } + } + max_power *= 1.1f; + max_torque *= 1.1f; + max_rpm = std::max(max_rpm * 1.05f, 1000.0f); + + sf::RectangleShape background({graph_width, graph_height}); + background.setPosition({graph_x, graph_y}); + background.setFillColor({20, 20, 20, 220}); + background.setOutlineColor({70, 78, 90}); + background.setOutlineThickness(1.0f); + window.draw(background); + + constexpr float left_margin = 50.0f, right_margin = 48.0f; + constexpr float top_margin = 24.0f, bottom_margin = 38.0f; + const float plot_x = graph_x + left_margin; + const float plot_y = graph_y + top_margin; + const float plot_width = graph_width - left_margin - right_margin; + const float plot_height = graph_height - top_margin - bottom_margin; + const sf::Color grid_color(50, 50, 50); + for (int index = 0; index <= 9; ++index) { + const float x = plot_x + static_cast(index) / 9.0f * plot_width; + sf::RectangleShape grid({1.0f, plot_height}); + grid.setPosition({x, plot_y}); + grid.setFillColor(grid_color); + window.draw(grid); + } + for (int index = 0; index <= 5; ++index) { + const float y = plot_y + (1.0f - static_cast(index) / 5.0f) * plot_height; + sf::RectangleShape grid({plot_width, 1.0f}); + grid.setPosition({plot_x, y}); + grid.setFillColor(grid_color); + window.draw(grid); + draw_label(window, font, std::to_string(static_cast(max_power * index / 5.0f)), + {graph_x + 3.0f, y - 7.0f}, sf::Color(240, 100, 90), 9); + draw_label(window, font, std::to_string(static_cast(max_torque * index / 5.0f)), + {graph_x + graph_width - 42.0f, y - 7.0f}, sf::Color(100, 150, 255), 9); + } + for (int index = 0; index <= 9; ++index) { + const float rpm = max_rpm * static_cast(index) / 9.0f; + const float x = plot_x + static_cast(index) / 9.0f * plot_width; + std::ostringstream tick; + tick << std::fixed << std::setprecision(1) << rpm / 1000.0f << 'k'; + draw_label(window, font, tick.str(), {x - 12.0f, graph_y + graph_height - 31.0f}, + sf::Color(200, 205, 215), 9); + } + draw_label(window, font, "DYNO", {graph_x + 8.0f, graph_y + 4.0f}, sf::Color::White, 11); + draw_label(window, font, "kW", {graph_x + 8.0f, graph_y + 17.0f}, sf::Color(255, 100, 90), 9); + draw_label(window, font, "Nm", {graph_x + graph_width - 30.0f, graph_y + 17.0f}, + sf::Color(110, 160, 255), 9); + draw_label(window, font, "RPM", {graph_x + graph_width * 0.5f - 12.0f, + graph_y + graph_height - 15.0f}, sf::Color::White, 10); + + const sf::Color power_color(255, 27, 27), torque_color(9, 90, 255); + for (std::size_t index = 1; index < bins.size(); ++index) { + const float rpm0 = (static_cast(index - 1) + 0.5f) * DynoCurve::rpm_bin_size; + const float rpm1 = (static_cast(index) + 0.5f) * DynoCurve::rpm_bin_size; + if (rpm0 > max_rpm || rpm1 > max_rpm) break; + const float x0 = plot_x + rpm0 / max_rpm * plot_width; + const float x1 = plot_x + rpm1 / max_rpm * plot_width; + if (bins[index - 1].power_kw > 0.0f && bins[index].power_kw > 0.0f) { + const float y0 = plot_y + (1.0f - bins[index - 1].power_kw / max_power) * plot_height; + const float y1 = plot_y + (1.0f - bins[index].power_kw / max_power) * plot_height; + draw_graph_line(window, {x0, y0}, {x1, y1}, power_color, 2.0f); + } + if (bins[index - 1].torque_nm > 0.0f && bins[index].torque_nm > 0.0f) { + const float y0 = plot_y + (1.0f - bins[index - 1].torque_nm / max_torque) * plot_height; + const float y1 = plot_y + (1.0f - bins[index].torque_nm / max_torque) * plot_height; + draw_graph_line(window, {x0, y0}, {x1, y1}, torque_color, 2.0f); + } + } + if (current_rpm > 0.0f && current_rpm <= max_rpm && current_power_kw > 0.0f) { + const float x = plot_x + current_rpm / max_rpm * plot_width; + const float y = plot_y + (1.0f - current_power_kw / max_power) * plot_height; + sf::CircleShape dot(3.0f); + dot.setPosition({x - 3.0f, y - 3.0f}); + dot.setFillColor(sf::Color::White); + window.draw(dot); + } +} + +void draw_scene(sf::RenderWindow& window, const sf::Font* font, const ScenarioSnapshot& state, + float throttle_target, float brake, double playback_speed, std::size_t buffer_fill, + std::size_t buffer_capacity, bool time_warp, DynoCurve& dyno_curve) { + constexpr float intake_y = 310.0f; + constexpr float upper_y = 152.0f; + constexpr float lower_y = 400.0f; + window.clear({9, 12, 17}); + const float power_kw = state.average_power / 1000.0f; + dyno_curve.update(state.rpm, power_kw, state.average_torque); + + // Intake: ambient pipe → throttle → plenum → two cylinder runners. + draw_pipe(window, state, 0, intake_y, 36.0f, 165.0f); + sf::RectangleShape throttle({8.0f, 28.0f}); + throttle.setPosition({169.0f, intake_y - 14.0f}); + throttle.setFillColor(sf::Color::Yellow); + window.draw(throttle); + draw_volume(window, state.intake_pressure, {{195.0f, intake_y - 45.0f}, {62.0f, 90.0f}}); + draw_pipe(window, state, 1, upper_y + 71.0f, 263.0f, 458.0f); + draw_pipe(window, state, 2, lower_y + 71.0f, 263.0f, 458.0f); + draw_cylinder(window, font, state.cylinder_pressure1, state.cylinder_piston_fraction1, + upper_y, "CYL 1", state.valve1_intake_area > 1e-7f, state.valve1_exhaust_area > 1e-7f); + draw_cylinder(window, font, state.cylinder_pressure2, state.cylinder_piston_fraction2, + lower_y, "CYL 2", state.valve2_intake_area > 1e-7f, state.valve2_exhaust_area > 1e-7f); + + // Two independent exhaust headers enter a shared collector and tailpipe. + draw_pipe(window, state, 3, upper_y + 71.0f, 556.0f, 820.0f); + draw_pipe(window, state, 4, lower_y + 71.0f, 556.0f, 820.0f); + draw_volume(window, state.exhaust_pressure, {{838.0f, 305.0f}, {70.0f, 82.0f}}); + draw_pipe(window, state, 5, 346.0f, 915.0f, 1236.0f); + + draw_label(window, font, "AMBIENT", {37.0f, intake_y + 22.0f}, {180, 190, 200}, 12); + draw_label(window, font, "THROTTLE", {151.0f, intake_y + 24.0f}, {240, 210, 90}, 11); + draw_label(window, font, "INTAKE PLENUM", {188.0f, intake_y - 62.0f}, sf::Color::White, 11); + draw_label(window, font, "SHARED EXHAUST", {818.0f, 285.0f}, sf::Color::White, 11); + draw_label(window, font, "TAILPIPE", {1084.0f, 362.0f}, {180, 190, 200}, 11); + + std::ostringstream status; + status << (time_warp ? "[Space] normal speed" : "[Space] real time") + << " Sim " << std::fixed << std::setprecision(3) << playback_speed << "x" + << " Audio buffer " << buffer_fill << "/" << buffer_capacity + << "\nBike " << std::setprecision(1) << state.speed_kmh << " km/h Gear " + << (state.gear == 0 ? "N" : std::to_string(state.gear)) + << " Throttle " << state.throttle * 100.0f << "% / target " << throttle_target * 100.0f << "% [W, up/down]" + << " Clutch " << state.clutch * 100.0f << "% [C] Brake " << brake * 100.0f << "% [left/right]" + << " Shift [A/D], Neutral [Z] Mechanical [M]: " + << (state.mechanical_sound_enabled ? "ON" : "OFF"); + if (font) { + sf::RectangleShape panel({1260.0f, 92.0f}); + panel.setPosition({10.0f, 10.0f}); + panel.setFillColor({12, 18, 26, 236}); + panel.setOutlineColor({80, 110, 130}); + panel.setOutlineThickness(1.0f); + window.draw(panel); + draw_label(window, font, status.str(), {22.0f, 20.0f}, sf::Color::White, 14); + draw_label(window, font, + "RPM " + std::to_string(static_cast(state.rpm)) + + " Power " + std::to_string(state.average_power / 1000.0f) + " kW" + + " Bike kinetic energy " + std::to_string(static_cast(state.kinetic_energy)) + " J", + {24.0f, 122.0f}, {220, 230, 240}, 17); + const sf::Color rev = state.rev_limit1 || state.rev_limit2 ? sf::Color::Red : sf::Color(120, 220, 120); + draw_label(window, font, (state.rev_limit1 || state.rev_limit2) ? "REV CUT" : "REV LIMIT 10000 rpm", + {24.0f, 151.0f}, rev, 14); + draw_label(window, font, + "Cylinder pressure " + std::to_string(state.cylinder_pressure1 / 1e5f) + " / " + + std::to_string(state.cylinder_pressure2 / 1e5f) + " bar Intake " + + std::to_string(state.intake_pressure / 1000.0f) + " kPa Exhaust " + + std::to_string(state.exhaust_pressure / 1000.0f) + " kPa", + {24.0f, 176.0f}, {205, 215, 225}, 13); + } + + draw_dyno_curve(window, font, dyno_curve, 860.0f, 455.0f, 400.0f, 250.0f, + state.rpm, power_kw); + window.display(); +} + +} // namespace + +int main(int argc, char** argv) { + using namespace fluidsim::port; + sf::RenderWindow window(sf::VideoMode({window_width, window_height}), "FluidSim C++ - parallel twin"); + window.setVerticalSyncEnabled(false); + window.setFramerateLimit(60); + window.setKeyRepeatEnabled(false); + + sf::Font font; + std::filesystem::path font_path{"fonts/FiraCodeNerdFont-Medium.ttf"}; + if (!std::filesystem::exists(font_path) && argc > 0) { + const auto executable_font = std::filesystem::path(argv[0]).parent_path() / + "fonts/FiraCodeNerdFont-Medium.ttf"; + if (std::filesystem::exists(executable_font)) font_path = executable_font; + } + const bool has_font = font.openFromFile(font_path.string()); + const sf::Font* font_ptr = has_font ? &font : nullptr; + + FourStrokeTwinScenario scenario(44100); + SimulationRingBuffer ring(ring_capacity); + for (std::size_t i = 0; i < preroll_samples; ++i) { + const AudioFrame audio = scenario.process(); + ring.write(audio.left, audio.right); + } + ring.reset_consumer(); + + SoundEngine audio(ring); + audio.set_volume(100.0f); + audio.start(); + std::jthread simulation([&](std::stop_token stop) { + while (ring.wait_for_space(target_fill, stop)) { + const AudioFrame frame = scenario.process(); + if (!ring.write(frame.left, frame.right, stop)) break; + } + }); + + double desired_speed = 0.001; + double current_speed = desired_speed; + double last_normal_speed = 0.1; + bool time_warp = false; + bool real_time = false; + float throttle_target = 1.0f, throttle_current = 0.0f; + float brake_target = 0.0f, brake_current = 0.0f; + float clutch_current = 0.0f; + bool throttle_held = false, clutch_held = false; + DynoCurve dyno_curve; + auto previous_frame = std::chrono::steady_clock::now(); + + while (window.isOpen()) { + while (const std::optional event = window.pollEvent()) { + if (event->is()) window.close(); + if (const auto* key = event->getIf()) { + switch (key->code) { + case sf::Keyboard::Key::Space: + time_warp = !time_warp; + if (!time_warp) { desired_speed = last_normal_speed; real_time = false; } + break; + case sf::Keyboard::Key::W: throttle_held = true; break; + case sf::Keyboard::Key::C: clutch_held = true; break; + case sf::Keyboard::Key::Up: throttle_target = std::min(1.0f, throttle_target + 0.05f); break; + case sf::Keyboard::Key::Down: throttle_target = std::max(0.0f, throttle_target - 0.05f); break; + case sf::Keyboard::Key::Left: brake_target = std::max(0.0f, brake_target - 0.05f); break; + case sf::Keyboard::Key::Right: brake_target = std::min(1.0f, brake_target + 0.05f); break; + case sf::Keyboard::Key::D: scenario.shift_up(); break; + case sf::Keyboard::Key::A: scenario.shift_down(); break; + case sf::Keyboard::Key::Z: scenario.set_neutral(); break; + case sf::Keyboard::Key::M: scenario.toggle_mechanical_sound(); break; + default: break; + } + } + if (const auto* key = event->getIf()) { + if (key->code == sf::Keyboard::Key::W) throttle_held = false; + if (key->code == sf::Keyboard::Key::C) clutch_held = false; + } + if (const auto* wheel = event->getIf(); wheel && !time_warp) { + desired_speed *= wheel->delta > 0.0f ? 1.1 : (wheel->delta < 0.0f ? 1.0 / 1.1 : 1.0); + desired_speed = std::clamp(desired_speed, 0.0001, 1.0); + last_normal_speed = desired_speed; + real_time = std::abs(desired_speed - 1.0) < 1e-6; + } + } + + const auto now = std::chrono::steady_clock::now(); + const float frame_dt = std::chrono::duration(now - previous_frame).count(); + previous_frame = now; + const double target_speed = time_warp ? 1.0 : desired_speed; + current_speed += (target_speed - current_speed) * (1.0 - std::exp(-8.0 * frame_dt)); + audio.set_speed(current_speed); + + const float wanted_throttle = throttle_held ? throttle_target : 0.0f; + if (wanted_throttle == 0.0f) throttle_current = 0.0f; + else throttle_current += (wanted_throttle - throttle_current) * (1.0f - std::exp(-throttle_lerp_rate * frame_dt)); + brake_current += (brake_target - brake_current) * (1.0f - std::exp(-throttle_lerp_rate * frame_dt)); + const float wanted_clutch = clutch_held ? 1.0f : 0.0f; + clutch_current += (wanted_clutch - clutch_current) * (1.0f - std::exp(-throttle_lerp_rate * frame_dt)); + scenario.set_controls(throttle_current, brake_current, clutch_current); + + draw_scene(window, font_ptr, scenario.snapshot(), throttle_target, brake_current, + current_speed, ring.available_samples(), ring_capacity, time_warp, dyno_curve); + } + + simulation.request_stop(); + audio.stop(); + return 0; +} diff --git a/src/engine.cpp b/src/engine.cpp new file mode 100644 index 0000000..328422c --- /dev/null +++ b/src/engine.cpp @@ -0,0 +1,261 @@ +#include "FluidSim/Engine.hpp" + +#include +#include + +namespace fluidsim::port { + +void Crankshaft::add_torque(float torque) { + external_torque_ += torque; +} + +void Crankshaft::step(float dt) { + if (!std::isfinite(omega)) omega = 0.0f; + if (!std::isfinite(external_torque_)) external_torque_ = 0.0f; + previous_angle = angle; + const float sign = omega > 0.0f ? 1.0f : omega < 0.0f ? -1.0f : 0.0f; + const float friction = friction_constant * sign + friction_viscous * omega; + const float net_torque = external_torque_ - friction; + const float total_torque = net_torque - std::max(load_torque, 0.0f) + drivetrain_torque; + omega = std::max(omega + total_torque / std::max(inertia, 1e-6f) * dt, 0.0f); + angle = std::fmod(angle + omega * dt, 4.0f * pi); + const float instant_power = total_torque * omega; + if (power_count_ == power_buffer_.size()) power_sum_ -= power_buffer_[power_index_]; + else ++power_count_; + power_buffer_[power_index_] = instant_power; + power_sum_ += instant_power; + power_index_ = (power_index_ + 1) % power_buffer_.size(); + average_power = power_sum_ / static_cast(power_count_); + + if (torque_count_ == torque_buffer_.size()) torque_sum_ -= torque_buffer_[torque_index_]; + else ++torque_count_; + torque_buffer_[torque_index_] = total_torque; + torque_sum_ += total_torque; + torque_index_ = (torque_index_ + 1) % torque_buffer_.size(); + average_torque = torque_sum_ / static_cast(torque_count_); + external_torque_ = 0.0f; + drivetrain_torque = 0.0f; +} + +float EngineCylinder::swept_volume() const { + return pi * 0.25f * bore * bore * stroke; +} + +float EngineCylinder::clearance_volume() const { + return swept_volume() / (compression - 1.0f); +} + +float EngineCylinder::total_mass() const { + return air_mass + exhaust_mass; +} + +float EngineCylinder::temperature() const { + const float density = total_mass() / std::max(volume, 1e-12f); + return pressure / std::max(density * gas_constant, 1e-12f); +} + +float EngineCylinder::air_fraction() const { + return air_mass / std::max(total_mass(), 1e-12f); +} + +float EngineCylinder::compute_volume(float theta) const { + const float radius = stroke * 0.5f; + const float obliquity = radius / rod; + const float sine = std::sin(theta), cosine = std::cos(theta); + const float term = std::sqrt(std::max(1.0f - obliquity * obliquity * sine * sine, 0.0f)); + const float piston_travel = radius * (1.0f - cosine) + rod * (1.0f - term); + return clearance_volume() + pi * 0.25f * bore * bore * piston_travel; +} + +void EngineCylinder::initialize(const Crankshaft& /*crank*/) { + // The C# EngineCylinder constructor starts every cylinder at TDC with its + // clearance volume. PhaseOffset is applied on the first PreStep, after the + // conserved state has been initialized; using phase here would silently + // change the initial state of phased cylinders (notably the second twin). + volume = clearance_volume(); + air_mass = 1.225f * volume; + exhaust_mass = 0.0f; + energy = 101325.0f * volume / (gamma - 1.0f); + pressure = 101325.0f; + intake.owner_air_mass = exhaust.owner_air_mass = &air_mass; + intake.owner_exhaust_mass = exhaust.owner_exhaust_mass = &exhaust_mass; + sync_ports(); +} + +float FourStrokeCylinder::valve_area(float cycle_degrees, bool inlet) const { + const float opens = inlet ? intake_open : exhaust_open; + float closes = inlet ? intake_close : exhaust_close; + const float peak_lift = inlet ? intake_valve_lift : exhaust_valve_lift; + const float diameter = inlet ? intake_valve_diameter : exhaust_valve_diameter; + if (closes < opens) closes += 720.0f; + float mapped = std::fmod(cycle_degrees, 720.0f); + if (mapped < 0.0f) mapped += 720.0f; + if (mapped < opens) mapped += 720.0f; + const float duration = closes - opens; + if (mapped < opens || mapped > closes || duration <= 0.0f) return 0.0f; + + const float ramp = duration * 0.25f; + const float hold = duration - 2.0f * ramp; + float lift = 0.0f; + if (mapped < opens + ramp) { + const float t = (mapped - opens) / ramp; + lift = peak_lift * t * t * (3.0f - 2.0f * t); + } else if (mapped < opens + ramp + hold) { + lift = peak_lift; + } else if (ramp > 0.0f) { + const float t = (mapped - (opens + ramp + hold)) / ramp; + lift = peak_lift * (1.0f - t) * (1.0f - t) * (1.0f + 2.0f * t); + } + return pi * diameter * lift; +} + +void FourStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle*/) { + const float rpm = crank.omega * 60.0f / (2.0f * pi); + spark_advance = std::clamp(10.0f + rpm * 0.002f, 5.0f, 40.0f); + update_rev_limiter(rpm, dt); + + const float theta = crank.angle + phase; + const float previous_volume = volume; + volume = compute_volume(theta); + pressure = (gamma - 1.0f) * energy / std::max(volume, 1e-12f); + const float dvolume = volume - previous_volume; + const float radius = stroke * 0.5f, obliquity = radius / rod; + const float sine = std::sin(theta), cosine = std::cos(theta); + const float term = std::sqrt(std::max(1.0f - obliquity * obliquity * sine * sine, 1e-12f)); + const float dx_dtheta = radius * sine * (1.0f + obliquity * cosine / term); + const float piston_area = pi * 0.25f * bore * bore; + crank.add_torque((pressure - 101325.0f) * piston_area * dx_dtheta); + energy -= pressure * dvolume; + + const float previous_degrees = wrap_degrees((crank.previous_angle + phase) * 180.0f / pi); + const float current_degrees = wrap_degrees(theta * 180.0f / pi); + if (crossed_window_close(previous_degrees, current_degrees, intake_open, intake_close)) { + fuel_mass = air_mass / std::max(afr, 1e-6f); + fuel_injected = true; + } + + const float spark_angle = wrap_degrees(720.0f - spark_advance); + const bool crossed_spark = (previous_degrees < spark_angle && current_degrees >= spark_angle) || + (previous_degrees > spark_angle && current_degrees < spark_angle); + if (crossed_spark && fuel_injected && !combustion_active && !rev_limit_active && random_unit() >= misfire_probability) { + combustion_active = true; + burn_fraction = 0.0f; + energy_factor_ = 1.0f + energy_variation * (2.0f * random_unit() - 1.0f); + } + + if (combustion_active) { + float angle_since_spark = current_degrees - spark_angle; + if (angle_since_spark < 0.0f) angle_since_spark += 720.0f; + const float phi = (angle_since_spark - burn_start) / std::max(burn_duration, 1e-6f); + const float burned = phi <= 0.0f ? 0.0f : 1.0f - std::exp(-wiebe_a * std::pow(phi, wiebe_m + 1.0f)); + float fraction = std::max(burn_fraction, std::clamp(burned, 0.0f, 1.0f)); + if (fraction >= 1.0f || angle_since_spark > burn_duration + burn_start + spark_advance) { + fraction = 1.0f; + combustion_active = false; + exhaust_mass += air_mass; + air_mass = 0.0f; + } + fuel_injected = false; + const float increment = fraction - burn_fraction; + if (increment > 0.0f) { + energy += fuel_mass * fuel_lhv * energy_factor_ * increment; + exhaust_mass += fuel_mass * increment; + burn_fraction = fraction; + } + } + + if (rev_limit_active) { + fuel_mass = 0.0f; + fuel_injected = false; + combustion_active = false; + } + const float heat_loss = heat_transfer * wall_area * (temperature() - ambient_temperature) * dt; + energy -= heat_loss; + sync_ports(); +} + +void EngineCylinder::update_state(float dt) { + float delta_air = 0.0f, delta_exhaust = 0.0f, delta_energy = 0.0f; + for (Port* port : {&intake, &exhaust}) { + const float mass_flow = port->mass_flow; + const float fraction = mass_flow >= 0.0f ? port->air_fraction : air_fraction(); + delta_air += mass_flow * fraction * dt; + delta_exhaust += mass_flow * (1.0f - fraction) * dt; + delta_energy += mass_flow * port->enthalpy * dt; + } + air_mass += delta_air; + exhaust_mass += delta_exhaust; + energy += delta_energy; + + const float current_volume = std::max(volume, 1e-12f); + float current_pressure = (gamma - 1.0f) * energy / current_volume; + if (current_pressure > 200.0e5f) { + energy = 200.0e5f * current_volume / (gamma - 1.0f); + current_pressure = 200.0e5f; + } + const float current_density = total_mass() / current_volume; + const float current_temperature = current_pressure / std::max(current_density * gas_constant, 1e-12f); + if (current_temperature > 3500.0f) { + const float pressure_at_limit = current_density * gas_constant * 3500.0f; + energy = pressure_at_limit * current_volume / (gamma - 1.0f); + } + if (total_mass() < 1e-9f) { + air_mass = 1e-9f; + exhaust_mass = 0.0f; + energy = 101325.0f * current_volume / (gamma - 1.0f); + } else if (energy < 0.0f) { + energy = 101325.0f * current_volume / (gamma - 1.0f); + } + air_mass = std::max(air_mass, 0.0f); + exhaust_mass = std::max(exhaust_mass, 0.0f); + pressure = (gamma - 1.0f) * energy / current_volume; + sync_ports(); +} + +float EngineCylinder::random_unit() { + return unit_(random_); +} + +float EngineCylinder::wrap_degrees(float value, float cycle_degrees) { + value = std::fmod(value, cycle_degrees); + return value < 0.0f ? value + cycle_degrees : value; +} + +bool EngineCylinder::crossed_window_close(float previous, float current, float open, float close) { + return previous >= open && previous < close && current >= close; +} + +void EngineCylinder::update_rev_limiter(float rpm, float dt) { + if (!rev_limit_enabled) { + rev_limit_active = false; + rev_limit_timer_ = 0.0f; + return; + } + if (rev_limit_active) { + rev_limit_timer_ -= dt; + if (rev_limit_timer_ <= 0.0f) { + if (rpm >= rev_limit_rpm - rev_limit_hysteresis) rev_limit_timer_ = rev_limit_cut_time; + else rev_limit_active = false; + } + } else if (rpm >= rev_limit_rpm) { + rev_limit_active = true; + rev_limit_timer_ = rev_limit_cut_time; + } +} + +void EngineCylinder::sync_ports() { + const float density = total_mass() / std::max(volume, 1e-12f); + pressure = (gamma - 1.0f) * energy / std::max(volume, 1e-12f); + const float temp = pressure / std::max(density * gas_constant, 1e-12f); + const float fraction = air_fraction(); + const float enthalpy = gamma / (gamma - 1.0f) * pressure / std::max(density, 1e-12f); + for (Port* port : {&intake, &exhaust}) { + port->pressure = pressure; + port->density = density; + port->temperature = temp; + port->air_fraction = fraction; + port->enthalpy = enthalpy; + } +} + +} // namespace fluidsim::port diff --git a/src/helmholtz_scenario.cpp b/src/helmholtz_scenario.cpp new file mode 100644 index 0000000..6b050ad --- /dev/null +++ b/src/helmholtz_scenario.cpp @@ -0,0 +1,25 @@ +#include "FluidSim/HelmholtzScenario.hpp" + +#include + +namespace fluidsim::port { + +HelmholtzScenario::HelmholtzScenario(int sample_rate) + : pipes_{}, dt_(1.0f / static_cast(std::max(sample_rate, 1))), + neck_area_(pi * 0.25f * 0.020f * 0.020f), cavity_orifice_area_(neck_area_), + boundaries_(pipes_), solver_(pipes_, boundaries_), sound_(sample_rate, 1.0f) { + pipes_.emplace_back(neck_area_, 0.050f, 20); + cavity_port_ = &cavity_.create_port(); + boundaries_.add_orifice_with_inertance(*cavity_port_, 0, true, cavity_orifice_area_, + 0.9f, 0.050f); + boundaries_.add_open_end(0, false, 101325.0f, neck_area_); + solver_.minimum_substeps = 8; + solver_.add_component(cavity_); +} + +float HelmholtzScenario::process() { + solver_.step(dt_); + return sound_.process(boundaries_.open_end_mass_flow(0)); +} + +} // namespace fluidsim::port diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..8cacea4 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,95 @@ +#include "FluidSim/Scenario.hpp" +#include "FluidSim/HelmholtzScenario.hpp" +#include "FluidSim/PulseTestScenario.hpp" +#include "FluidSim/TwoStrokeScenario.hpp" + +#include +#include +#include +#include +#include +#include + +int main(int argc, char** argv) { + using namespace fluidsim::port; + float duration_seconds = 1.0f; + float throttle = 0.35f; + std::string scenario_name = "twin"; + for (int i = 1; i + 1 < argc; ++i) { + const std::string option = argv[i]; + if (option == "--seconds") duration_seconds = std::max(0.0f, std::strtof(argv[++i], nullptr)); + else if (option == "--throttle") throttle = std::clamp(std::strtof(argv[++i], nullptr), 0.0f, 1.0f); + else if (option == "--scenario") scenario_name = argv[++i]; + } + + constexpr int sample_rate = 44100; + const auto total_samples = static_cast(duration_seconds * sample_rate); + const auto report_interval = static_cast(sample_rate / 10); + if (scenario_name == "twin") { + FourStrokeTwinScenario scenario; + scenario.throttle = throttle; + std::cout << "FluidSimCpp | 600 cc parallel twin | throttle " << std::fixed << std::setprecision(0) + << throttle * 100.0f << "%\n"; + for (std::uint64_t i = 0; i < total_samples; ++i) { + scenario.process(); + if ((i + 1) % report_interval == 0 || i + 1 == total_samples) { + std::cout << std::fixed << std::setprecision(1) + << "t=" << static_cast(i + 1) / sample_rate << " s" + << " rpm=" << scenario.crankshaft.omega * 60.0f / (2.0f * pi) + << " speed=" << scenario.vehicle.speed * 3.6f << " km/h" + << " gear=" << scenario.vehicle.current_gear + << " cyl=" << scenario.cylinder1.pressure / 1000.0f << "/" + << scenario.cylinder2.pressure / 1000.0f << " kPa" + << " intake=" << scenario.intake_plenum.pressure() / 1000.0f << " kPa" + << " exhaust=" << scenario.exhaust_collector.pressure() / 1000.0f << " kPa\n"; + } + } + } else if (scenario_name == "two-stroke") { + TwoStrokeScenario scenario; + scenario.set_throttle(throttle); + std::cout << "FluidSimCpp | 125 cc crankcase-scavenged two-stroke | throttle " + << std::fixed << std::setprecision(0) << throttle * 100.0f << "%\n"; + for (std::uint64_t i = 0; i < total_samples; ++i) { + scenario.process(); + if ((i + 1) % report_interval == 0 || i + 1 == total_samples) { + const auto state = scenario.snapshot(); + std::cout << std::fixed << std::setprecision(1) + << "t=" << static_cast(i + 1) / sample_rate << " s rpm=" << state.rpm + << " cyl=" << state.cylinder_pressure / 100000.0f << " bar" + << " crankcase=" << state.crankcase_pressure / 1000.0f << " kPa" + << " intake=" << state.intake_pressure / 1000.0f << " kPa" + << " reed=" << (state.reed_open ? "open" : "closed") << '\n'; + } + } + } else if (scenario_name == "helmholtz") { + HelmholtzScenario scenario; + std::cout << "FluidSimCpp | Helmholtz resonator | initial cavity 1.1 bar\n"; + for (std::uint64_t i = 0; i < total_samples; ++i) { + const float sample = scenario.process(); + if ((i + 1) % report_interval == 0 || i + 1 == total_samples) + std::cout << std::fixed << std::setprecision(1) + << "t=" << static_cast(i + 1) / sample_rate << " s" + << " cavity=" << scenario.cavity_pressure() / 1000.0f << " kPa" + << " flow=" << scenario.open_end_flow() << " kg/s" + << " audio=" << sample << '\n'; + } + } else if (scenario_name == "pulse") { + PulseTestScenario scenario; + std::cout << "FluidSimCpp | pressure pulse in a 2 m pipe\n"; + for (std::uint64_t i = 0; i < total_samples; ++i) { + const float sample = scenario.process(); + if ((i + 1) % report_interval == 0 || i + 1 == total_samples) + std::cout << std::fixed << std::setprecision(1) + << "t=" << static_cast(i + 1) / sample_rate << " s" + << " left=" << scenario.pipe().state(0).pressure / 1000.0f << " kPa" + << " right=" << scenario.pipe().state(199).pressure / 1000.0f << " kPa" + << " flow=" << scenario.pipe().state(199).velocity << " m/s" + << " audio=" << sample << '\n'; + } + } else { + std::cerr << "Unknown scenario '" << scenario_name + << "'. Choose twin, two-stroke, helmholtz, or pulse.\n"; + return 2; + } + return 0; +} diff --git a/src/mechanical_sound_processor.cpp b/src/mechanical_sound_processor.cpp new file mode 100644 index 0000000..a4fcaee --- /dev/null +++ b/src/mechanical_sound_processor.cpp @@ -0,0 +1,132 @@ +#include "FluidSim/MechanicalSoundProcessor.hpp" + +#include +#include + +namespace fluidsim::port { +namespace { +constexpr float two_pi = 2.0f * pi; +constexpr float closed_threshold = 1e-7f; +} + +MechanicalSoundProcessor::Resonator::Resonator(float frequency_hz, float q, float sample_rate) { + const float omega = two_pi * frequency_hz / sample_rate; + const float radius = std::exp(-omega / (2.0f * q)); + b1_ = 2.0f * radius * std::cos(omega); + b2_ = -radius * radius; + gain_ = std::sin(omega); +} +float MechanicalSoundProcessor::Resonator::process(float input) { + const float output = input * gain_ + b1_ * y1_ + b2_ * y2_; + y2_ = y1_; + y1_ = output; + return output; +} + +MechanicalSoundProcessor::OnePoleLowPass::OnePoleLowPass(float cutoff_hz, float sample_rate) { + alpha_ = 1.0f - std::exp(-(two_pi * cutoff_hz / sample_rate)); +} +float MechanicalSoundProcessor::OnePoleLowPass::process(float input) { + state_ += alpha_ * (input - state_); + return state_; +} + +MechanicalSoundProcessor::MechanicalSoundProcessor(int sample_rate) + : dt_(1.0f / static_cast(std::max(sample_rate, 1))), + intake_tick_(4200.0f, 20.0f, static_cast(std::max(sample_rate, 1))), + exhaust_tick_(3400.0f, 20.0f, static_cast(std::max(sample_rate, 1))), + slap_filter_(450.0f, static_cast(std::max(sample_rate, 1))), + combustion_filter_(600.0f, static_cast(std::max(sample_rate, 1))), + chain_filter_(1500.0f, static_cast(std::max(sample_rate, 1))) { + slap_decay_ = std::exp(-dt_ / 0.018f); + combustion_decay_ = std::exp(-dt_ / 0.025f); + chain_decay_ = std::exp(-dt_ / 0.005f); +} + +float MechanicalSoundProcessor::process(const Crankshaft& crank, + const EngineCylinder& cylinder, float load) { + const float rpm = crank.omega * 60.0f / two_pi; + const float rpm_norm = std::clamp(std::abs(rpm) / 7000.0f, 0.0f, 1.5f); + const float normalized_load = std::clamp(load, 0.0f, 1.0f); + const float degrees = std::fmod((crank.angle + cylinder.phase) * 180.0f / pi, 720.0f); + const float intake_area = cylinder.valve_area(degrees, true); + const float exhaust_area = cylinder.valve_area(degrees, false); + const float pressure = cylinder.pressure; + const float pressure_rate = (pressure - previous_pressure_) / dt_; + previous_pressure_ = pressure; + + float cycle_position = (crank.angle + cylinder.phase) / (4.0f * pi); + cycle_position -= std::floor(cycle_position); + const float previous_cycle = previous_cycle_position_ < 0.0f ? cycle_position : previous_cycle_position_; + float crank_fraction = crank.angle / two_pi; + crank_fraction -= std::floor(crank_fraction); + const float previous_crank = previous_crank_fraction_ < 0.0f ? crank_fraction : previous_crank_fraction_; + + float sample = 0.0f; + float intake_impulse = 0.0f, exhaust_impulse = 0.0f; + const bool intake_opened = previous_intake_area_ <= closed_threshold && intake_area > closed_threshold; + const bool intake_closed = previous_intake_area_ > closed_threshold && intake_area <= closed_threshold; + const bool exhaust_opened = previous_exhaust_area_ <= closed_threshold && exhaust_area > closed_threshold; + const bool exhaust_closed = previous_exhaust_area_ > closed_threshold && exhaust_area <= closed_threshold; + const float open_amplitude = 0.20f + 0.30f * rpm_norm; + const float close_amplitude = 0.45f + 0.85f * rpm_norm; + if (intake_opened) intake_impulse += random_sample() * open_amplitude; + if (intake_closed) intake_impulse += random_sample() * close_amplitude * 1.15f; + if (exhaust_opened) exhaust_impulse += random_sample() * open_amplitude * 1.10f; + if (exhaust_closed) exhaust_impulse += random_sample() * close_amplitude * 1.30f; + sample += intake_tick_.process(intake_impulse) * valve_tick_gain; + sample += exhaust_tick_.process(exhaust_impulse) * valve_tick_gain; + + if (previous_cycle_position_ >= 0.0f) { + float amplitude = 0.0f; + if (crossed_fraction(previous_cycle, cycle_position, 0.0f)) + amplitude = (0.35f + 0.65f * rpm_norm) * (0.4f + 0.6f * normalized_load); + else if (crossed_fraction(previous_cycle, cycle_position, 0.5f)) + amplitude = (0.10f + 0.20f * rpm_norm) * 0.7f; + else if (crossed_fraction(previous_cycle, cycle_position, 0.25f) || + crossed_fraction(previous_cycle, cycle_position, 0.75f)) + amplitude = (0.08f + 0.18f * rpm_norm) * 0.8f; + slap_envelope_ = std::max(slap_envelope_, amplitude); + } + sample += slap_filter_.process(random_sample() * slap_envelope_) * piston_slap_gain; + slap_envelope_ *= slap_decay_; + + const bool power_stroke = cycle_position < 0.25f || cycle_position > 0.95f; + if (power_stroke && pressure_rate > 2e5f) { + const float drive = std::clamp(pressure_rate * 1.2e-7f, 0.0f, 3.0f); + combustion_envelope_ = std::max(combustion_envelope_, drive * (0.3f + 0.7f * normalized_load)); + } + sample += combustion_filter_.process(random_sample() * combustion_envelope_) * combustion_gain; + combustion_envelope_ *= combustion_decay_; + + if (previous_crank_fraction_ >= 0.0f && + (crossed_fraction(previous_crank, crank_fraction, 0.0f) || + crossed_fraction(previous_crank, crank_fraction, 0.5f))) + chain_envelope_ += 0.15f + 0.25f * rpm_norm; + sample += chain_filter_.process(random_sample() * chain_envelope_) * cam_chain_gain; + chain_envelope_ *= chain_decay_; + + const float crank_hz = std::abs(rpm) / 60.0f; + gear_phase1_ = std::fmod(gear_phase1_ + two_pi * crank_hz * 15.0f * dt_, two_pi); + gear_phase2_ = std::fmod(gear_phase2_ + two_pi * crank_hz * 30.2f * dt_, two_pi); + gear_phase3_ = std::fmod(gear_phase3_ + two_pi * crank_hz * 45.4f * dt_, two_pi); + const float gear_tone = std::sin(gear_phase1_) * 0.5f + + std::sin(gear_phase2_) * 0.3f + std::sin(gear_phase3_) * 0.2f; + sample += gear_tone * gear_whine_gain * rpm_norm; + + previous_intake_area_ = intake_area; + previous_exhaust_area_ = exhaust_area; + previous_cycle_position_ = cycle_position; + previous_crank_fraction_ = crank_fraction; + return std::tanh(sample * master_gain); +} + +float MechanicalSoundProcessor::random_sample() { return unit_(random_); } + +bool MechanicalSoundProcessor::crossed_fraction(float previous, float current, float target) { + if (current == previous) return false; + return current > previous ? previous < target && target <= current + : previous < target || target <= current; +} + +} // namespace fluidsim::port diff --git a/src/outdoor_exhaust_reverb.cpp b/src/outdoor_exhaust_reverb.cpp new file mode 100644 index 0000000..2674e51 --- /dev/null +++ b/src/outdoor_exhaust_reverb.cpp @@ -0,0 +1,127 @@ +#include "FluidSim/OutdoorExhaustReverb.hpp" +#include "FluidSim/Engine.hpp" + +#include +#include +#include + +namespace fluidsim::port { + +OutdoorExhaustReverb::DelayLine::DelayLine(std::size_t size) : buffer_(std::max(size, 1)) {} +float OutdoorExhaustReverb::DelayLine::read() const { return buffer_[position_]; } +float OutdoorExhaustReverb::DelayLine::read_write(float value) { + const float previous = buffer_[position_]; + buffer_[position_] = value; + position_ = (position_ + 1) % buffer_.size(); + return previous; +} +void OutdoorExhaustReverb::DelayLine::write(float value) { + buffer_[position_] = value; + position_ = (position_ + 1) % buffer_.size(); +} + +OutdoorExhaustReverb::OnePoleLowPass::OnePoleLowPass(float cutoff_hz, int sample_rate) { + const float angular = 2.0f * pi * cutoff_hz / static_cast(std::max(sample_rate, 1)); + alpha_ = 1.0f - std::exp(-angular); +} +float OutdoorExhaustReverb::OnePoleLowPass::process(float value) { + state_ += alpha_ * (value - state_); + return state_; +} + +OutdoorExhaustReverb::OutdoorExhaustReverb(int sample_rate) + : sample_rate_(std::max(sample_rate, 1)), dt_(1.0f / static_cast(std::max(sample_rate, 1))), + pre_delay_(delay_samples(sample_rate_, pre_delay_ms)), + ground_left_(delay_samples(sample_rate_, ground_delay_ms)), + ground_right_(delay_samples(sample_rate_, ground_delay_ms + 1.3f)), + wall_left_(delay_samples(sample_rate_, wall_delay_ms)), + wall_right_(delay_samples(sample_rate_, wall_delay_ms + 2.2f)), + ground_filter_left_(1400.0f, sample_rate_), ground_filter_right_(1400.0f, sample_rate_), + wall_filter_left_(3800.0f, sample_rate_), wall_filter_right_(3800.0f, sample_rate_) { + constexpr std::array lengths{3203, 4027, 5521, 7027, 8521, 10007, 11503, 13009}; + const int maximum_delay = std::max(4, sample_rate_ / 4); + float average_delay_seconds = 0.0f; + fdn_.reserve(lengths.size()); + low_filters_.reserve(lengths.size()); + mid_filters_.reserve(lengths.size()); + for (const int requested : lengths) { + const int size = std::min(requested, maximum_delay); + fdn_.emplace_back(static_cast(size)); + low_filters_.emplace_back(band_low_hz, sample_rate_); + mid_filters_.emplace_back(band_high_hz, sample_rate_); + average_delay_seconds += static_cast(size) / static_cast(sample_rate_); + } + average_delay_seconds /= static_cast(lengths.size()); + gain_low_ = std::pow(10.0f, -3.0f * average_delay_seconds / low_decay_seconds); + gain_mid_ = std::pow(10.0f, -3.0f * average_delay_seconds / mid_decay_seconds); + gain_high_ = std::pow(10.0f, -3.0f * average_delay_seconds / high_decay_seconds); + read_buffer_.resize(lengths.size()); + mix_buffer_.resize(lengths.size()); + + fast_attack_ = 1.0f - std::exp(-dt_ / 0.0015f); + fast_release_ = 1.0f - std::exp(-dt_ / 0.015f); + slow_attack_ = 1.0f - std::exp(-dt_ / 0.040f); + slow_release_ = 1.0f - std::exp(-dt_ / 0.300f); + bloom_decay_ = std::exp(-dt_ / (transient_decay_ms * 0.001f)); +} + +AudioFrame OutdoorExhaustReverb::process_stereo(float input) { + const float magnitude = std::abs(input); + fast_envelope_ += (magnitude > fast_envelope_ ? fast_attack_ : fast_release_) * (magnitude - fast_envelope_); + slow_envelope_ += (magnitude > slow_envelope_ ? slow_attack_ : slow_release_) * (magnitude - slow_envelope_); + const float bang = std::tanh(std::max(0.0f, fast_envelope_ - slow_envelope_) * 4.0f); + bloom_ = std::max(bloom_ * bloom_decay_, bang); + + const float pre = pre_delay_.read_write(input); + const float ground_left = ground_filter_left_.process(ground_left_.read_write(pre)); + const float ground_right = ground_filter_right_.process(ground_right_.read_write(pre)); + const float wall_left = wall_filter_left_.process(wall_left_.read_write(pre)); + const float wall_right = wall_filter_right_.process(wall_right_.read_write(pre)); + const float early_left = (ground_left * 0.85f + wall_left * 0.45f) * early_mix; + const float early_right = (ground_right * 0.85f + wall_right * 0.45f) * early_mix; + + float sum = 0.0f; + for (std::size_t i = 0; i < fdn_.size(); ++i) { + read_buffer_[i] = fdn_[i].read(); + sum += read_buffer_[i]; + } + const float mix_factor = 2.0f / static_cast(fdn_.size()); + for (std::size_t i = 0; i < fdn_.size(); ++i) + mix_buffer_[i] = mix_factor * sum - read_buffer_[i]; + + const float feedback_boost = 1.0f + transient_feedback_boost * bloom_; + const float hf_boost = 1.0f + transient_hf_boost * bloom_; + const float effective_low = std::min(gain_low_ * feedback_boost, 0.97f); + const float effective_mid = std::min(gain_mid_ * feedback_boost, 0.97f); + const float effective_high = std::min(gain_high_ * hf_boost, 0.90f); + const float injection = pre * 0.35355339f; + const float drive = std::max(saturation_drive, 1e-5f); + for (std::size_t i = 0; i < fdn_.size(); ++i) { + const float value = mix_buffer_[i]; + const float low = low_filters_[i].process(value); + const float low_mid = mid_filters_[i].process(value); + const float mid = low_mid - low; + const float high = value - low_mid; + const float feedback = low * effective_low + mid * effective_mid + high * effective_high; + fdn_[i].write(std::tanh(feedback * drive) / drive + injection); + } + + const float tail_left = (read_buffer_[0] + read_buffer_[2] + read_buffer_[4] + read_buffer_[6]) * 0.25f * tail_mix; + const float tail_right = (read_buffer_[1] + read_buffer_[3] + read_buffer_[5] + read_buffer_[7]) * 0.25f * tail_mix; + return {input * dry_mix + early_left + tail_left, + input * dry_mix + early_right + tail_right}; +} + +float OutdoorExhaustReverb::process(float input) { + const AudioFrame result = process_stereo(input); + const float value = 0.5f * (result.left + result.right); + const float magnitude = std::abs(value); + if (magnitude <= 1.0f) return value; + return std::copysign(1.0f + (1.0f - std::exp(-(magnitude - 1.0f))), value); +} + +std::size_t OutdoorExhaustReverb::delay_samples(int sample_rate, float milliseconds) { + return static_cast(std::max(1, static_cast(std::max(sample_rate, 1) * milliseconds * 0.001f))); +} + +} // namespace fluidsim::port diff --git a/src/pulse_test_scenario.cpp b/src/pulse_test_scenario.cpp new file mode 100644 index 0000000..3f146ab --- /dev/null +++ b/src/pulse_test_scenario.cpp @@ -0,0 +1,34 @@ +#include "FluidSim/PulseTestScenario.hpp" + +#include + +namespace fluidsim::port { + +PulseTestScenario::PulseTestScenario(int sample_rate) + : pipes_{}, dt_(1.0f / static_cast(std::max(sample_rate, 1))), + pipe_area_(pi * 0.25f * 0.020f * 0.020f), + boundaries_(pipes_), solver_(pipes_, boundaries_), sound_(sample_rate, 1.0f) { + pipes_.emplace_back(pipe_area_, 2.0f, 200); + auto& pipe = pipes_.front(); + pipe.damping_multiplier = 0.0f; + pipe.energy_relaxation_rate = 0.0f; + pipe.ambient_pressure = 101325.0f; + + const float bubble_pressure = 10.0f * 101325.0f; + const float bubble_temperature = 2000.0f; + const float bubble_density = bubble_pressure / (gas_constant * bubble_temperature); + for (int i = 0; i <= 10; ++i) + pipe.cells.at(static_cast(i)) = + {bubble_density, 0.0f, bubble_pressure / (gamma - 1.0f), 1.0f}; + + boundaries_.add_closed_end(0, true); + boundaries_.add_open_end(0, false, 101325.0f, pipe_area_); + solver_.minimum_substeps = 3; +} + +float PulseTestScenario::process() { + solver_.step(dt_); + return sound_.process(boundaries_.open_end_mass_flow(0)); +} + +} // namespace fluidsim::port diff --git a/src/scenario.cpp b/src/scenario.cpp new file mode 100644 index 0000000..396351d --- /dev/null +++ b/src/scenario.cpp @@ -0,0 +1,222 @@ +#include "FluidSim/Scenario.hpp" + +#include +#include + +namespace fluidsim::port { + +FourStrokeTwinScenario::FourStrokeTwinScenario(int sample_rate) + : pipes(create_pipes()), + dt_(1.0f / static_cast(std::max(sample_rate, 1))), + boundaries_(pipes), solver_(pipes, boundaries_), + intake_radiation_(sample_rate), exhaust_radiation_(sample_rate), + reverb_(sample_rate), mechanical_sound_(sample_rate) { + crankshaft.omega = 2000.0f * 2.0f * pi / 60.0f; + crankshaft.inertia = 0.04f; + crankshaft.friction_constant = 5.0f; + crankshaft.friction_viscous = 0.003f; + + vehicle.set_neutral(); + vehicle.clutch_input = clutch; + + cylinder1.intake_valve_diameter = 0.030f; + cylinder1.intake_valve_lift = 0.0080f; + cylinder1.exhaust_valve_diameter = 0.026f; + cylinder1.exhaust_valve_lift = 0.0075f; + cylinder2.phase = pi; + cylinder2.intake_valve_diameter = 0.030f; + cylinder2.intake_valve_lift = 0.0080f; + cylinder2.exhaust_valve_diameter = 0.026f; + cylinder2.exhaust_valve_lift = 0.0075f; + cylinder1.rev_limit_enabled = cylinder2.rev_limit_enabled = true; + cylinder1.rev_limit_rpm = cylinder2.rev_limit_rpm = 10000.0f; + cylinder1.rev_limit_cut_time = cylinder2.rev_limit_cut_time = 0.020f; + cylinder1.rev_limit_hysteresis = cylinder2.rev_limit_hysteresis = 250.0f; + cylinder1.initialize(crankshaft); + cylinder2.initialize(crankshaft); + + const float intake_area = area_for_diameter(0.0275f); + const float exhaust_area = area_for_diameter(0.028f); + runner_area_ = intake_area; + header_area_ = tail_area_ = exhaust_area; + for (auto& pipe : pipes) { + pipe.damping_multiplier = 1.0f; + pipe.energy_relaxation_rate = 0.5f; + pipe.ambient_pressure = 101325.0f; + } + + throttle_port_ = &intake_plenum.create_port(); + runner_feed1_ = &intake_plenum.create_port(); + runner_feed2_ = &intake_plenum.create_port(); + collector_in1_ = &exhaust_collector.create_port(); + collector_in2_ = &exhaust_collector.create_port(); + collector_out_ = &exhaust_collector.create_port(); + + intake_radiation_.pipe_end = {-0.20f, 0.0f, 0.0f}; + intake_radiation_.pipe_axis = {-1.0f, 0.0f, 0.0f}; + intake_radiation_.gain = 30.0f; + intake_radiation_.directivity_balance = 0.30f; + intake_radiation_.flow_tau = 0.005f; + intake_radiation_.smooth_tau = 0.020f; + exhaust_radiation_.pipe_end = {0.55f, 0.0f, 0.0f}; + exhaust_radiation_.pipe_axis = {1.0f, 0.0f, 0.0f}; + exhaust_radiation_.gain = 25.0f; + exhaust_radiation_.directivity_balance = 0.30f; + exhaust_radiation_.flow_tau = 0.005f; + exhaust_radiation_.smooth_tau = 0.020f; + + reverb_.dry_mix = 0.0f; + reverb_.early_mix = 0.30f; + reverb_.tail_mix = 0.22f; + reverb_.transient_feedback_boost = 0.40f; + reverb_.transient_hf_boost = 3.5f; + reverb_.saturation_drive = 1.5f; + + solver_.add_component(cylinder1); + solver_.add_component(cylinder2); + solver_.add_component(intake_plenum); + solver_.add_component(exhaust_collector); + + boundaries_.add_open_end(0, true, 101325.0f, intake_area); + boundaries_.add_open_end(5, false, 101325.0f, exhaust_area); + boundaries_.add_orifice(*throttle_port_, 0, false, throttle_area_, 0.7f); + boundaries_.add_orifice(*runner_feed1_, 1, true, runner_area_); + boundaries_.add_orifice(*runner_feed2_, 2, true, runner_area_); + boundaries_.add_orifice(cylinder1.intake, 1, false, valve1_in_); + boundaries_.add_orifice(cylinder2.intake, 2, false, valve2_in_); + boundaries_.add_orifice(cylinder1.exhaust, 3, true, valve1_out_); + boundaries_.add_orifice(cylinder2.exhaust, 4, true, valve2_out_); + boundaries_.add_orifice(*collector_in1_, 3, false, header_area_); + boundaries_.add_orifice(*collector_in2_, 4, false, header_area_); + boundaries_.add_orifice(*collector_out_, 5, true, tail_area_); +} + +AudioFrame FourStrokeTwinScenario::process() { + const std::lock_guard guard(mutex_); + crankshaft.load_torque = std::clamp(load, 0.0f, 1.0f) * 30.0f; + vehicle.clutch_input = std::clamp(clutch, 0.0f, 1.0f); + crankshaft.step(dt_); + crankshaft.drivetrain_torque = vehicle.update(crankshaft.omega, crankshaft.inertia, dt_); + cylinder1.pre_step(crankshaft, dt_, throttle); + cylinder2.pre_step(crankshaft, dt_, throttle); + + const float throttle_fraction = std::clamp(throttle, 0.0f, 1.0f); + throttle_area_ = throttle_fraction <= 0.001f ? 0.0f : area_for_diameter(0.026f) * throttle_fraction; + const auto cycle_degrees = [this](float phase) { + float value = std::fmod((crankshaft.angle + phase) * 180.0f / pi, 720.0f); + return value < 0.0f ? value + 720.0f : value; + }; + valve1_in_ = cylinder1.valve_area(cycle_degrees(cylinder1.phase), true); + valve2_in_ = cylinder2.valve_area(cycle_degrees(cylinder2.phase), true); + valve1_out_ = cylinder1.valve_area(cycle_degrees(cylinder1.phase), false); + valve2_out_ = cylinder2.valve_area(cycle_degrees(cylinder2.phase), false); + + solver_.step(dt_); + intake_flow = boundaries_.open_end_mass_flow(0); + exhaust_flow = boundaries_.open_end_mass_flow(1); + const AudioFrame intake_audio = intake_radiation_.process(intake_flow); + const AudioFrame exhaust_audio = exhaust_radiation_.process(exhaust_flow); + const float mechanical = mechanical_sound_enabled + ? mechanical_sound_.process(crankshaft, cylinder1, load) * mechanical_mix : 0.0f; + const float reverb_input = intake_radiation_.last_mono_pressure() * intake_send + + exhaust_radiation_.last_mono_pressure() * exhaust_send + + mechanical * mechanical_send; + const AudioFrame wet = reverb_.process_stereo(reverb_input); + const float flow_left = std::tanh(intake_audio.left + exhaust_audio.left + wet.left * reverb_send_gain); + const float flow_right = std::tanh(intake_audio.right + exhaust_audio.right + wet.right * reverb_send_gain); + last_mix_l = std::isfinite(flow_left + mechanical) + ? std::clamp(flow_left + mechanical, -1.0f, 1.0f) : 0.0f; + last_mix_r = std::isfinite(flow_right + mechanical) + ? std::clamp(flow_right + mechanical, -1.0f, 1.0f) : 0.0f; + ++sample_count; + return {last_mix_l, last_mix_r}; +} + +void FourStrokeTwinScenario::set_controls(float throttle_value, float load_value, float clutch_value) { + const std::lock_guard guard(mutex_); + throttle = std::clamp(throttle_value, 0.0f, 1.0f); + load = std::clamp(load_value, 0.0f, 1.0f); + clutch = std::clamp(clutch_value, 0.0f, 1.0f); +} + +void FourStrokeTwinScenario::toggle_mechanical_sound() { + const std::lock_guard guard(mutex_); + mechanical_sound_enabled = !mechanical_sound_enabled; +} + +void FourStrokeTwinScenario::shift_up() { + const std::lock_guard guard(mutex_); + vehicle.shift_up(); +} + +void FourStrokeTwinScenario::shift_down() { + const std::lock_guard guard(mutex_); + vehicle.shift_down(); +} + +void FourStrokeTwinScenario::set_neutral() { + const std::lock_guard guard(mutex_); + vehicle.set_neutral(); +} + +void FourStrokeTwinScenario::set_first_gear() { + const std::lock_guard guard(mutex_); + vehicle.set_first_gear(); +} + +ScenarioSnapshot FourStrokeTwinScenario::snapshot() const { + const std::lock_guard guard(mutex_); + ScenarioSnapshot result; + result.rpm = crankshaft.omega * 60.0f / (2.0f * pi); + result.speed_kmh = vehicle.speed * 3.6f; + result.kinetic_energy = vehicle.kinetic_energy(); + result.average_power = crankshaft.average_power; + result.average_torque = crankshaft.average_torque; + result.cylinder_pressure1 = cylinder1.pressure; + result.cylinder_pressure2 = cylinder2.pressure; + result.cylinder_temperature1 = cylinder1.temperature(); + result.cylinder_temperature2 = cylinder2.temperature(); + result.cylinder_piston_fraction1 = (cylinder1.volume - cylinder1.clearance_volume()) / cylinder1.swept_volume(); + result.cylinder_piston_fraction2 = (cylinder2.volume - cylinder2.clearance_volume()) / cylinder2.swept_volume(); + result.valve1_intake_area = valve1_in_; + result.valve1_exhaust_area = valve1_out_; + result.valve2_intake_area = valve2_in_; + result.valve2_exhaust_area = valve2_out_; + result.intake_pressure = intake_plenum.pressure(); + result.exhaust_pressure = exhaust_collector.pressure(); + result.throttle = throttle; + result.clutch = clutch; + result.intake_flow = intake_flow; + result.exhaust_flow = exhaust_flow; + result.gear = vehicle.current_gear; + result.rev_limit1 = cylinder1.rev_limit_active; + result.rev_limit2 = cylinder2.rev_limit_active; + result.mechanical_sound_enabled = mechanical_sound_enabled; + for (std::size_t pipe_index = 0; pipe_index < pipes.size(); ++pipe_index) { + auto& states = result.pipe_states[pipe_index]; + states.reserve(pipes[pipe_index].cells.size()); + for (std::size_t cell = 0; cell < pipes[pipe_index].cells.size(); ++cell) + states.push_back(pipes[pipe_index].state(static_cast(cell))); + } + return result; +} + +float FourStrokeTwinScenario::area_for_diameter(float diameter) { + return pi * diameter * diameter * 0.25f; +} + +std::vector FourStrokeTwinScenario::create_pipes() { + const float intake_area = area_for_diameter(0.0275f); + const float exhaust_area = area_for_diameter(0.028f); + std::vector result; + result.reserve(6); + result.emplace_back(intake_area, 0.06f, 10); + result.emplace_back(intake_area, 0.12f, 10); + result.emplace_back(intake_area, 0.12f, 10); + result.emplace_back(exhaust_area, 0.20f, 6); + result.emplace_back(exhaust_area, 0.20f, 6); + result.emplace_back(exhaust_area, 0.40f, 10); + return result; +} + +} // namespace fluidsim::port diff --git a/src/simulation_ring_buffer.cpp b/src/simulation_ring_buffer.cpp new file mode 100644 index 0000000..d46f66c --- /dev/null +++ b/src/simulation_ring_buffer.cpp @@ -0,0 +1,84 @@ +#include "FluidSim/SimulationRingBuffer.hpp" + +#include +#include +#include + +namespace fluidsim::port { + +SimulationRingBuffer::SimulationRingBuffer(std::size_t capacity) + : left_(capacity), right_(capacity), mask_(capacity - 1) { + if (capacity < 2 || (capacity & (capacity - 1)) != 0) + throw std::invalid_argument("simulation audio ring capacity must be a power of two"); +} + +std::size_t SimulationRingBuffer::available_samples() const { + const auto written = write_head_.load(std::memory_order_acquire); + const auto read = read_head_.load(std::memory_order_acquire); + return static_cast(written - read); +} + +std::size_t SimulationRingBuffer::free_space() const { + return left_.size() - available_samples(); +} + +bool SimulationRingBuffer::write(float left, float right, std::stop_token stop) { + while (free_space() == 0) { + if (stop.stop_requested()) return false; + std::unique_lock lock(wait_mutex_); + space_available_.wait_for(lock, std::chrono::milliseconds(10)); + } + if (stop.stop_requested()) return false; + const auto write = write_head_.load(std::memory_order_relaxed); + left_[static_cast(write) & mask_] = left; + right_[static_cast(write) & mask_] = right; + write_head_.store(write + 1, std::memory_order_release); + return true; +} + +bool SimulationRingBuffer::wait_for_space(std::size_t target_fill, std::stop_token stop) { + while (available_samples() > target_fill) { + if (stop.stop_requested()) return false; + std::unique_lock lock(wait_mutex_); + space_available_.wait_for(lock, std::chrono::milliseconds(10)); + } + return !stop.stop_requested(); +} + +void SimulationRingBuffer::reset_consumer() { + consumer_initialized_ = false; +} + +std::size_t SimulationRingBuffer::read_interpolated(std::span dst_left, + std::span dst_right, + double speed) { + const std::size_t count = std::min(dst_left.size(), dst_right.size()); + if (!consumer_initialized_) { + read_position_ = static_cast(read_head_.load(std::memory_order_acquire)); + consumer_initialized_ = true; + } + const auto write = write_head_.load(std::memory_order_acquire); + std::size_t produced = 0; + for (; produced < count; ++produced) { + const auto floor_position = static_cast(read_position_); + const auto ceil_position = floor_position + 1; + if (ceil_position >= write) break; + const float left0 = left_[static_cast(floor_position) & mask_]; + const float left1 = left_[static_cast(ceil_position) & mask_]; + const float right0 = right_[static_cast(floor_position) & mask_]; + const float right1 = right_[static_cast(ceil_position) & mask_]; + const float fraction = static_cast(read_position_ - static_cast(floor_position)); + dst_left[produced] = left0 + (left1 - left0) * fraction; + dst_right[produced] = right0 + (right1 - right0) * fraction; + read_position_ += std::clamp(speed, 0.0001, 1.0); + } + const auto new_read_head = static_cast(read_position_); + const auto previous_read_head = read_head_.load(std::memory_order_relaxed); + if (new_read_head > previous_read_head) { + read_head_.store(new_read_head, std::memory_order_release); + space_available_.notify_one(); + } + return produced; +} + +} // namespace fluidsim::port diff --git a/src/solver.cpp b/src/solver.cpp new file mode 100644 index 0000000..92a4b87 --- /dev/null +++ b/src/solver.cpp @@ -0,0 +1,27 @@ +#include "FluidSim/Solver.hpp" + +#include + +namespace fluidsim::port { + +Solver::Solver(std::vector& pipes, BoundarySystem& boundaries) + : pipes_(pipes), boundaries_(boundaries) {} + +void Solver::add_component(Component& component) { + components_.push_back(&component); +} + +void Solver::step(float dt) { + if (dt <= 0.0f) return; + int substeps = std::max(minimum_substeps, 1); + for (const auto& pipe : pipes_) substeps = std::max(substeps, pipe.required_substeps(dt, 0.8f)); + const float sub_dt = dt / static_cast(substeps); + for (int i = 0; i < substeps; ++i) { + boundaries_.resolve(sub_dt); + for (auto& pipe : pipes_) pipe.advance(sub_dt); + } + for (Component* component : components_) component->update_state(dt); + ++step_count_; +} + +} // namespace fluidsim::port diff --git a/src/sound_engine.cpp b/src/sound_engine.cpp new file mode 100644 index 0000000..daa2824 --- /dev/null +++ b/src/sound_engine.cpp @@ -0,0 +1,39 @@ +#include "FluidSim/SoundEngine.hpp" + +namespace fluidsim::port { + +SoundEngine::SoundEngine(SimulationRingBuffer& source_buffer) : stream_(source_buffer) {} + +SoundEngine::~SoundEngine() { + stop(); +} + +void SoundEngine::start() { + if (is_playing_) return; + stream_.play(); + is_playing_ = true; +} + +void SoundEngine::stop() { + if (!is_playing_) return; + stream_.stop(); + is_playing_ = false; +} + +void SoundEngine::set_speed(double speed_value) { + stream_.set_speed(speed_value); +} + +double SoundEngine::speed() const { + return stream_.speed(); +} + +void SoundEngine::set_volume(float volume_value) { + stream_.setVolume(volume_value); +} + +float SoundEngine::volume() const { + return stream_.getVolume(); +} + +} // namespace fluidsim::port diff --git a/src/two_stroke_cylinder.cpp b/src/two_stroke_cylinder.cpp new file mode 100644 index 0000000..2313d2c --- /dev/null +++ b/src/two_stroke_cylinder.cpp @@ -0,0 +1,123 @@ +#include "FluidSim/TwoStrokeCylinder.hpp" + +#include +#include +#include + +namespace fluidsim::port { + +TwoStrokeCylinder::TwoStrokeCylinder(float bore_value, float stroke_value, + float connecting_rod_length, float compression_ratio, + float transfer_duration_degrees, float exhaust_duration_degrees) + : transfer_duration_(transfer_duration_degrees), exhaust_duration_(exhaust_duration_degrees) { + if (transfer_duration_ <= 0.0f || transfer_duration_ >= 360.0f || + exhaust_duration_ <= 0.0f || exhaust_duration_ >= 360.0f) + throw std::invalid_argument("two-stroke port durations must be between 0 and 360 degrees"); + bore = bore_value; + stroke = stroke_value; + rod = connecting_rod_length; + compression = compression_ratio; + if (exhaust_open() >= transfer_open()) + throw std::invalid_argument("two-stroke exhaust must open before the transfer port"); +} + +float TwoStrokeCylinder::lift(float angle, float opens, float closes, float peak_lift) { + float mapped = std::fmod(angle, 360.0f); + if (mapped < 0.0f) mapped += 360.0f; + float effective_close = closes; + if (effective_close < opens) effective_close += 360.0f; + if (mapped < opens) mapped += 360.0f; + const float duration = effective_close - opens; + if (mapped < opens || mapped > effective_close || duration <= 0.0f) return 0.0f; + + const float ramp = duration * 0.25f; + const float hold = duration - 2.0f * ramp; + if (mapped < opens + ramp) { + const float t = (mapped - opens) / ramp; + return peak_lift * t * t * (3.0f - 2.0f * t); + } + if (mapped < opens + ramp + hold) return peak_lift; + const float t = (mapped - (opens + ramp + hold)) / ramp; + return peak_lift * (1.0f - t) * (1.0f - t) * (1.0f + 2.0f * t); +} + +float TwoStrokeCylinder::valve_area(float cycle_degrees, bool inlet) const { + const float opens = inlet ? transfer_open() : exhaust_open(); + const float closes = inlet ? transfer_close() : exhaust_close(); + const float diameter = inlet ? intake_valve_diameter : exhaust_valve_diameter; + const float peak_lift = inlet ? intake_valve_lift : exhaust_valve_lift; + return pi * diameter * lift(cycle_degrees, opens, closes, peak_lift); +} + +void TwoStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle*/) { + const float rpm = crank.omega * 60.0f / (2.0f * pi); + spark_advance = std::clamp(10.0f + rpm * 0.002f, 5.0f, 40.0f); + update_rev_limiter(rpm, dt); + + const float theta = crank.angle + phase; + const float previous_volume = volume; + volume = compute_volume(theta); + pressure = (gamma - 1.0f) * energy / std::max(volume, 1e-12f); + const float dvolume = volume - previous_volume; + const float radius = stroke * 0.5f; + const float obliquity = radius / rod; + const float sine = std::sin(theta), cosine = std::cos(theta); + const float term = std::sqrt(std::max(1.0f - obliquity * obliquity * sine * sine, 1e-12f)); + const float dx_dtheta = radius * sine * (1.0f + obliquity * cosine / term); + const float piston_area = pi * 0.25f * bore * bore; + const float back_pressure = crankcase_ ? crankcase_->pressure() : 101325.0f; + crank.add_torque((pressure - back_pressure) * piston_area * dx_dtheta); + energy -= pressure * dvolume; + + const float previous_degrees = wrap_degrees((crank.previous_angle + phase) * 180.0f / pi, 360.0f); + const float current_degrees = wrap_degrees(theta * 180.0f / pi, 360.0f); + handle_cycle_events(previous_degrees, current_degrees); + + if (rev_limit_active) { + fuel_mass = 0.0f; + fuel_injected = false; + combustion_active = false; + } + energy -= heat_transfer * wall_area * (temperature() - ambient_temperature) * dt; + sync_ports(); +} + +void TwoStrokeCylinder::handle_cycle_events(float previous_degrees, float current_degrees) { + if (previous_degrees >= transfer_open() && previous_degrees < transfer_close() && + current_degrees >= transfer_close()) { + fuel_mass = air_mass / std::max(afr, 1e-6f); + fuel_injected = true; + } + + const float spark_angle = wrap_degrees(360.0f - spark_advance, 360.0f); + const bool crossed_spark = (previous_degrees < spark_angle && current_degrees >= spark_angle) || + (previous_degrees > spark_angle && current_degrees < spark_angle); + if (crossed_spark && fuel_injected && !combustion_active && !rev_limit_active && + random_unit() >= misfire_probability) { + combustion_active = true; + burn_fraction = 0.0f; + energy_factor_ = 1.0f + energy_variation * (2.0f * random_unit() - 1.0f); + } + + if (!combustion_active) return; + float angle_since_spark = current_degrees - spark_angle; + if (angle_since_spark < 0.0f) angle_since_spark += 360.0f; + const float phi = (angle_since_spark - burn_start) / std::max(burn_duration, 1e-6f); + float fraction = phi <= 0.0f ? 0.0f : 1.0f - std::exp(-wiebe_a * std::pow(phi, wiebe_m + 1.0f)); + fraction = std::max(burn_fraction, std::clamp(fraction, 0.0f, 1.0f)); + if (fraction >= 1.0f || angle_since_spark > burn_duration + burn_start + spark_advance) { + fraction = 1.0f; + combustion_active = false; + exhaust_mass += air_mass; + air_mass = 0.0f; + } + fuel_injected = false; + const float increment = fraction - burn_fraction; + if (increment > 0.0f) { + energy += fuel_mass * fuel_lhv * energy_factor_ * increment; + exhaust_mass += fuel_mass * increment; + burn_fraction = fraction; + } +} + +} // namespace fluidsim::port diff --git a/src/two_stroke_scenario.cpp b/src/two_stroke_scenario.cpp new file mode 100644 index 0000000..b685f9b --- /dev/null +++ b/src/two_stroke_scenario.cpp @@ -0,0 +1,146 @@ +#include "FluidSim/TwoStrokeScenario.hpp" + +#include +#include + +namespace fluidsim::port { + +TwoStrokeScenario::TwoStrokeScenario(int sample_rate) + : cylinder(0.054f, 0.0545f, 0.110f, 7.2f, 140.0f, 195.0f), + crankcase(crankshaft, 0.0545f * 0.5f, 0.110f, 0.054f, 150e-6f), + pipes(create_pipes()), dt_(1.0f / static_cast(std::max(sample_rate, 1))), + max_throttle_area_(area_for_diameter(0.042f)), + reed_in_area_(area_for_diameter(0.040f)), + reed_out_area_(0.0f), transfer_pipe_area_(area_for_diameter(0.040f)), + transfer_valve_area_(0.0f), exhaust_valve_area_(0.0f), boundaries_(pipes), solver_(pipes, boundaries_), + intake_sound_(sample_rate), exhaust_sound_(sample_rate), reverb_(sample_rate) { + crankshaft.omega = 3000.0f * 2.0f * pi / 60.0f; + crankshaft.inertia = 0.01f; + crankshaft.friction_constant = 1.0f; + crankshaft.friction_viscous = 0.002f; + + cylinder.intake_valve_diameter = 0.030f; + cylinder.intake_valve_lift = 0.010f; + cylinder.exhaust_valve_diameter = 0.040f; + cylinder.exhaust_valve_lift = 0.010f; + cylinder.initialize(crankshaft); + cylinder.set_crankcase(crankcase); + + for (auto& pipe : pipes) { + pipe.damping_multiplier = 0.8f; + pipe.energy_relaxation_rate = 0.4f; + pipe.ambient_pressure = 101325.0f; + } + plenum_inlet_ = &intake_plenum.create_port(); + plenum_outlet_ = &intake_plenum.create_port(); + + boundaries_.add_open_end(0, true, 101325.0f, pipes[0].area); + boundaries_.add_open_end(3, false, 101325.0f, pipes[3].cell_area(static_cast(pipes[3].cells.size()) - 1)); + boundaries_.add_orifice(*plenum_inlet_, 0, false, max_throttle_area_, 0.72f); + boundaries_.add_orifice(*plenum_outlet_, 1, true, reed_in_area_); + boundaries_.add_orifice(crankcase.intake_port, 1, false, reed_out_area_, 0.9f); + boundaries_.add_orifice(crankcase.transfer_port, 2, true, transfer_pipe_area_); + boundaries_.add_orifice(cylinder.intake, 2, false, transfer_valve_area_); + boundaries_.add_orifice(cylinder.exhaust, 3, true, exhaust_valve_area_); + solver_.minimum_substeps = 4; + solver_.add_component(intake_plenum); + solver_.add_component(crankcase); + solver_.add_component(cylinder); + + intake_sound_.gain = 4.5f; + exhaust_sound_.gain = 4.5f; +} + +AudioFrame TwoStrokeScenario::process() { + const float reed_margin = 200.0f; + if (crankcase.pressure() < intake_plenum.pressure() - reed_margin) reed_open_ = true; + else if (crankcase.pressure() > intake_plenum.pressure() + reed_margin) reed_open_ = false; + + const float throttle_fraction = std::clamp(throttle, 0.0f, 1.0f); + max_throttle_area_ = throttle_fraction < 0.001f ? 0.0f : area_for_diameter(0.042f) * throttle_fraction; + reed_out_area_ = reed_open_ ? reed_in_area_ : 0.0f; + if (step_count_ < 20000) crankshaft.add_torque(5.0f); + crankshaft.step(dt_); + crankcase.pre_step(dt_); + cylinder.pre_step(crankshaft, dt_, throttle_fraction); + transfer_valve_area_ = cylinder.valve_area(crankshaft.angle * 180.0f / pi, true); + exhaust_valve_area_ = cylinder.valve_area(crankshaft.angle * 180.0f / pi, false); + solver_.step(dt_); + ++step_count_; + + const float intake = intake_sound_.process(boundaries_.open_end_mass_flow(0)); + const float exhaust = exhaust_sound_.process(boundaries_.open_end_mass_flow(1)); + const float reverb_input = 0.5f * (intake + exhaust); + return reverb_.process_stereo(reverb_input); +} + +void TwoStrokeScenario::set_throttle(float value) { + throttle = std::clamp(value, 0.0f, 1.0f); +} + +TwoStrokeSnapshot TwoStrokeScenario::snapshot() const { + return { + crankshaft.omega * 60.0f / (2.0f * pi), cylinder.pressure, + crankcase.pressure(), intake_plenum.pressure(), throttle, + cylinder.valve_area(crankshaft.angle * 180.0f / pi, true), + cylinder.valve_area(crankshaft.angle * 180.0f / pi, false), reed_open_ + }; +} + +float TwoStrokeScenario::area_for_diameter(float diameter) { + return pi * diameter * diameter * 0.25f; +} + +std::vector TwoStrokeScenario::create_pipes() { + const float intake_area = area_for_diameter(0.042f); + const float reed_area = area_for_diameter(0.040f); + const float transfer_area = area_for_diameter(0.040f); + std::vector result; + result.reserve(4); + result.emplace_back(intake_area, 0.100f, 8); + result.emplace_back(reed_area, 0.030f, 4); + result.emplace_back(transfer_area, 0.200f, 8); + + constexpr int exhaust_cells = 60; + constexpr float header_length = 0.130f, diffuser_length = 0.250f; + constexpr float belly_length = 0.220f, convergent_length = 0.160f; + constexpr float stinger_length = 0.080f; + constexpr float header_diameter = 0.040f, belly_diameter = 0.070f; + constexpr float stinger_diameter = 0.028f; + constexpr float total_length = header_length + diffuser_length + belly_length + + convergent_length + stinger_length; + const int header_cells = static_cast(exhaust_cells * header_length / total_length); + const int diffuser_cells = static_cast(exhaust_cells * diffuser_length / total_length); + const int belly_cells = static_cast(exhaust_cells * belly_length / total_length); + const int convergent_cells = static_cast(exhaust_cells * convergent_length / total_length); + const int stinger_cells = exhaust_cells - header_cells - diffuser_cells - belly_cells - convergent_cells; + std::vector exhaust_areas; + std::vector exhaust_lengths; + exhaust_areas.reserve(exhaust_cells); + exhaust_lengths.reserve(exhaust_cells); + const auto append_uniform = [&](int count, float diameter, float length) { + for (int i = 0; i < count; ++i) { + exhaust_areas.push_back(area_for_diameter(diameter)); + exhaust_lengths.push_back(length / static_cast(count)); + } + }; + append_uniform(header_cells, header_diameter, header_length); + for (int i = 0; i < diffuser_cells; ++i) { + const float t = diffuser_cells > 1 ? static_cast(i) / (diffuser_cells - 1) : 0.0f; + const float diameter = header_diameter + (belly_diameter - header_diameter) * t; + exhaust_areas.push_back(area_for_diameter(diameter)); + exhaust_lengths.push_back(diffuser_length / static_cast(diffuser_cells)); + } + append_uniform(belly_cells, belly_diameter, belly_length); + for (int i = 0; i < convergent_cells; ++i) { + const float t = convergent_cells > 1 ? static_cast(i) / (convergent_cells - 1) : 0.0f; + const float diameter = belly_diameter + (stinger_diameter - belly_diameter) * t; + exhaust_areas.push_back(area_for_diameter(diameter)); + exhaust_lengths.push_back(convergent_length / static_cast(convergent_cells)); + } + append_uniform(stinger_cells, stinger_diameter, stinger_length); + result.emplace_back(std::move(exhaust_areas), std::move(exhaust_lengths)); + return result; +} + +} // namespace fluidsim::port diff --git a/src/vehicle.cpp b/src/vehicle.cpp new file mode 100644 index 0000000..85f2a81 --- /dev/null +++ b/src/vehicle.cpp @@ -0,0 +1,61 @@ +#include "FluidSim/Vehicle.hpp" + +#include +#include + +namespace fluidsim::port { + +float Vehicle::total_gear_ratio() const { + return current_gear == 0 ? 0.0f : gear_ratios.at(static_cast(current_gear - 1)) * final_drive_ratio * primary_reduction; +} + +float Vehicle::kinetic_energy() const { + return 0.5f * (mass + rotational_inertia / (wheel_radius * wheel_radius)) * speed * speed; +} + +float Vehicle::update(float engine_omega, float engine_inertia, float dt) { + if (dt <= 0.0f || engine_inertia <= 0.0f || wheel_radius <= 0.0f || mass <= 0.0f) return 0.0f; + const float target = current_gear == 0 || shift_cut_timer_ > 0.0f + ? 0.0f : std::clamp(1.0f - clutch_input, 0.0f, 1.0f); + engagement += (target - engagement) * (1.0f - std::exp(-clutch_engagement_rate * dt)); + shift_cut_timer_ = std::max(0.0f, shift_cut_timer_ - dt); + const float ratio = total_gear_ratio(); + if (ratio <= 0.0f || engagement <= 0.0f) { + apply_road_resistance(dt); + return 0.0f; + } + + const float vehicle_mass = mass + rotational_inertia / (wheel_radius * wheel_radius); + const float wheel_inertia = vehicle_mass * wheel_radius * wheel_radius; + const float reflected_engine = engine_inertia * ratio * ratio; + const float wheel_omega = speed / wheel_radius; + const float slip = engine_omega / ratio - wheel_omega; + const float sync_impulse = slip * (wheel_inertia * reflected_engine / (wheel_inertia + reflected_engine)); + const float max_impulse = clutch_capacity * engagement * dt * ratio; + float wheel_impulse = std::clamp(sync_impulse, -max_impulse, max_impulse); + wheel_impulse = std::max(wheel_impulse, -wheel_omega * wheel_inertia); + speed += wheel_impulse / (vehicle_mass * wheel_radius); + apply_road_resistance(dt); + return -wheel_impulse / (ratio * dt); +} + +void Vehicle::shift_up() { + if (current_gear < static_cast(gear_ratios.size())) { ++current_gear; shift_cut_timer_ = 0.15f; } +} +void Vehicle::shift_down() { + if (current_gear > 1) { --current_gear; shift_cut_timer_ = 0.15f; } +} +void Vehicle::set_neutral() { current_gear = 0; shift_cut_timer_ = 0.0f; } +void Vehicle::set_first_gear() { + if (current_gear == 0) { current_gear = 1; shift_cut_timer_ = 0.15f; } +} + +void Vehicle::apply_road_resistance(float dt) { + const float v = std::max(speed, 0.0f); + const float rolling = rolling_resistance * mass * gravity; + const float aerodynamic = 0.5f * air_density * drag_coefficient * frontal_area * v * v; + const float effective_mass = mass + rotational_inertia / (wheel_radius * wheel_radius); + speed = std::max(0.0f, speed - (rolling + aerodynamic) / std::max(effective_mass, 1e-6f) * dt); +} + +} // namespace fluidsim::port diff --git a/tests/physics_tests.cpp b/tests/physics_tests.cpp new file mode 100644 index 0000000..da7699f --- /dev/null +++ b/tests/physics_tests.cpp @@ -0,0 +1,205 @@ +#include "FluidSim/Scenario.hpp" +#include "FluidSim/Atmosphere.hpp" +#include "FluidSim/Crankcase.hpp" +#include "FluidSim/DynoCurve.hpp" +#include "FluidSim/TwoStrokeCylinder.hpp" +#include "FluidSim/TwoStrokeScenario.hpp" +#include "FluidSim/OutdoorExhaustReverb.hpp" +#include "FluidSim/HelmholtzScenario.hpp" +#include "FluidSim/PulseTestScenario.hpp" +#include "FluidSim/ThreadLoadTracker.hpp" +#include "FluidSim/Units.hpp" + +#include +#include +#include + +namespace { +void require(bool condition, const char* message) { + if (!condition) { + std::cerr << "FAIL: " << message << '\n'; + throw std::runtime_error(message); + } +} +} + +int main() { + using namespace fluidsim::port; + + require(std::abs(units::celsius(25.0) - 298.15) < 1e-9 && + std::abs(units::area_from_diameter(0.02) - static_cast(pi) * 1e-4) < 1e-12, + "unit helpers must preserve the C# conversions"); + DynoCurve dyno_curve; + dyno_curve.update(1000.0f, 10.0f, 50.0f); + dyno_curve.update(1150.0f, 22.0f, 80.0f); + require(dyno_curve.bins().size() == 24 && + std::abs(dyno_curve.bins()[21].power_kw - 14.0f) < 1e-5f && + std::abs(dyno_curve.bins()[22].torque_nm - 70.0f) < 1e-5f, + "dyno history must interpolate skipped 50-RPM bins like the C# scenario base"); + dyno_curve.update(1150.0f, 5.0f, 40.0f); + require(std::abs(dyno_curve.bins()[23].power_kw - 22.0f) < 1e-5f && + std::abs(dyno_curve.bins()[23].torque_nm - 80.0f) < 1e-5f, + "dyno bins must retain their peak power and torque"); + dyno_curve.reset(); + require(dyno_curve.bins().empty(), "dyno reset must clear accumulated history"); + ThreadLoadTracker thread_load; + thread_load.record(1.0, 2.0); + require(std::abs(thread_load.load_percent() - 2.5) < 1e-9 && + std::abs(thread_load.last_step_microseconds() - 1000.0) < 1e-9, + "thread load tracker must report EMA budget use and last-step time"); + + Crankshaft crank; + FourStrokeCylinder cylinder; + cylinder.initialize(crank); + require(std::abs(cylinder.compute_volume(0.0f) - cylinder.clearance_volume()) < 1e-8f, + "cylinder volume at TDC must equal clearance volume"); + require(std::abs(cylinder.compute_volume(pi) - cylinder.clearance_volume() - cylinder.swept_volume()) < 1e-8f, + "cylinder volume at BDC must include the swept volume"); + + TwoStrokeCylinder two_stroke{0.054f, 0.0545f, 0.110f, 7.2f, 140.0f, 195.0f}; + two_stroke.initialize(crank); + two_stroke.intake_valve_diameter = 0.030f; + two_stroke.intake_valve_lift = 0.010f; + two_stroke.exhaust_valve_diameter = 0.040f; + two_stroke.exhaust_valve_lift = 0.010f; + require(two_stroke.valve_area(180.0f, true) > 0.0f && + two_stroke.valve_area(0.0f, true) == 0.0f, + "two-stroke transfer port must open around BDC only"); + require(two_stroke.valve_area(180.0f, false) > 0.0f && + two_stroke.valve_area(0.0f, false) == 0.0f, + "two-stroke exhaust port must open around BDC only"); + + TwoStrokeScenario two_stroke_scenario; + two_stroke_scenario.set_throttle(0.25f); + for (int i = 0; i < 2000; ++i) two_stroke_scenario.process(); + const auto two_stroke_state = two_stroke_scenario.snapshot(); + require(std::isfinite(two_stroke_state.rpm) && std::isfinite(two_stroke_state.cylinder_pressure) && + std::isfinite(two_stroke_state.crankcase_pressure) && std::isfinite(two_stroke_state.intake_pressure), + "two-stroke scenario must remain finite while running"); + + Crankcase crankcase{crank, 0.02725f, 0.110f, 0.054f, 150e-6f}; + two_stroke.set_crankcase(crankcase); + crankcase.pre_step(1.0f / 44100.0f); + require(std::isfinite(crankcase.pressure()) && crankcase.pressure() > 0.0f && + crankcase.volume() > 0.0f, + "crankcase must retain a finite thermodynamic state"); + + Atmosphere atmosphere; + require(std::abs(atmosphere.port.pressure - 101325.0f) < 1e-3f && + std::abs(atmosphere.port.air_fraction - 1.0f) < 1e-6f, + "atmosphere must expose ambient air through its port"); + + OutdoorExhaustReverb reverb{44100}; + bool reverb_tail = false; + for (int i = 0; i < 12000; ++i) { + const AudioFrame output = reverb.process_stereo(i == 0 ? 0.8f : 0.0f); + require(std::isfinite(output.left) && std::isfinite(output.right), + "outdoor reverb must remain finite after a transient"); + if (i > 300 && (std::abs(output.left) > 1e-5f || std::abs(output.right) > 1e-5f)) reverb_tail = true; + } + require(reverb_tail, "outdoor reverb must produce a delayed reflection tail"); + SoundProcessor flow_sound{44100, 1.0f}; + for (int i = 0; i < 2000; ++i) + require(std::isfinite(flow_sound.process(i < 1000 ? 0.0f : 0.01f)), + "flow sound processor must remain finite for a flow step"); + + HelmholtzScenario helmholtz; + float resonator_peak = 0.0f; + for (int i = 0; i < 3000; ++i) { + const float sample = helmholtz.process(); + require(std::isfinite(sample) && std::isfinite(helmholtz.cavity_pressure()), + "Helmholtz resonator must remain finite while discharging"); + resonator_peak = std::max(resonator_peak, std::abs(sample)); + } + require(resonator_peak > 0.0f, "overpressurized Helmholtz cavity must excite the neck"); + + PulseTestScenario pulse_test; + for (int i = 0; i < 3000; ++i) { + const float sample = pulse_test.process(); + require(std::isfinite(sample) && std::isfinite(pulse_test.pipe().state(0).pressure), + "pulse test pipe must remain finite while the pressure pulse travels"); + } + + Volume0D chamber{1e-3f}; + Port& chamber_port = chamber.create_port(); + std::vector pipes; + pipes.emplace_back(1e-3f, 0.1f, 4); + BoundarySystem boundary{pipes}; + float closed_area = 0.0f; + boundary.add_orifice(chamber_port, 0, false, closed_area); + const float initial_pressure = chamber.pressure(); + boundary.resolve(1e-5f); + require(chamber_port.mass_flow == 0.0f, "closed orifice must have zero mass flow"); + require(std::abs(chamber.pressure() - initial_pressure) < 1e-2f, + "resolving a closed orifice must not alter the chamber state"); + + Volume0D equal_pressure_volume{1e-3f}; + Port& equal_pressure_port = equal_pressure_volume.create_port(); + equal_pressure_port.density = 1.0f; + std::vector equal_pressure_pipes; + equal_pressure_pipes.emplace_back(1e-3f, 0.1f, 4); + BoundarySystem equal_pressure_boundary{equal_pressure_pipes}; + float equal_pressure_area = 1e-4f; + equal_pressure_boundary.add_orifice(equal_pressure_port, 0, false, equal_pressure_area); + equal_pressure_boundary.resolve(1e-6f); + require(std::abs(equal_pressure_pipes.front().right_ghost.rho - equal_pressure_port.density) < 1e-6f, + "equal-pressure orifices must retain the volume-side density in the ghost state"); + + Volume0D configurable_volume{1e-3f, 150000.0f, 350.0f, 300.0f, 1.3f}; + auto& configurable_port = configurable_volume.create_port(); + require(std::abs(configurable_volume.temperature() - 350.0f) < 0.01f && + std::abs(configurable_port.enthalpy - 1.3f / 0.3f * 300.0f * 350.0f) < 1.0f, + "0-D volumes must honor configurable gamma and gas constant"); + configurable_volume.set_pressure(200000.0f, 400.0f); + require(std::abs(configurable_volume.pressure() - 200000.0f) < 0.1f && + std::abs(configurable_volume.temperature() - 400.0f) < 1e-3f, + "setting volume pressure and temperature must update conserved mass and energy"); + configurable_volume.energy_relaxation_rate = 0.0f; + configurable_volume.volume_change_rate = 0.01f; + const float before_volume_work = configurable_volume.pressure(); + configurable_volume.update_state(1e-3f); + require(configurable_volume.pressure() < before_volume_work && + std::abs(configurable_volume.pressure() - before_volume_work * 0.997f) < 1.0f, + "moving-volume work must remove p*dV energy using the configured volume rate"); + + Pipe1D profiled_pipe{{0.001f, 0.002f, 0.003f}, {0.1f, 0.2f, 0.3f}}; + require(std::abs(profiled_pipe.cell_mass(1) - 1.225f * 0.002f * 0.2f) < 1e-8f && + std::abs(profiled_pipe.minimum_cell_length() - 0.1f) < 1e-8f, + "profiled pipes must retain each cell's area and length for mass/CFL calculations"); + + Vehicle bike; + bike.set_first_gear(); + const float rpm = 10000.0f; + const float first_gear_speed = (rpm * 2.0f * pi / 60.0f) / bike.total_gear_ratio() * bike.wheel_radius * 3.6f; + require(first_gear_speed > 55.0f && first_gear_speed < 65.0f, + "first gear should be geared for about 60 km/h at the rev limit"); + Vehicle neutral; + neutral.update(500.0f, 0.04f, 1e-3f); + require(neutral.speed == 0.0f, "neutral must isolate the driveline"); + bike.clutch_input = 0.0f; + for (int i = 0; i < 2000; ++i) bike.update(500.0f, 0.04f, 1e-3f); + require(bike.speed > 0.0f && bike.kinetic_energy() > 0.0f, + "clutch impulse must transfer engine angular momentum into the vehicle"); + + FourStrokeTwinScenario scenario; + scenario.throttle = 0.35f; + for (int i = 0; i < 10000; ++i) scenario.process(); + require(std::isfinite(scenario.crankshaft.omega) && scenario.crankshaft.omega > 0.0f, + "twin scenario must retain a finite running crank speed"); + require(std::isfinite(scenario.intake_plenum.pressure()) && std::isfinite(scenario.exhaust_collector.pressure()), + "shared intake and exhaust volumes must remain finite"); + require(scenario.cylinder1.pressure > 0.0f && scenario.cylinder2.pressure > 0.0f, + "both cylinders must retain positive pressure"); + + scenario.toggle_mechanical_sound(); + for (int i = 0; i < 1000; ++i) { + const AudioFrame frame = scenario.process(); + require(std::isfinite(frame.left) && std::isfinite(frame.right), + "the integrated mechanical/reverb audio path must remain finite"); + } + require(scenario.snapshot().mechanical_sound_enabled, + "mechanical-noise toggle must be reflected in the scenario snapshot"); + + std::cout << "FluidSim C++ physics tests passed\n"; + return 0; +} diff --git a/tools/parity_trace.cpp b/tools/parity_trace.cpp new file mode 100644 index 0000000..06abff9 --- /dev/null +++ b/tools/parity_trace.cpp @@ -0,0 +1,83 @@ +#include "FluidSim/Scenario.hpp" + +#include +#include +#include + +int main() { + using namespace fluidsim::port; + constexpr int sample_rate = 44100; + constexpr int sample_count = sample_rate; + std::vector checkpoints; + checkpoints.reserve(106); + for (int step = 1; step <= 100; ++step) checkpoints.push_back(step); + checkpoints.insert(checkpoints.end(), {300, 500, 1000, 4410, 11025, 22050, 33075, 44100}); + + FourStrokeTwinScenario scenario(sample_rate); + scenario.set_controls(0.35f, 0.0f, 0.0f); + scenario.cylinder1.energy_variation = 0.0f; + scenario.cylinder2.energy_variation = 0.0f; + + std::cout << "step,rpm,angle,cylinder1_pa,cylinder2_pa,intake_pa,exhaust_pa,speed_kmh," + "valve1_in,valve2_in,intake_flow,exhaust_flow,pipe0_left_pa,pipe0_right_pa," + "pipe1_left_pa,pipe1_right_pa,pipe2_left_pa,pipe2_right_pa," + "pipe3_left_pa,pipe3_right_pa,pipe4_left_pa,pipe4_right_pa," + "throttle_port_flow,runner1_port_flow,runner2_port_flow," + "cyl1_intake_flow,cyl2_intake_flow,cyl1_exhaust_flow,cyl2_exhaust_flow," + "pipe0_right_rho,pipe0_right_u,pipe1_left_rho,pipe1_left_u,pipe1_right_rho,pipe1_right_u," + "pipe2_left_rho,pipe2_left_u,pipe2_right_rho,pipe2_right_u," + "intake_rho,intake_temperature,pipe0_left_rho,pipe0_left_u\n"; + std::cout << std::setprecision(10); + std::size_t checkpoint = 0; + for (int step = 1; step <= sample_count; ++step) { + scenario.process(); + if (checkpoint < checkpoints.size() && step == checkpoints[checkpoint]) { + const auto state = scenario.snapshot(); + std::cout << step << ',' + << scenario.crankshaft.omega * 60.0f / (2.0f * pi) << ',' + << scenario.crankshaft.angle << ',' + << scenario.cylinder1.pressure << ',' + << scenario.cylinder2.pressure << ',' + << scenario.intake_plenum.pressure() << ',' + << scenario.exhaust_collector.pressure() << ',' + << scenario.vehicle.speed * 3.6f << ',' + << state.valve1_intake_area << ',' + << state.valve2_intake_area << ',' + << state.intake_flow << ',' + << state.exhaust_flow << ',' + << state.pipe_states[0].front().pressure << ',' + << state.pipe_states[0].back().pressure << ',' + << state.pipe_states[1].front().pressure << ',' + << state.pipe_states[1].back().pressure << ',' + << state.pipe_states[2].front().pressure << ',' + << state.pipe_states[2].back().pressure << ',' + << state.pipe_states[3].front().pressure << ',' + << state.pipe_states[3].back().pressure << ',' + << state.pipe_states[4].front().pressure << ',' + << state.pipe_states[4].back().pressure << ',' + << scenario.intake_plenum.ports()[0]->mass_flow << ',' + << scenario.intake_plenum.ports()[1]->mass_flow << ',' + << scenario.intake_plenum.ports()[2]->mass_flow << ',' + << scenario.cylinder1.intake.mass_flow << ',' + << scenario.cylinder2.intake.mass_flow << ',' + << scenario.cylinder1.exhaust.mass_flow << ',' + << scenario.cylinder2.exhaust.mass_flow << ',' + << state.pipe_states[0].back().rho << ',' + << state.pipe_states[0].back().velocity << ',' + << state.pipe_states[1].front().rho << ',' + << state.pipe_states[1].front().velocity << ',' + << state.pipe_states[1].back().rho << ',' + << state.pipe_states[1].back().velocity << ',' + << state.pipe_states[2].front().rho << ',' + << state.pipe_states[2].front().velocity << ',' + << state.pipe_states[2].back().rho << ',' + << state.pipe_states[2].back().velocity << ',' + << scenario.intake_plenum.density() << ',' + << scenario.intake_plenum.temperature() << ',' + << state.pipe_states[0].front().rho << ',' + << state.pipe_states[0].front().velocity << '\n'; + ++checkpoint; + } + } + return 0; +} diff --git a/verification/CSharpParityTrace.csproj b/verification/CSharpParityTrace.csproj new file mode 100644 index 0000000..15e98c4 --- /dev/null +++ b/verification/CSharpParityTrace.csproj @@ -0,0 +1,11 @@ + + + Exe + net10.0 + enable + enable + + + + + diff --git a/verification/Compare-ScenarioTrace.ps1 b/verification/Compare-ScenarioTrace.ps1 new file mode 100644 index 0000000..f827616 --- /dev/null +++ b/verification/Compare-ScenarioTrace.ps1 @@ -0,0 +1,72 @@ +[CmdletBinding()] +param( + [string]$BuildDirectory = (Join-Path $PSScriptRoot '..\build-msvc'), + [double]$Tolerance = 0.01, + [switch]$ShowAllDifferences +) + +$cppRunner = Join-Path $BuildDirectory 'fluidsim_parity_trace.exe' +$csProject = Join-Path $PSScriptRoot 'CSharpParityTrace.csproj' +if (-not (Test-Path -LiteralPath $cppRunner)) { throw "C++ trace runner not found: $cppRunner" } + +$cppText = & $cppRunner +if ($LASTEXITCODE -ne 0) { throw "C++ trace runner failed with exit code $LASTEXITCODE." } +$csText = & dotnet run --project $csProject --configuration Release --no-launch-profile --verbosity quiet +if ($LASTEXITCODE -ne 0) { throw "C# reference trace failed with exit code $LASTEXITCODE." } + +$cppRows = @($cppText | Where-Object { $_ -match '^(step,|\d+,)' } | ConvertFrom-Csv) +$csRows = @($csText | Where-Object { $_ -match '^(step,|\d+,)' } | ConvertFrom-Csv) +if ($cppRows.Count -ne $csRows.Count) { throw "Checkpoint count differs: C++=$($cppRows.Count), C#=$($csRows.Count)." } + +$columns = @('rpm', 'angle', 'cylinder1_pa', 'cylinder2_pa', 'intake_pa', 'exhaust_pa', 'speed_kmh', + 'valve1_in', 'valve2_in', 'intake_flow', 'exhaust_flow', 'pipe0_left_pa', 'pipe0_right_pa', + 'pipe1_left_pa', 'pipe1_right_pa', 'pipe2_left_pa', 'pipe2_right_pa', + 'pipe3_left_pa', 'pipe3_right_pa', 'pipe4_left_pa', 'pipe4_right_pa', + 'throttle_port_flow', 'runner1_port_flow', 'runner2_port_flow', + 'cyl1_intake_flow', 'cyl2_intake_flow', 'cyl1_exhaust_flow', 'cyl2_exhaust_flow', + 'pipe0_right_rho', 'pipe0_right_u', 'pipe1_left_rho', 'pipe1_left_u', 'pipe1_right_rho', 'pipe1_right_u', + 'pipe2_left_rho', 'pipe2_left_u', 'pipe2_right_rho', 'pipe2_right_u', + 'intake_rho', 'intake_temperature', 'pipe0_left_rho', 'pipe0_left_u') +$largestDifference = 0.0 +$largestField = '' +$largestStep = 0 +$differences = [System.Collections.Generic.List[object]]::new() +for ($rowIndex = 0; $rowIndex -lt $cppRows.Count; ++$rowIndex) { + if ($cppRows[$rowIndex].step -ne $csRows[$rowIndex].step) { throw 'Checkpoint steps differ.' } + foreach ($column in $columns) { + $cppValue = [double]::Parse($cppRows[$rowIndex].$column, [Globalization.CultureInfo]::InvariantCulture) + $csValue = [double]::Parse($csRows[$rowIndex].$column, [Globalization.CultureInfo]::InvariantCulture) + $difference = [Math]::Abs($cppValue - $csValue) + $scale = [Math]::Max(1.0, [Math]::Abs($csValue)) + $relativeDifference = $difference / $scale + if ($relativeDifference -gt $largestDifference) { + $largestDifference = $relativeDifference + $largestField = $column + $largestStep = [int]$cppRows[$rowIndex].step + } + if ($relativeDifference -gt 0.0) { + $differences.Add([pscustomobject]@{ + Step = [int]$cppRows[$rowIndex].step + Field = $column + CSharp = $csValue + Cpp = $cppValue + RelativeDifference = $relativeDifference + }) + } + if ($ShowAllDifferences) { + "step={0} {1}: C#={2:G9}, C++={3:G9}, rel={4:G6}" -f ` + $cppRows[$rowIndex].step, $column, $csValue, $cppValue, $relativeDifference + } + } +} + +"Compared {0} checkpoints across {1} fields; tolerance {2:G6}." -f $cppRows.Count, $columns.Count, $Tolerance +"Largest relative difference: {0:G6} in {1} at step {2}." -f $largestDifference, $largestField, $largestStep +if ($largestDifference -gt $Tolerance) { + 'Largest remaining differences:' + $differences | Sort-Object RelativeDifference -Descending | Select-Object -First 10 | ForEach-Object { + " step={0} {1}: C#={2:G9}, C++={3:G9}, rel={4:G6}" -f ` + $_.Step, $_.Field, $_.CSharp, $_.Cpp, $_.RelativeDifference + } +} +if ($largestDifference -gt $Tolerance) { exit 1 } diff --git a/verification/Program.cs b/verification/Program.cs new file mode 100644 index 0000000..4650a4f --- /dev/null +++ b/verification/Program.cs @@ -0,0 +1,94 @@ +using System.Globalization; +using System.Linq; +using System.Reflection; +using FluidSim.Core; +using FluidSim.Components; +using FluidSim.Tests; + +static T PrivateField(object target, string name) +{ + var field = target.GetType().GetField(name, BindingFlags.Instance | BindingFlags.NonPublic) + ?? throw new InvalidOperationException($"Could not find scenario field '{name}'."); + return (T)(field.GetValue(target) ?? throw new InvalidOperationException($"Field '{name}' is null.")); +} + +const int sampleRate = 44100; +const int sampleCount = sampleRate; +int[] checkpoints = [.. Enumerable.Range(1, 100), 300, 500, 1000, 4410, 11025, 22050, 33075, 44100]; +var scenario = new SingleCylScenario(); +var originalOutput = Console.Out; +Console.SetOut(TextWriter.Null); +scenario.Initialize(sampleRate); +scenario.Throttle = 0.35f; +scenario.Load = 0.0f; +scenario.Clutch = 0.0f; + +var crankshaft = PrivateField(scenario, "crankshaft"); +var cylinder1 = PrivateField(scenario, "cylinder"); +var cylinder2 = PrivateField(scenario, "cylinder2"); +var intakePlenum = PrivateField(scenario, "intakePlenum"); +var exhaustCollector = PrivateField(scenario, "exhaustCollector"); +var vehicle = PrivateField(scenario, "vehicle"); +var boundaries = PrivateField(scenario, "boundaries"); +var pipeSystem = PrivateField(scenario, "pipeSystem"); +cylinder1.EnergyVariationFraction = 0.0f; +cylinder2.EnergyVariationFraction = 0.0f; + +Console.SetOut(originalOutput); +Console.WriteLine("step,rpm,angle,cylinder1_pa,cylinder2_pa,intake_pa,exhaust_pa,speed_kmh,valve1_in,valve2_in,intake_flow,exhaust_flow,pipe0_left_pa,pipe0_right_pa,pipe1_left_pa,pipe1_right_pa,pipe2_left_pa,pipe2_right_pa,pipe3_left_pa,pipe3_right_pa,pipe4_left_pa,pipe4_right_pa,throttle_port_flow,runner1_port_flow,runner2_port_flow,cyl1_intake_flow,cyl2_intake_flow,cyl1_exhaust_flow,cyl2_exhaust_flow,pipe0_right_rho,pipe0_right_u,pipe1_left_rho,pipe1_left_u,pipe1_right_rho,pipe1_right_u,pipe2_left_rho,pipe2_left_u,pipe2_right_rho,pipe2_right_u,intake_rho,intake_temperature,pipe0_left_rho,pipe0_left_u"); +var checkpointIndex = 0; +Console.SetOut(TextWriter.Null); +for (var step = 1; step <= sampleCount; ++step) +{ + scenario.Process(); + if (checkpointIndex < checkpoints.Length && step == checkpoints[checkpointIndex]) + { + Console.SetOut(originalOutput); + Console.WriteLine(string.Join(',', + step.ToString(CultureInfo.InvariantCulture), + (crankshaft.AngularVelocity * 60f / (2f * MathF.PI)).ToString("R", CultureInfo.InvariantCulture), + crankshaft.CrankAngle.ToString("R", CultureInfo.InvariantCulture), + cylinder1.Pressure.ToString("R", CultureInfo.InvariantCulture), + cylinder2.Pressure.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Pressure.ToString("R", CultureInfo.InvariantCulture), + exhaustCollector.Pressure.ToString("R", CultureInfo.InvariantCulture), + vehicle.SpeedKmh.ToString("R", CultureInfo.InvariantCulture), + cylinder1.IntakeValveArea.ToString("R", CultureInfo.InvariantCulture), + cylinder2.IntakeValveArea.ToString("R", CultureInfo.InvariantCulture), + boundaries.GetOpenEndMassFlow(0).ToString("R", CultureInfo.InvariantCulture), + boundaries.GetOpenEndMassFlow(1).ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(0).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(0).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(1).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(1).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(2).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(2).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(3).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(3).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(4).p.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(4).p.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Ports[0].MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Ports[2].MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Ports[3].MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + cylinder1.IntakePort.MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + cylinder2.IntakePort.MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + cylinder1.ExhaustPort.MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + cylinder2.ExhaustPort.MassFlowRate.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(0).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(0).u.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(1).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(1).u.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(1).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(1).u.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(2).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(2).u.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(2).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateRight(2).u.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Density.ToString("R", CultureInfo.InvariantCulture), + intakePlenum.Temperature.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(0).rho.ToString("R", CultureInfo.InvariantCulture), + pipeSystem.GetInteriorStateLeft(0).u.ToString("R", CultureInfo.InvariantCulture))); + Console.SetOut(TextWriter.Null); + ++checkpointIndex; + } +}