1175 lines
64 KiB
C++
1175 lines
64 KiB
C++
#pragma once
|
|
|
|
#include "FluidSim/Engine.hpp"
|
|
#include "FluidSim/OpenEndRadiation.hpp"
|
|
#include "FluidSim/PipeSolver.hpp"
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <span>
|
|
#include <cmath>
|
|
#include <deque>
|
|
#include <limits>
|
|
#include <stdexcept>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
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)) {
|
|
if (!std::isfinite(volume_m3) || volume_m3 <= 0 || !std::isfinite(pressure) || pressure <= 0 ||
|
|
!std::isfinite(temperature) || temperature <= 0 || !std::isfinite(gas_constant_value) ||
|
|
gas_constant_value <= 0 || !std::isfinite(gamma_value) || gamma_value <= 1)
|
|
throw std::invalid_argument("volume thermodynamic state must be finite and positive, gamma > 1");
|
|
}
|
|
|
|
float gamma{1.4f};
|
|
float gas_constant{287.0f};
|
|
float ambient_pressure{101325.0f};
|
|
float ambient_temperature{300.0f};
|
|
float energy_relaxation_rate{0.0f};
|
|
float volume_change_rate{};
|
|
|
|
Port &create_port() {
|
|
ports_.emplace_back();
|
|
auto &p = ports_.back();
|
|
p.owner_energy = &energy_;
|
|
p.owner_air_mass = &air_mass_;
|
|
p.owner_exhaust_mass = &exhaust_mass_;
|
|
p.owner_fuel_mass = &fuel_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_ + fuel_mass_; }
|
|
float internal_energy() const { return energy_; }
|
|
float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); }
|
|
float fuel_fraction() const { return fuel_mass_ / std::max(mass(), 1e-12f); }
|
|
float fuel_mass() const { return fuel_mass_; }
|
|
float fuel_equivalent_fraction()const{return fuel_equivalent_mass_/std::max(mass(),1e-12f);}
|
|
void add_fuel_mass(float value, float temperature = 300.0f) {
|
|
if (!std::isfinite(value) || value < 0.0f || !std::isfinite(temperature) || temperature <= 0.0f)
|
|
throw std::invalid_argument("invalid fuel addition");
|
|
fuel_mass_ += value;
|
|
fuel_equivalent_mass_ += value;
|
|
energy_ += value * gas_constant / (gamma - 1.0f) * temperature;
|
|
sync_ports();
|
|
}
|
|
void set_pressure(float value, float temp = -1.0f) {
|
|
if (!std::isfinite(value) || value <= 0 || !std::isfinite(temp))
|
|
throw std::invalid_argument("invalid volume pressure/temperature");
|
|
const float target_t = temp > 0.0f ? temp : temperature();
|
|
const float total = value * volume_ / (gas_constant * target_t);
|
|
const float fraction = air_fraction();
|
|
const float fuel = fuel_fraction();
|
|
const float equivalent=fuel_equivalent_fraction();
|
|
air_mass_ = total * fraction;
|
|
fuel_mass_ = total * fuel;
|
|
exhaust_mass_ = total * std::max(1.0f - fraction - fuel, 0.0f);
|
|
fuel_equivalent_mass_=total*equivalent;equivalent_roundoff_=0;
|
|
energy_ = value * volume_ / (gamma - 1.0f);
|
|
air_roundoff_ = exhaust_roundoff_ = fuel_roundoff_ = energy_roundoff_ = 0;
|
|
sync_ports();
|
|
}
|
|
void update(float dt) {
|
|
if (!std::isfinite(dt) || dt < 0)
|
|
throw std::invalid_argument("invalid volume timestep");
|
|
if (dt == 0)
|
|
return;
|
|
double d_air = 0, d_exhaust = 0, d_fuel = 0, d_energy = 0, d_equivalent=0;
|
|
const float fraction = air_fraction();
|
|
const float fuel = fuel_fraction();
|
|
for (const Port &p : ports_) {
|
|
const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction;
|
|
const float flow_fuel = p.mass_flow >= 0.0f ? p.fuel_fraction : fuel;
|
|
d_air += double(p.mass_flow) * flow_fraction;
|
|
d_fuel += double(p.mass_flow) * flow_fuel;
|
|
d_exhaust += double(p.mass_flow) * std::max(1.0f - flow_fraction - flow_fuel, 0.0f);
|
|
d_energy += double(p.mass_flow) * p.enthalpy;
|
|
d_equivalent+=double(p.mass_flow)*(p.mass_flow>=0?std::max(p.fuel_equivalent_fraction,p.fuel_fraction):fuel_equivalent_fraction());
|
|
}
|
|
accumulate(air_mass_, air_roundoff_, d_air * dt);
|
|
accumulate(exhaust_mass_, exhaust_roundoff_, d_exhaust * dt);
|
|
accumulate(fuel_mass_, fuel_roundoff_, d_fuel * dt);
|
|
accumulate(energy_, energy_roundoff_, d_energy * dt);
|
|
accumulate(fuel_equivalent_mass_,equivalent_roundoff_,d_equivalent*dt);
|
|
const float new_volume = volume_ + volume_change_rate * dt;
|
|
if (!(new_volume > 0) || !std::isfinite(new_volume))
|
|
throw std::runtime_error("invalid moving volume");
|
|
if (new_volume != volume_) {
|
|
const double scale = std::pow(double(volume_) / new_volume, gamma - 1.0f);
|
|
const double precise = (double(energy_)+energy_roundoff_)*scale;
|
|
energy_ = float(precise);
|
|
energy_roundoff_ = precise-energy_;
|
|
}
|
|
volume_ = new_volume;
|
|
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);
|
|
energy_roundoff_ = 0;
|
|
}
|
|
if (!(mass() > 0) || !(energy_ > 0) || !std::isfinite(mass()) || !std::isfinite(energy_))
|
|
throw std::runtime_error(
|
|
"nonphysical volume state; reduce timestep or check external port fluxes");
|
|
if (air_mass_ < -1e-12f || exhaust_mass_ < -1e-12f || fuel_mass_ < -1e-12f)
|
|
throw std::runtime_error("negative volume species mass");
|
|
air_mass_ = std::max(air_mass_, 0.0f);
|
|
exhaust_mass_ = std::max(exhaust_mass_, 0.0f);
|
|
fuel_mass_ = std::max(fuel_mass_, 0.0f);
|
|
fuel_equivalent_mass_=std::clamp(fuel_equivalent_mass_,fuel_mass_,mass());
|
|
if(fuel_equivalent_mass_<std::numeric_limits<float>::min()){fuel_equivalent_mass_=0;equivalent_roundoff_=0;}
|
|
sync_ports();
|
|
}
|
|
void update_state(float dt) override { update(dt); }
|
|
std::vector<Port *> ports() override {
|
|
std::vector<Port *> result;
|
|
result.reserve(ports_.size());
|
|
for (Port &port : ports_)
|
|
result.push_back(&port);
|
|
return result;
|
|
}
|
|
const std::deque<Port> &ports() const { return ports_; }
|
|
|
|
private:
|
|
float volume_, air_mass_{}, exhaust_mass_{}, fuel_mass_{}, energy_;
|
|
float fuel_equivalent_mass_{};double equivalent_roundoff_{};
|
|
// Preserve sub-ulp transfers without changing the float port/budget ABI.
|
|
double air_roundoff_{}, exhaust_roundoff_{}, fuel_roundoff_{}, energy_roundoff_{};
|
|
static void accumulate(float& value, double& roundoff, double increment) {
|
|
const double precise = double(value)+roundoff+increment;
|
|
value = float(precise);
|
|
roundoff = precise-value;
|
|
}
|
|
std::deque<Port> ports_;
|
|
void sync_port(Port &p) const {
|
|
p.gamma_value = gamma;
|
|
p.gas_constant_value = gas_constant;
|
|
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();
|
|
p.fuel_fraction = fuel_fraction();
|
|
p.fuel_equivalent_fraction=fuel_equivalent_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};
|
|
float fuel{};
|
|
float fuel_equivalent{};
|
|
};
|
|
struct GasState {
|
|
float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f};
|
|
float fuel{};
|
|
float fuel_equivalent{};
|
|
};
|
|
|
|
class Pipe1D {
|
|
public:
|
|
std::vector<float> linear_resistance, quadratic_resistance;
|
|
// Integrated face fluxes: kg/s, N, W, kg/s of fresh air and fuel, positive along +x.
|
|
struct Flux {
|
|
float m{}, p{}, e{}, y{}, fuel{}, fuel_equivalent{};
|
|
};
|
|
struct RoeState { double root_density{}, enthalpy{}; };
|
|
PipeSolver solver{PipeSolver::Rusanov};
|
|
std::uint64_t hllc_fallbacks{};
|
|
float area{}, dx{};
|
|
float damping_multiplier{0.0f}, energy_relaxation_rate{0.0f}, ambient_pressure{101325.0f};
|
|
std::vector<float> cell_areas, cell_lengths;
|
|
std::vector<GasCell> cells;
|
|
GasState left_ghost{}, right_ghost{};
|
|
std::uint64_t positivity_repairs{};
|
|
std::size_t scratch_bytes() const {
|
|
return next_cells_.capacity()*sizeof(GasCell)+states_.capacity()*sizeof(GasState)+
|
|
reconstructed_left_.capacity()*sizeof(GasState)+reconstructed_right_.capacity()*sizeof(GasState)+
|
|
waves_.capacity()*sizeof(float)+faces_.capacity()*sizeof(Flux)+roe_states_.capacity()*sizeof(RoeState);
|
|
}
|
|
Pipe1D(float a, float length, int count)
|
|
: Pipe1D(uniform(a, count), uniform(length / std::max(count, 1), count)) {}
|
|
Pipe1D(std::vector<float> areas, std::vector<float> lengths)
|
|
: area(areas.empty() ? 0 : areas.front()), dx(lengths.empty() ? 0 : lengths.front()),
|
|
cell_areas(std::move(areas)), cell_lengths(std::move(lengths)), cells(cell_areas.size()),
|
|
next_cells_(cells.size()), states_(cells.size()), waves_(cells.size()), faces_(cells.size() + 1) {
|
|
if (cells.empty() || cell_lengths.size() != cells.size())
|
|
throw std::invalid_argument("pipe profiles must be nonempty and have equal length");
|
|
for (std::size_t i = 0; i < cells.size(); ++i)
|
|
if (!std::isfinite(cell_areas[i]) || !std::isfinite(cell_lengths[i]) || cell_areas[i] <= 0 ||
|
|
cell_lengths[i] <= 0)
|
|
throw std::invalid_argument("pipe geometry must be finite and positive");
|
|
if(std::any_of(cell_areas.begin(),cell_areas.end(),[&](float a){return a!=cell_areas.front();})){
|
|
geometry_cfl_lengths_.resize(cells.size());
|
|
for(std::size_t i=0;i<cells.size();++i){const float a=cell_areas[i],left=i?.5f*(cell_areas[i-1]+a):a,right=i+1<cells.size()?.5f*(a+cell_areas[i+1]):a;
|
|
geometry_cfl_lengths_[i]=cell_lengths[i]*(a/std::max({a,left,right}));}
|
|
}
|
|
}
|
|
GasState state(int i) const { return primitive(cells.at(static_cast<std::size_t>(i))); }
|
|
float cell_mass(int i) const { return cells.at(i).rho * cell_areas.at(i) * cell_lengths.at(i); }
|
|
float cell_area(int i) const { return cell_areas.at(i); }
|
|
float cell_length(int i) const { return cell_lengths.at(i); }
|
|
float minimum_cell_length() const { return *std::min_element(cell_lengths.begin(), cell_lengths.end()); }
|
|
float stable_timestep(float cfl = 0.8f) const {
|
|
float dt = std::numeric_limits<float>::max();
|
|
for (std::size_t i = 0; i < cells.size(); ++i) {
|
|
auto s = primitive(cells[i]);
|
|
if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
|
|
throw std::runtime_error("nonphysical pipe state in CFL calculation");
|
|
float speed=wave(s),length=cell_lengths[i];
|
|
if(!geometry_cfl_lengths_.empty()){length=geometry_cfl_lengths_[i];if(i)speed=std::max(speed,wave(primitive(cells[i-1])));if(i+1<cells.size())speed=std::max(speed,wave(primitive(cells[i+1])));}
|
|
dt = std::min(dt, cfl * length / speed);
|
|
}
|
|
return dt;
|
|
}
|
|
// Solver preparation is valid until advance; boundary resolution changes
|
|
// fluxes/ghosts, not cells. Standalone advance still prepares its own cache.
|
|
float prepare_step(float cfl = 0.8f) {
|
|
if(!scratch_solver_initialized_ || scratch_solver_!=solver) {
|
|
if(solver==PipeSolver::Muscl) {
|
|
if(reconstructed_left_.size()!=cells.size()) reconstructed_left_.resize(cells.size());
|
|
if(reconstructed_right_.size()!=cells.size()) reconstructed_right_.resize(cells.size());
|
|
} else if(solver==PipeSolver::Hlle || solver==PipeSolver::Hllc) {
|
|
if(roe_states_.size()!=cells.size()) roe_states_.resize(cells.size());
|
|
}
|
|
scratch_solver_=solver;
|
|
scratch_solver_initialized_=true;
|
|
}
|
|
float dt = std::numeric_limits<float>::max();
|
|
for (std::size_t i = 0; i < cells.size(); ++i) {
|
|
states_[i] = primitive(cells[i]);
|
|
waves_[i] = wave(states_[i]);
|
|
}
|
|
// Keep the arithmetic loop free of exceptional control flow so SIMD
|
|
// can evaluate independent cells, then validate before any flux update.
|
|
for (std::size_t i = 0; i < cells.size(); ++i) {
|
|
const auto &s = states_[i];
|
|
if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
|
|
throw std::runtime_error("nonphysical pipe state in preparation");
|
|
dt = std::min(dt, cfl * cell_lengths[i] / waves_[i]);
|
|
}
|
|
// A narrow cell exchanges flux through faces larger than its own area.
|
|
// Limit that face/volume ratio and use each face's fastest adjacent wave.
|
|
// Uniform pipes retain the original fast preparation path.
|
|
if(!geometry_cfl_lengths_.empty())for(std::size_t i=0;i<cells.size();++i){
|
|
float speed=waves_[i];if(i)speed=std::max(speed,waves_[i-1]);if(i+1<cells.size())speed=std::max(speed,waves_[i+1]);
|
|
dt=std::min(dt,cfl*geometry_cfl_lengths_[i]/speed);
|
|
}
|
|
if (solver == PipeSolver::Hlle || solver == PipeSolver::Hllc)
|
|
for (std::size_t i=0; i<cells.size(); ++i) {
|
|
const auto& s=states_[i];
|
|
const double energy=s.pressure*pipe_gm1_inverse+.5*double(s.rho)*s.velocity*s.velocity;
|
|
roe_states_[i]={std::sqrt(s.rho),(energy+s.pressure)/s.rho};
|
|
}
|
|
prepared_ = true;
|
|
return dt;
|
|
}
|
|
int required_substeps(float dt, float cfl = 0.8f) const {
|
|
if (!std::isfinite(dt) || dt < 0 || !std::isfinite(cfl) || cfl <= 0 || cfl > 1)
|
|
throw std::invalid_argument("invalid timestep/CFL");
|
|
const double n = std::ceil(double(dt) / stable_timestep(cfl));
|
|
if (n > 1000000)
|
|
throw std::runtime_error("excessive pipe substeps");
|
|
return std::max(1, static_cast<int>(n));
|
|
}
|
|
void set_boundary_flux(bool left, Flux f) {
|
|
if (left) {
|
|
left_flux_ = f;
|
|
left_override_ = true;
|
|
} else {
|
|
right_flux_ = f;
|
|
right_override_ = true;
|
|
}
|
|
}
|
|
static GasState primitive(const GasCell &c) {
|
|
const float u = c.momentum / std::max(c.rho, 1e-12f);
|
|
const float air = std::clamp(c.air, 0.0f, 1.0f);
|
|
return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), air,
|
|
std::clamp(c.fuel, 0.0f, 1.0f - air),std::clamp(std::max(c.fuel,c.fuel_equivalent),0.f,1.f)};
|
|
}
|
|
static float wave(const GasState &s) {
|
|
return std::abs(s.velocity) + std::sqrt(gamma * s.pressure / s.rho);
|
|
}
|
|
static Flux riemann(const GasState &a, const GasState &b) {
|
|
return riemann(a, b, std::max(wave(a), wave(b)));
|
|
}
|
|
// HLLE/HLLC use bounds enclosing endpoint and Roe acoustic speeds.
|
|
template<PipeSolver method>
|
|
static Flux flux_kernel(const GasState& a, const GasState& b,
|
|
float ca, float cb, bool* fallback = nullptr,
|
|
const RoeState* left_roe = nullptr, const RoeState* right_roe = nullptr) {
|
|
if constexpr (method == PipeSolver::Rusanov)
|
|
return riemann(a,b,std::max(std::abs(a.velocity)+ca,std::abs(b.velocity)+cb));
|
|
const double ra=left_roe?left_roe->root_density:std::sqrt(a.rho);
|
|
const double rb=right_roe?right_roe->root_density:std::sqrt(b.rho), sum=ra+rb;
|
|
const double ea=a.pressure*pipe_gm1_inverse+.5*double(a.rho)*a.velocity*a.velocity;
|
|
const double eb=b.pressure*pipe_gm1_inverse+.5*double(b.rho)*b.velocity*b.velocity;
|
|
const double ua=a.velocity, ub=b.velocity;
|
|
const double uroe=(ra*ua+rb*ub)/sum;
|
|
const double ha=left_roe?left_roe->enthalpy:(ea+a.pressure)/a.rho;
|
|
const double hb=right_roe?right_roe->enthalpy:(eb+b.pressure)/b.rho;
|
|
const double hroe=(ra*ha+rb*hb)/sum;
|
|
const double croe=std::sqrt(std::max(0.,double(gamma-1)*(hroe-.5*uroe*uroe)));
|
|
const double sl=std::min({ua-ca,ub-cb,uroe-croe});
|
|
const double sr=std::max({ua+ca,ub+cb,uroe+croe});
|
|
const auto physical=[](const GasState& q,double energy) {
|
|
const float m=q.rho*q.velocity;
|
|
return Flux{m,m*q.velocity+q.pressure,float((energy+q.pressure)*q.velocity),m*q.air,m*q.fuel,m*q.fuel_equivalent};
|
|
};
|
|
const Flux fa=physical(a,ea), fb=physical(b,eb);
|
|
if(sl>=0) return fa;
|
|
if(sr<=0) return fb;
|
|
const auto hlle=[&] {
|
|
const double inv=1/(sr-sl), diff=sl*sr;
|
|
return Flux{float((sr*fa.m-sl*fb.m+diff*(b.rho-a.rho))*inv),
|
|
float((sr*fa.p-sl*fb.p+diff*(b.rho*ub-a.rho*ua))*inv),
|
|
float((sr*fa.e-sl*fb.e+diff*(eb-ea))*inv),
|
|
float((sr*fa.y-sl*fb.y+diff*(b.rho*b.air-a.rho*a.air))*inv),
|
|
float((sr*fa.fuel-sl*fb.fuel+diff*(b.rho*b.fuel-a.rho*a.fuel))*inv),float((sr*fa.fuel_equivalent-sl*fb.fuel_equivalent+diff*(b.rho*b.fuel_equivalent-a.rho*a.fuel_equivalent))*inv)};
|
|
};
|
|
if constexpr(method==PipeSolver::Hlle) return hlle();
|
|
const double dl=a.rho*(sl-ua), dr=b.rho*(sr-ub);
|
|
const double sm=(b.pressure-a.pressure+dl*ua-dr*ub)/(dl-dr);
|
|
const double ps=a.pressure+dl*(sm-ua);
|
|
// Reject nonphysical HLLC intermediate states without modifying conserved cells.
|
|
const bool left=sm>=0;
|
|
const auto& q=left?a:b;
|
|
const auto f=left?fa:fb;
|
|
const double speed=left?sl:sr;
|
|
const double inv=1/(speed-sm);
|
|
const double density=(left?dl:dr)*inv;
|
|
const double energy=((speed-q.velocity)*(left?ea:eb)-q.pressure*q.velocity+ps*sm)*inv;
|
|
if(!std::isfinite(sm) || !(sl<sm && sm<sr) || !(ps>0) || !(density>0) ||
|
|
!(energy>.5*density*sm*sm)) {
|
|
if(fallback) *fallback=true;
|
|
return hlle();
|
|
}
|
|
const double mass=f.m+speed*(density-q.rho);
|
|
return {float(mass),float(f.p+speed*(density*sm-q.rho*q.velocity)),
|
|
float(f.e+speed*(energy-(left?ea:eb))),float(mass*q.air),float(mass*q.fuel),float(mass*q.fuel_equivalent)};
|
|
}
|
|
static Flux numerical_flux(const GasState& a,const GasState& b,PipeSolver solver,
|
|
float ca,float cb,bool* fallback=nullptr) {
|
|
if(solver==PipeSolver::Rusanov) return flux_kernel<PipeSolver::Rusanov>(a,b,ca,cb,fallback);
|
|
if(solver==PipeSolver::Hlle) return flux_kernel<PipeSolver::Hlle>(a,b,ca,cb,fallback);
|
|
return flux_kernel<PipeSolver::Hllc>(a,b,ca,cb,fallback);
|
|
}
|
|
static Flux numerical_flux(const GasState& a,const GasState& b,PipeSolver solver) {
|
|
return numerical_flux(a,b,solver,std::sqrt(gamma*a.pressure/a.rho),std::sqrt(gamma*b.pressure/b.rho));
|
|
}
|
|
void set_acoustic_neighbor(bool left,GasState state) {
|
|
if(left) { left_acoustic_neighbor_=state; has_left_acoustic_neighbor_=true; }
|
|
else { right_acoustic_neighbor_=state; has_right_acoustic_neighbor_=true; }
|
|
}
|
|
void prepare_acoustic_reconstruction(float dt) {
|
|
if(solver!=PipeSolver::Muscl) return;
|
|
if(!prepared_) prepare_step();
|
|
reconstruct(dt); reconstructed_prepared_=true;
|
|
}
|
|
GasState acoustic_edge(bool left) const {
|
|
return left?reconstructed_left_.front():reconstructed_right_.back();
|
|
}
|
|
void advance(float dt) {
|
|
if (dt <= 0)
|
|
return;
|
|
if (!std::isfinite(dt))
|
|
throw std::invalid_argument("nonfinite pipe timestep");
|
|
const std::size_t n = cells.size();
|
|
if (!prepared_)
|
|
prepare_step();
|
|
faces_[0] = left_override_ ? left_flux_ : scaled(numerical_flux(left_ghost, states_[0], solver), cell_areas[0]);
|
|
faces_[n] =
|
|
right_override_ ? right_flux_ : scaled(numerical_flux(states_[n - 1], right_ghost, solver), cell_areas[n - 1]);
|
|
if(solver==PipeSolver::Muscl) {
|
|
if(!reconstructed_prepared_) reconstruct(dt);
|
|
for(std::size_t i=1;i<n;++i) {
|
|
bool fallback=false;
|
|
const auto& a=reconstructed_right_[i-1];
|
|
const auto& b=reconstructed_left_[i];
|
|
faces_[i]=scaled(numerical_flux(a,b,PipeSolver::Hllc,
|
|
std::sqrt(gamma*a.pressure/a.rho),std::sqrt(gamma*b.pressure/b.rho),&fallback),
|
|
.5f*(cell_areas[i-1]+cell_areas[i]));
|
|
hllc_fallbacks+=fallback;
|
|
}
|
|
} else if (solver == PipeSolver::Rusanov) {
|
|
for (std::size_t i = 1; i < n; ++i)
|
|
faces_[i] = scaled(riemann(states_[i-1], states_[i], std::max(waves_[i-1], waves_[i])),
|
|
.5f*(cell_areas[i-1]+cell_areas[i]));
|
|
} else {
|
|
const auto fill_faces=[&]<PipeSolver method>() {
|
|
for (std::size_t i = 1; i < n; ++i) {
|
|
bool fallback=false;
|
|
faces_[i] = scaled(flux_kernel<method>(states_[i-1], states_[i],
|
|
waves_[i-1]-std::abs(states_[i-1].velocity), waves_[i]-std::abs(states_[i].velocity),
|
|
&fallback, &roe_states_[i-1], &roe_states_[i]),
|
|
.5f*(cell_areas[i-1]+cell_areas[i]));
|
|
hllc_fallbacks += fallback;
|
|
}
|
|
};
|
|
if(solver==PipeSolver::Hlle) fill_faces.template operator()<PipeSolver::Hlle>();
|
|
else fill_faces.template operator()<PipeSolver::Hllc>();
|
|
}
|
|
const float relaxation = energy_relaxation_rate > 0 ? std::exp(-energy_relaxation_rate * dt) : 1.0f;
|
|
for (std::size_t i = 0; i < n; ++i) {
|
|
const auto fl = faces_[i], fr = faces_[i + 1];
|
|
const float factor = dt / (cell_areas[i] * cell_lengths[i]);
|
|
const float al = i ? .5f * (cell_areas[i - 1] + cell_areas[i]) : cell_areas[i];
|
|
const float ar = i + 1 < n ? .5f * (cell_areas[i] + cell_areas[i + 1]) : cell_areas[i];
|
|
float rho = cells[i].rho - factor * (fr.m - fl.m);
|
|
float mom = cells[i].momentum - factor * (fr.p - fl.p - states_[i].pressure * (ar - al));
|
|
float en = cells[i].energy - factor * (fr.e - fl.e);
|
|
const float air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y);
|
|
const float fuel = cells[i].rho * cells[i].fuel - factor * (fr.fuel - fl.fuel);
|
|
const float equivalent=cells[i].rho*std::max(cells[i].fuel_equivalent,cells[i].fuel)-factor*(fr.fuel_equivalent-fl.fuel_equivalent);
|
|
if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) ||
|
|
!std::isfinite(air) || !std::isfinite(fuel))
|
|
throw std::runtime_error("nonfinite finite-volume update");
|
|
if (rho < 1e-12f) {
|
|
rho = 1e-12f;
|
|
++positivity_repairs;
|
|
}
|
|
if (!linear_resistance.empty()) {
|
|
const float speed = std::abs(mom/rho);
|
|
mom /= 1 + dt*(linear_resistance[i] + .5f*quadratic_resistance[i]*speed);
|
|
// Retain total energy: dissipated kinetic energy becomes heat.
|
|
}
|
|
if (damping_multiplier > 0)
|
|
mom *= std::exp(-.1f * damping_multiplier * dt / cells[i].rho);
|
|
if (energy_relaxation_rate > 0)
|
|
en = ambient_pressure * pipe_gm1_inverse +
|
|
(en - ambient_pressure * pipe_gm1_inverse) * relaxation;
|
|
const float minimum = 1e-6f * pipe_gm1_inverse + .5f * mom * mom / rho;
|
|
if (en < minimum) {
|
|
en = minimum;
|
|
++positivity_repairs;
|
|
}
|
|
// Purged tracer tails below the smallest normal float carry less
|
|
// than 1.18e-38 of the gas mass. Avoid repeatedly transporting
|
|
// subnormal concentrations; total gas mass and energy are unchanged.
|
|
const float fraction = std::clamp(air / rho, 0.0f, 1.0f);
|
|
const float fuel_fraction = std::clamp(fuel / rho, 0.0f, 1.0f - fraction);
|
|
const float equivalent_fraction=std::clamp(equivalent/rho,fuel_fraction,1.f);
|
|
next_cells_[i] = {rho, mom, en,
|
|
fraction < std::numeric_limits<float>::min() ? 0.0f : fraction,
|
|
fuel_fraction < std::numeric_limits<float>::min() ? 0.0f : fuel_fraction,equivalent_fraction<std::numeric_limits<float>::min()?0.f:equivalent_fraction};
|
|
}
|
|
cells.swap(next_cells_);
|
|
prepared_ = false;
|
|
reconstructed_prepared_=has_left_acoustic_neighbor_=has_right_acoustic_neighbor_=false;
|
|
left_override_ = right_override_ = false;
|
|
}
|
|
static Flux scaled(Flux f, float a) {
|
|
return {f.m * a, f.p * a, f.e * a, f.y * a, f.fuel * a,f.fuel_equivalent*a};
|
|
}
|
|
|
|
private:
|
|
std::vector<GasCell> next_cells_;
|
|
std::vector<float> geometry_cfl_lengths_;
|
|
std::vector<GasState> states_;
|
|
std::vector<GasState> reconstructed_left_, reconstructed_right_;
|
|
std::vector<float> waves_;
|
|
std::vector<Flux> faces_;
|
|
std::vector<RoeState> roe_states_;
|
|
PipeSolver scratch_solver_{PipeSolver::Rusanov};
|
|
bool scratch_solver_initialized_{};
|
|
Flux left_flux_{}, right_flux_{};
|
|
bool left_override_{}, right_override_{}, prepared_{};
|
|
GasState left_acoustic_neighbor_{}, right_acoustic_neighbor_{};
|
|
bool has_left_acoustic_neighbor_{}, has_right_acoustic_neighbor_{}, reconstructed_prepared_{};
|
|
// Primitive minmod MUSCL-Hancock on uniform sections. Component interfaces
|
|
// and area/profile transitions retain the existing first-order coupling.
|
|
void reconstruct(float dt) {
|
|
const auto minmod=[](float a,float b) {
|
|
return a*b<=0 ? 0.f : std::copysign(std::min(std::abs(a),std::abs(b)),a);
|
|
};
|
|
for(std::size_t i=0;i<states_.size();++i) {
|
|
const auto q=states_[i];
|
|
reconstructed_left_[i]=reconstructed_right_[i]=q;
|
|
if((i==0 && !has_left_acoustic_neighbor_) ||
|
|
(i+1==states_.size() && !has_right_acoustic_neighbor_) ||
|
|
(i>0 && (cell_areas[i-1]!=cell_areas[i] || cell_lengths[i-1]!=cell_lengths[i])) ||
|
|
(i+1<states_.size() && (cell_areas[i+1]!=cell_areas[i] || cell_lengths[i+1]!=cell_lengths[i]))) continue;
|
|
const auto a=i==0?left_acoustic_neighbor_:states_[i-1];
|
|
const auto b=i+1==states_.size()?right_acoustic_neighbor_:states_[i+1];
|
|
const float dr=minmod(q.rho-a.rho,b.rho-q.rho);
|
|
const float du=minmod(q.velocity-a.velocity,b.velocity-q.velocity);
|
|
const float dp=minmod(q.pressure-a.pressure,b.pressure-q.pressure);
|
|
const float da=minmod(q.air-a.air,b.air-q.air);
|
|
const float df=minmod(q.fuel-a.fuel,b.fuel-q.fuel);
|
|
const float de=minmod(q.fuel_equivalent-a.fuel_equivalent,b.fuel_equivalent-q.fuel_equivalent);
|
|
const float h=.5f*dt/cell_lengths[i];
|
|
const GasState half{q.rho-h*(q.velocity*dr+q.rho*du),
|
|
q.velocity-h*(q.velocity*du+dp/q.rho),
|
|
q.pressure-h*(q.velocity*dp+gamma*q.pressure*du),
|
|
q.air-h*q.velocity*da,q.fuel-h*q.velocity*df,q.fuel_equivalent-h*q.velocity*de};
|
|
const auto edge=[&](float sign) { return GasState{half.rho+sign*dr,
|
|
half.velocity+sign*du,half.pressure+sign*dp,half.air+sign*da,half.fuel+sign*df,half.fuel_equivalent+sign*de}; };
|
|
const auto left=edge(-.5f), right=edge(.5f);
|
|
const auto physical=[](const GasState& v) {
|
|
return v.rho>1e-12f && v.pressure>1e-6f && std::isfinite(v.pressure) &&
|
|
std::isfinite(v.rho) && std::isfinite(v.velocity) && v.air>=0 && v.fuel>=0 && v.air+v.fuel<=1.000001f;
|
|
};
|
|
if(physical(left) && physical(right)) {
|
|
reconstructed_left_[i]=left; reconstructed_right_[i]=right;
|
|
}
|
|
}
|
|
}
|
|
static std::vector<float> uniform(float v, int n) {
|
|
if (n <= 0)
|
|
throw std::invalid_argument("pipe count must be positive");
|
|
return std::vector<float>(n, v);
|
|
}
|
|
static Flux riemann(const GasState &a, const GasState &b, float c) {
|
|
const float ma = a.rho * a.velocity, mb = b.rho * b.velocity;
|
|
const float ea = a.pressure * pipe_gm1_inverse + .5f * ma * a.velocity;
|
|
const float eb = b.pressure * pipe_gm1_inverse + .5f * mb * b.velocity;
|
|
return {.5f * (ma + mb - c * (b.rho - a.rho)),
|
|
.5f * (ma * a.velocity + a.pressure + mb * b.velocity + b.pressure - c * (mb - ma)),
|
|
.5f * ((ea + a.pressure) * a.velocity + (eb + b.pressure) * b.velocity - c * (eb - ea)),
|
|
.5f * (ma * a.air + mb * b.air - c * (b.rho * b.air - a.rho * a.air)),
|
|
.5f * (ma * a.fuel + mb * b.fuel - c * (b.rho * b.fuel - a.rho * a.fuel)),.5f*(ma*a.fuel_equivalent+mb*b.fuel_equivalent-c*(b.rho*b.fuel_equivalent-a.rho*a.fuel_equivalent))};
|
|
}
|
|
};
|
|
|
|
class BoundarySystem {
|
|
public:
|
|
explicit BoundarySystem(std::vector<Pipe1D> &pipes) : pipes_(pipes) {}
|
|
void add_orifice(Port &port, int pipe, bool left, float &area, float cd = 1.0f) {
|
|
if (std::abs(port.gamma_value - gamma) > 1e-6f ||
|
|
std::abs(port.gas_constant_value - gas_constant) > .001f)
|
|
throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
|
|
if (!std::isfinite(cd) || cd < 0 || cd > 1)
|
|
throw std::invalid_argument("orifice Cd must be in [0,1]");
|
|
for (auto &o : orifices_)
|
|
if (o.port == &port)
|
|
throw std::invalid_argument("each connection requires a distinct component port");
|
|
claim(pipe, left);
|
|
Orifice o;
|
|
o.port = &port;
|
|
o.pipe = pipe;
|
|
o.left = left;
|
|
o.area = &area;
|
|
o.cd = cd;
|
|
if (port.owner_air_mass) {
|
|
auto it = std::find_if(budgets_.begin(), budgets_.end(),
|
|
[&](const Budget &b) { return b.air == port.owner_air_mass; });
|
|
if (it == budgets_.end()) {
|
|
o.budget = budgets_.size();
|
|
budgets_.push_back({port.owner_air_mass, port.owner_exhaust_mass,
|
|
port.owner_fuel_mass, port.owner_energy});
|
|
} else
|
|
o.budget = static_cast<std::size_t>(it - budgets_.begin());
|
|
}
|
|
orifices_.push_back(o);
|
|
}
|
|
void add_orifice_with_inertance(Port &p, int pipe, bool left, float &a, float cd, float length,
|
|
float loss = 0) {
|
|
if (!std::isfinite(length) || length <= 0 || !std::isfinite(loss) || loss < 0)
|
|
throw std::invalid_argument("invalid inertance");
|
|
add_orifice(p, pipe, left, a, cd);
|
|
auto &o = orifices_.back();
|
|
o.length = length;
|
|
o.loss = loss;
|
|
}
|
|
void add_closed_end(int pipe, bool left) {
|
|
claim(pipe, left);
|
|
Orifice o;
|
|
o.pipe = pipe;
|
|
o.left = left;
|
|
orifices_.push_back(o);
|
|
}
|
|
void add_open_end(int pipe, bool left, float p, float a) {
|
|
if (!std::isfinite(p) || p <= 0 || !std::isfinite(a) || a <= 0)
|
|
throw std::invalid_argument("invalid open boundary");
|
|
claim(pipe, left);
|
|
opens_.push_back({pipe, left, p, a});
|
|
}
|
|
// Full circular mouth radiating into free space; no flange or mean-flow correction.
|
|
void add_radiating_open_end(int pipe, bool left, float pressure, float area, float temperature = 300) {
|
|
if (!std::isfinite(temperature) || temperature <= 0)
|
|
throw std::invalid_argument("invalid ambient temperature");
|
|
if (pipe < 0 || static_cast<std::size_t>(pipe) >= pipes_.size())
|
|
throw std::invalid_argument("invalid boundary pipe index");
|
|
const float face_area = pipes_[pipe].cell_area(end(pipes_[pipe], left));
|
|
if (std::abs(area-face_area) > 1e-6f*face_area)
|
|
throw std::invalid_argument("radiating mouth must match pipe face area");
|
|
add_open_end(pipe, left, pressure, area);
|
|
opens_.back().radiating = true;
|
|
opens_.back().temperature = temperature;
|
|
}
|
|
// An ideal, zero-volume junction between two duct faces. Orientation is explicit.
|
|
void add_connection(int a, bool left_a, int b, bool left_b) {
|
|
if (a < 0 || b < 0 || static_cast<std::size_t>(a) >= pipes_.size() ||
|
|
static_cast<std::size_t>(b) >= pipes_.size())
|
|
throw std::invalid_argument("invalid connection pipe index");
|
|
const float aa = pipes_[a].cell_area(end(pipes_[a], left_a));
|
|
const float ab = pipes_[b].cell_area(end(pipes_[b], left_b));
|
|
if (std::abs(aa - ab) > 1e-6f * std::max(aa, ab))
|
|
throw std::invalid_argument("unequal-area junction needs a profiled pipe transition");
|
|
if (a == b && left_a == left_b)
|
|
throw std::invalid_argument("cannot join a face to itself");
|
|
check_face(a, left_a);
|
|
check_face(b, left_b);
|
|
claim(a, left_a);
|
|
claim(b, left_b);
|
|
connections_.push_back({a, b, left_a, left_b});
|
|
}
|
|
struct JunctionFace {
|
|
int pipe; bool left;
|
|
// Unit vector pointing OUT of the node into this branch, in world axes.
|
|
std::array<float,3> direction{};
|
|
};
|
|
enum class JunctionMethod { Characteristic, LegacyRusanov };
|
|
void add_junction(std::initializer_list<JunctionFace> ends,
|
|
JunctionMethod method=JunctionMethod::Characteristic,float turning_loss=0) {
|
|
add_junction(std::span(ends.begin(),ends.size()),method,turning_loss);
|
|
}
|
|
void add_junction(std::span<const JunctionFace> ends,
|
|
JunctionMethod method=JunctionMethod::Characteristic,float turning_loss=0) {
|
|
if(ends.size()<2 || !std::isfinite(turning_loss) || turning_loss<0 || turning_loss>10)
|
|
throw std::invalid_argument("invalid junction size or turning loss");
|
|
JunctionNode node; node.method=method; node.turning_loss=turning_loss;
|
|
for(auto e:ends) {
|
|
check_face(e.pipe,e.left);
|
|
for(auto previous:node.ends)
|
|
if(previous.pipe==e.pipe && previous.left==e.left) throw std::invalid_argument("duplicate junction face");
|
|
double norm=0;
|
|
for(float v:e.direction) { if(!std::isfinite(v)) throw std::invalid_argument("invalid junction direction"); norm+=double(v)*v; }
|
|
if(norm>0) for(float& v:e.direction) v/=std::sqrt(norm);
|
|
else if(turning_loss>0) throw std::invalid_argument("direction required for turning loss");
|
|
node.ends.push_back(e);
|
|
}
|
|
for(auto e:ends) claim(e.pipe,e.left);
|
|
node.states.resize(ends.size()); node.faces.resize(ends.size());
|
|
nodes_.push_back(std::move(node));
|
|
}
|
|
float open_end_mass_flow(std::size_t i) const { return opens_.at(i).last_mass_flow; }
|
|
float open_end_volume_flow(std::size_t i) const { return opens_.at(i).last_volume_flow; }
|
|
float open_end_pressure(std::size_t i) const { return opens_.at(i).last_face_pressure; }
|
|
std::uint64_t limited_fluxes{}, nonlinear_iterations{}, junction_fallbacks{}, junction_iterations{}, junction_limited_fluxes{};
|
|
void validate_topology() const {
|
|
if (claims_.size() != 2 * pipes_.size())
|
|
throw std::runtime_error("every pipe endpoint needs an explicit boundary or connection");
|
|
}
|
|
void resolve(float dt) {
|
|
if (!std::isfinite(dt) || dt <= 0)
|
|
throw std::invalid_argument("boundary timestep must be positive");
|
|
for (auto &b : budgets_) {
|
|
b.mass = .25f * (*b.air + (b.exhaust ? *b.exhaust : 0) + (b.fuel ? *b.fuel : 0));
|
|
b.heat = b.energy ? .25f * (*b.energy) : std::numeric_limits<float>::max();
|
|
}
|
|
// Supply matching direct-pipe neighbors before reconstruction, preserving
|
|
// the same stencil across artificial pipe splits and all orientations.
|
|
for(const auto& o:connections_) {
|
|
auto& a=pipes_[o.a]; auto& b=pipes_[o.b];
|
|
const int ia=end(a,o.left_a), ib=end(b,o.left_b);
|
|
if(a.solver==PipeSolver::Muscl && b.solver==PipeSolver::Muscl &&
|
|
a.cell_area(ia)==b.cell_area(ib) && a.cell_length(ia)==b.cell_length(ib)) {
|
|
auto sa=a.state(ia), sb=b.state(ib);
|
|
const float sign=o.left_a==o.left_b?-1.f:1.f;
|
|
sa.velocity*=sign; sb.velocity*=sign;
|
|
a.set_acoustic_neighbor(o.left_a,sb); b.set_acoustic_neighbor(o.left_b,sa);
|
|
}
|
|
}
|
|
for(auto& pipe:pipes_) pipe.prepare_acoustic_reconstruction(dt);
|
|
for (auto &o : orifices_)
|
|
resolve_orifice(o, dt);
|
|
for (auto &o : opens_)
|
|
resolve_open(o, dt);
|
|
for (auto &o : connections_)
|
|
resolve_connection(o);
|
|
for(auto& node:nodes_) resolve_junction(node,dt);
|
|
}
|
|
|
|
private:
|
|
struct JunctionNode {
|
|
std::vector<JunctionFace> ends;
|
|
std::vector<GasState> states;
|
|
struct Face { double pressure{},rho{},velocity{},mass{},enthalpy{},air{},fuel{},loss{},equivalent{}; };
|
|
std::vector<Face> faces;
|
|
JunctionMethod method{JunctionMethod::Characteristic};
|
|
double turning_loss{};
|
|
};
|
|
std::vector<JunctionNode> nodes_;
|
|
// Zero-storage Rusanov star state. Integrated mass, total energy and
|
|
// species fluxes sum to zero algebraically across every connected face.
|
|
// The node has no conserved state, timestep or storage volume of its own.
|
|
void resolve_junction(JunctionNode& node,float dt) {
|
|
if(node.method==JunctionMethod::Characteristic && characteristic_junction(node,dt)) return;
|
|
if(node.method==JunctionMethod::Characteristic) ++junction_fallbacks;
|
|
legacy_junction(node);
|
|
}
|
|
void legacy_junction(JunctionNode& node) {
|
|
double speed=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto e=node.ends[k];auto& p=pipes_[e.pipe];auto s=p.state(end(p,e.left));
|
|
s.velocity*=e.left?-1:1;node.states[k]=s;
|
|
speed=std::max(speed,2*(std::sqrt(gamma*s.pressure/s.rho)+std::abs(double(s.velocity))));
|
|
}
|
|
double area_sum=0,rho=0,energy=0,air=0,fuel=0,equivalent=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto e=node.ends[k];auto& p=pipes_[e.pipe];const auto s=node.states[k];
|
|
const double a=p.cell_area(end(p,e.left)), u=s.velocity;
|
|
const double en=s.pressure/pipe_gm1+.5*s.rho*u*u;
|
|
area_sum+=a;rho+=a*s.rho*(1+u/speed);
|
|
energy+=a*(en+(en+s.pressure)*u/speed);
|
|
air+=a*s.rho*s.air*(1+u/speed);fuel+=a*s.rho*s.fuel*(1+u/speed);
|
|
equivalent+=a*s.rho*s.fuel_equivalent*(1+u/speed);
|
|
}
|
|
rho/=area_sum;energy/=area_sum;air/=area_sum;fuel/=area_sum;
|
|
equivalent/=area_sum;
|
|
const double pressure=energy*pipe_gm1;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto e=node.ends[k];auto& p=pipes_[e.pipe];const auto s=node.states[k];
|
|
const double a=p.cell_area(end(p,e.left)), u=s.velocity;
|
|
const double en=s.pressure/pipe_gm1+.5*s.rho*u*u;
|
|
const double orientation=e.left?-1:1;
|
|
const double m=.5*a*(s.rho*u-speed*(rho-s.rho));
|
|
const double heat=.5*a*((en+s.pressure)*u-speed*(energy-en));
|
|
const double y=.5*a*(s.rho*s.air*u-speed*(air-s.rho*s.air));
|
|
const double f=.5*a*(s.rho*s.fuel*u-speed*(fuel-s.rho*s.fuel));
|
|
const double momentum=.5*a*(s.rho*u*u+s.pressure+pressure+speed*s.rho*u);
|
|
const double marker=.5*a*(s.rho*s.fuel_equivalent*u-speed*(equivalent-s.rho*s.fuel_equivalent));
|
|
p.set_boundary_flux(e.left,{float(orientation*m),float(momentum),float(orientation*heat),float(orientation*y),float(orientation*f),float(orientation*marker)});
|
|
}
|
|
}
|
|
// Nonlinear outgoing acoustic-wave curves (shock or rarefaction), with
|
|
// common node pressure, conservative stagnation-enthalpy/species mixing,
|
|
// and outgoing contact density set from the mixed enthalpy. This is a
|
|
// compact subsonic engineering closure, not a resolved 3D junction.
|
|
bool characteristic_junction(JunctionNode& node,float dt) {
|
|
double low=std::numeric_limits<double>::max(), high=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
|
|
auto state=pipe.state(end(pipe,e.left));
|
|
state.velocity*=e.left?-1:1; // positive INTO the node
|
|
const double c=std::sqrt(gamma*state.pressure/state.rho);
|
|
if(std::abs(state.velocity)>=.95*c) return false;
|
|
node.states[k]=state;
|
|
low=std::min(low,double(state.pressure)); high=std::max(high,double(state.pressure));
|
|
node.faces[k].loss=0;
|
|
}
|
|
const auto wave_curve=[](const GasState& q,double p) {
|
|
if(p>q.pressure) {
|
|
const double a=2/((gamma+1)*q.rho), b=(gamma-1)/(gamma+1)*q.pressure;
|
|
return (p-q.pressure)*std::sqrt(a/(p+b));
|
|
}
|
|
const double c=std::sqrt(gamma*q.pressure/q.rho);
|
|
return 2*c/(gamma-1)*(std::pow(p/q.pressure,(gamma-1)/(2*gamma))-1);
|
|
};
|
|
const auto density=[](const GasState& q,double p) {
|
|
const double ratio=p/q.pressure, b=(gamma-1)/(gamma+1);
|
|
return p>q.pressure ? q.rho*(ratio+b)/(b*ratio+1) : q.rho*std::pow(ratio,1/gamma);
|
|
};
|
|
// Each pressure evaluation writes scratch face states only, never pipes.
|
|
const auto residual=[&](double pressure) {
|
|
double incoming=0, heat=0, air=0, fuel=0,equivalent=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
const auto q=node.states[k]; auto& f=node.faces[k];
|
|
double p=pressure;
|
|
for(int iteration=0;iteration<(f.loss>0?8:1);++iteration) {
|
|
const double u=q.velocity-wave_curve(q,p), rho=density(q,p);
|
|
const double target=std::clamp(pressure+std::copysign(.5*f.loss*rho*u*u,u),pressure*.05,pressure*20);
|
|
p=.5*(p+target);
|
|
}
|
|
f.pressure=p; f.rho=density(q,p); f.velocity=q.velocity-wave_curve(q,p);
|
|
const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
|
|
f.mass=pipe.cell_area(end(pipe,e.left))*f.rho*f.velocity;
|
|
f.enthalpy=gamma/(gamma-1)*p/f.rho+.5*f.velocity*f.velocity;
|
|
f.air=q.air; f.fuel=q.fuel;
|
|
f.equivalent=q.fuel_equivalent;if(f.mass>0)equivalent+=f.mass*f.equivalent;
|
|
if(f.mass>0) { incoming+=f.mass; heat+=f.mass*f.enthalpy; air+=f.mass*q.air; fuel+=f.mass*q.fuel; }
|
|
}
|
|
double balance=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto& f=node.faces[k];
|
|
if(f.mass<0 && incoming>1e-20) {
|
|
const double h=heat/incoming, thermal=h-.5*f.velocity*f.velocity;
|
|
if(!(thermal>0)) return std::numeric_limits<double>::quiet_NaN();
|
|
f.rho=gamma/(gamma-1)*f.pressure/thermal;
|
|
const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
|
|
f.mass=pipe.cell_area(end(pipe,e.left))*f.rho*f.velocity;
|
|
f.enthalpy=h; f.air=air/incoming; f.fuel=fuel/incoming;
|
|
f.equivalent=equivalent/incoming;
|
|
}
|
|
balance+=f.mass;
|
|
}
|
|
return balance;
|
|
};
|
|
// Iterate the donor-weighted direction coefficients separately from the
|
|
// scalar pressure root; K=0 remains the ideal acoustic scattering case.
|
|
double node_pressure=.5*(low+high);
|
|
for(int geometry=0;geometry<(node.turning_loss>0?3:1);++geometry) {
|
|
double lo=low*.01, hi=high*100;
|
|
double fl=residual(lo), fh=residual(hi);
|
|
if(!std::isfinite(fl) || !std::isfinite(fh) || fl<0 || fh>0) return false;
|
|
node_pressure=std::clamp(node_pressure,lo,hi);
|
|
bool converged=false;
|
|
for(int iteration=0;iteration<40;++iteration) {
|
|
const double value=residual(node_pressure); ++junction_iterations;
|
|
if(!std::isfinite(value)) return false;
|
|
double magnitude=0; for(const auto& face:node.faces) magnitude+=std::abs(face.mass);
|
|
if(std::abs(value)<1e-14+1e-8*magnitude) { converged=true; break; }
|
|
if(value>0) lo=node_pressure; else hi=node_pressure;
|
|
const double delta=std::max(1e-3,node_pressure*1e-5);
|
|
const double derivative=(residual(node_pressure+delta)-value)/delta;
|
|
const double next=node_pressure-value/derivative;
|
|
node_pressure=std::isfinite(next) && next>lo && next<hi ? next : .5*(lo+hi);
|
|
}
|
|
if(!converged) return false;
|
|
residual(node_pressure);
|
|
if(geometry+1<(node.turning_loss>0?3:1)) {
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto& f=node.faces[k]; double weight=0, loss=0;
|
|
if(f.mass<0) for(std::size_t j=0;j<node.ends.size();++j) {
|
|
const double m=std::max(0.,node.faces[j].mass); double dot=0;
|
|
for(int axis=0;axis<3;++axis) dot+=node.ends[k].direction[axis]*node.ends[j].direction[axis];
|
|
weight+=m; loss+=m*.5*(1+std::clamp(dot,-1.,1.));
|
|
}
|
|
f.loss=weight>0?node.turning_loss*loss/weight:0;
|
|
}
|
|
}
|
|
}
|
|
double incoming=0, outgoing=0, heat=0, air=0, fuel=0, scale=1,equivalent=0;
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
const auto& f=node.faces[k];
|
|
if(!std::isfinite(f.mass) || !std::isfinite(f.enthalpy)) return false;
|
|
if(std::abs(f.velocity)>=.98*std::sqrt(gamma*f.pressure/f.rho)) return false;
|
|
if(f.mass>0) { incoming+=f.mass;heat+=f.mass*f.enthalpy;air+=f.mass*f.air;fuel+=f.mass*f.fuel; }
|
|
if(f.mass>0)equivalent+=f.mass*f.equivalent;
|
|
else outgoing-=f.mass;
|
|
}
|
|
if(incoming>1e-16 && outgoing<=1e-16) {
|
|
if(incoming>1e-11) return false;
|
|
// Residual tolerance at rest must not create a one-sided mass source.
|
|
for(auto& face:node.faces) face.mass=0;
|
|
incoming=heat=air=fuel=0;
|
|
equivalent=0;
|
|
}
|
|
// Remove root roundoff by scaling all outflows together. A shared donor
|
|
// budget scales advective fluxes together, preserving node conservation.
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
auto& f=node.faces[k];
|
|
if(f.mass<0) { f.mass*=outgoing>0?incoming/outgoing:0; f.enthalpy=incoming>0?heat/incoming:0; f.air=incoming>0?air/incoming:0; f.fuel=incoming>0?fuel/incoming:0; }
|
|
if(f.mass<0)f.equivalent=incoming>0?equivalent/incoming:0;
|
|
if(f.mass>0) {
|
|
const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];const int i=end(pipe,e.left);
|
|
scale=std::min(scale,.25*pipe.cell_mass(i)/(f.mass*dt));
|
|
scale=std::min(scale,.25*double(pipe.cells[i].energy)*pipe.cell_area(i)*pipe.cell_length(i)/(f.mass*f.enthalpy*dt));
|
|
}
|
|
}
|
|
if(scale<1) { ++limited_fluxes; ++junction_limited_fluxes; }
|
|
for(std::size_t k=0;k<node.ends.size();++k) {
|
|
const auto e=node.ends[k]; auto& pipe=pipes_[e.pipe]; const auto& f=node.faces[k];
|
|
const double orientation=e.left?-1:1, m=f.mass*scale;
|
|
const double momentum=m*f.velocity+f.pressure*pipe.cell_area(end(pipe,e.left));
|
|
pipe.set_boundary_flux(e.left,{float(orientation*m),float(momentum),float(orientation*m*f.enthalpy),float(orientation*m*f.air),float(orientation*m*f.fuel),float(orientation*m*f.equivalent)});
|
|
}
|
|
return true;
|
|
}
|
|
struct Budget {
|
|
float *air{}, *exhaust{}, *fuel{}, *energy{};
|
|
float mass{}, heat{};
|
|
};
|
|
struct Orifice {
|
|
Port *port{};
|
|
int pipe{};
|
|
bool left{};
|
|
float *area{};
|
|
float cd{1}, length{}, loss{}, mdot{};
|
|
std::size_t budget{std::numeric_limits<std::size_t>::max()};
|
|
double previous_pressure{};
|
|
};
|
|
struct Open {
|
|
int pipe;
|
|
bool left;
|
|
float pressure, area, last_mass_flow{}, last_face_pressure{}, last_volume_flow{};
|
|
bool radiating{};
|
|
float temperature{300};
|
|
OpenEndRadiation reflection{};
|
|
};
|
|
struct Connection {
|
|
int a, b;
|
|
bool left_a, left_b;
|
|
};
|
|
struct Claim {
|
|
int pipe;
|
|
bool left;
|
|
};
|
|
std::vector<Pipe1D> &pipes_;
|
|
std::vector<Orifice> orifices_;
|
|
std::vector<Open> opens_;
|
|
std::vector<Connection> connections_;
|
|
std::vector<Budget> budgets_;
|
|
std::vector<Claim> claims_;
|
|
void check_face(int pipe, bool left) const {
|
|
if (pipe < 0 || static_cast<std::size_t>(pipe) >= pipes_.size())
|
|
throw std::invalid_argument("invalid boundary pipe index");
|
|
for (auto c : claims_)
|
|
if (c.pipe == pipe && c.left == left)
|
|
throw std::invalid_argument("pipe face already connected");
|
|
}
|
|
void claim(int pipe, bool left) {
|
|
check_face(pipe, left);
|
|
claims_.push_back({pipe, left});
|
|
}
|
|
static int end(const Pipe1D &p, bool left) { return left ? 0 : static_cast<int>(p.cells.size()) - 1; }
|
|
static void ghost(Pipe1D &p, bool left, GasState s) {
|
|
if (left)
|
|
p.left_ghost = s;
|
|
else
|
|
p.right_ghost = s;
|
|
}
|
|
void resolve_orifice(Orifice &o, float dt) {
|
|
auto &p = pipes_[o.pipe];
|
|
int i = end(p, o.left);
|
|
auto s = p.state(i);
|
|
float face_area = p.cell_area(i);
|
|
if (o.area && (!std::isfinite(*o.area) || *o.area < 0))
|
|
throw std::invalid_argument("invalid dynamic orifice area");
|
|
if (o.port && (!(o.port->pressure > 0) || !(o.port->density > 0) ||
|
|
!std::isfinite(o.port->pressure) || !std::isfinite(o.port->density)))
|
|
throw std::runtime_error("nonphysical component port");
|
|
const float area = o.area ? std::min(std::max(*o.area, 0.f), face_area) : 0.f;
|
|
if (!o.port || area <= 0 || o.cd == 0) {
|
|
auto g = s;
|
|
g.velocity = -s.velocity;
|
|
ghost(p, o.left, g);
|
|
auto wall_flux =
|
|
Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), face_area);
|
|
// A stationary adiabatic wall transfers momentum, never mass or energy.
|
|
p.set_boundary_flux(o.left, {0, wall_flux.p, 0, 0});
|
|
if (o.port)
|
|
o.port->mass_flow = 0;
|
|
o.mdot = 0;
|
|
return;
|
|
}
|
|
auto &v = *o.port;
|
|
if (std::abs(v.gamma_value - gamma) > 1e-6f || std::abs(v.gas_constant_value - gas_constant) > .001f)
|
|
throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
|
|
const float dp = v.pressure - s.pressure;
|
|
const double g = gamma;
|
|
const double critical = std::pow(2.0 / (g + 1), g / (g - 1));
|
|
const double c = std::sqrt(g * s.pressure / s.rho);
|
|
const double ui = s.velocity * (o.left ? 1 : -1);
|
|
const auto face_density = [&](double pressure) {
|
|
return s.rho * std::pow(pressure / s.pressure, 1 / g);
|
|
};
|
|
const auto nozzle_flow = [&](double pressure) {
|
|
const bool from = v.pressure >= pressure;
|
|
const double up = from ? v.pressure : pressure;
|
|
const double down = from ? pressure : v.pressure;
|
|
const double density = from ? v.density : face_density(pressure);
|
|
const double ratio = std::max(down / up, critical);
|
|
const double r = std::pow(ratio, 1 / g);
|
|
const double speed2 = 2 * g / (g - 1) * up / density * (1 - std::pow(ratio, (g - 1) / g));
|
|
return (from ? 1 : -1) * density * r * std::sqrt(std::max(0.0, speed2)) * area * o.cd;
|
|
};
|
|
// Match nozzle discharge to the outgoing duct characteristic. Imposing
|
|
// sqrt(dp) against the *cell-centre* pressure is an unstable explicit
|
|
// boundary near equilibrium and ignores the duct's acoustic impedance.
|
|
const double zero_sound = c - .5 * (g - 1) * ui;
|
|
if (zero_sound <= 0)
|
|
throw std::runtime_error("orifice characteristic outside subsonic model");
|
|
const double zero_pressure = s.pressure * std::pow(zero_sound / c, 2 * g / (g - 1));
|
|
double low = std::min(zero_pressure, double(v.pressure));
|
|
double high = std::max(zero_pressure, double(v.pressure));
|
|
double duct_pressure = o.length == 0 && o.previous_pressure > low && o.previous_pressure < high
|
|
? o.previous_pressure : .5 * (low + high);
|
|
const double flow_sign = v.pressure >= zero_pressure ? 1 : -1;
|
|
const double flow_scale = area * o.cd;
|
|
const double density_ratio_scale = std::pow(double(s.pressure) / v.pressure, 1 / g);
|
|
const double critical_density_ratio = std::pow(critical, 1 / g);
|
|
// Squared-flow residual avoids the nozzle sqrt singularity at dp=0.
|
|
// Safeguarded analytic Newton; bisection remains the fallback.
|
|
if (o.length == 0 && high > low)
|
|
for (int iteration = 0; iteration < 40; ++iteration) {
|
|
++nonlinear_iterations;
|
|
const double density = face_density(duct_pressure);
|
|
const double sound = std::sqrt(g * duct_pressure / density);
|
|
const double velocity = ui + 2 * (sound - c) / (g - 1);
|
|
const double q = density * p.cell_area(i) * velocity;
|
|
const double dq = p.cell_area(i) * density / (g * duct_pressure) * (velocity + sound);
|
|
const bool from = v.pressure >= duct_pressure;
|
|
const double up = from ? v.pressure : duct_pressure;
|
|
const double down = from ? duct_pressure : v.pressure;
|
|
const double upstream_density = from ? v.density : density;
|
|
const double raw_ratio = down / up;
|
|
const double ratio = std::max(raw_ratio, critical);
|
|
const double root_ratio = raw_ratio <= critical
|
|
? critical_density_ratio
|
|
: (from ? density / s.rho * density_ratio_scale
|
|
: s.rho / (density * density_ratio_scale));
|
|
const double r1 = root_ratio * root_ratio, r2 = root_ratio * ratio;
|
|
const double coefficient = flow_scale * flow_scale * 2 * g / (g - 1) * upstream_density * up;
|
|
const double q2 = coefficient * (r1 - r2);
|
|
const double dr = raw_ratio > critical ? (from ? 1 / up : -ratio / up) : 0;
|
|
const double derivative = coefficient * ((2 / g * r1 - (g + 1) / g * r2) / ratio) * dr +
|
|
(from ? 0 : q2 * (1 + 1 / g) / duct_pressure);
|
|
const double residual = flow_sign * (q * q - q2);
|
|
if (residual > 0)
|
|
high = duct_pressure;
|
|
else
|
|
low = duct_pressure;
|
|
if (high - low < std::max(1e-9, duct_pressure * 1e-12))
|
|
break;
|
|
const double slope = flow_sign * (2 * q * dq - derivative);
|
|
const double candidate = duct_pressure - residual / slope;
|
|
if (std::abs(residual) < 4e-7 * std::max(q * q + q2, 1e-24))
|
|
break;
|
|
duct_pressure = (std::isfinite(candidate) && candidate > low && candidate < high)
|
|
? candidate
|
|
: .5 * (low + high);
|
|
}
|
|
o.previous_pressure = duct_pressure;
|
|
float mdot = static_cast<float>(nozzle_flow(duct_pressure));
|
|
const bool upstream_port = mdot >= 0;
|
|
const float pup = upstream_port ? v.pressure : static_cast<float>(duct_pressure);
|
|
const float ratio =
|
|
static_cast<float>(std::max((upstream_port ? duct_pressure : v.pressure) / pup, critical));
|
|
const float throat_rho =
|
|
(upstream_port ? v.density : static_cast<float>(face_density(duct_pressure))) *
|
|
std::pow(ratio, 1 / gamma);
|
|
if (o.length > 0) {
|
|
// d(mdot)/dt = A/L dp - K mdot |mdot|/(2 rho A L).
|
|
const float donor_rho = o.mdot >= 0 ? v.density : s.rho;
|
|
// Semi-implicit resistance remains dissipative for large dt.
|
|
mdot = (o.mdot + area / o.length * dp * dt) /
|
|
(1 + dt * o.loss * std::abs(o.mdot) / (2 * donor_rho * area * o.length));
|
|
}
|
|
const bool from = mdot >= 0;
|
|
const float h = from ? gamma / (gamma - 1) * v.pressure / v.density
|
|
: gamma / (gamma - 1) * s.pressure / s.rho + .5f * s.velocity * s.velocity;
|
|
const float requested_mdot = mdot;
|
|
if (from && o.budget < budgets_.size()) {
|
|
auto &b = budgets_[o.budget];
|
|
mdot = std::min(mdot, std::max(0.f, std::min(b.mass, b.heat / std::max(h, 1.f))) / dt);
|
|
b.mass -= mdot * dt;
|
|
b.heat -= mdot * h * dt;
|
|
} else if (!from) {
|
|
const float cell_energy = p.cells[i].energy * face_area * p.cell_length(i);
|
|
mdot = std::max(mdot, -.25f * std::min(p.cell_mass(i), cell_energy / std::max(h, 1.f)) / dt);
|
|
}
|
|
if (mdot != requested_mdot)
|
|
++limited_fluxes;
|
|
o.mdot = mdot;
|
|
const float donor_rho = from ? v.density : s.rho;
|
|
const float density = o.length > 0 ? donor_rho : throat_rho;
|
|
const float velocity = mdot / (density * area);
|
|
const float face_pressure =
|
|
o.length > 0 ? .5f * (v.pressure + s.pressure) : static_cast<float>(duct_pressure);
|
|
const float direction = o.left ? 1.f : -1.f;
|
|
const float species = from ? v.air_fraction : s.air;
|
|
const float fuel = from ? v.fuel_fraction : s.fuel;
|
|
const float equivalent=from?std::max(v.fuel_equivalent_fraction,v.fuel_fraction):s.fuel_equivalent;
|
|
const float signed_mass = direction * mdot;
|
|
p.set_boundary_flux(o.left, {signed_mass,
|
|
mdot * velocity + face_pressure * area + s.pressure * (face_area - area),
|
|
signed_mass * h, signed_mass * species, signed_mass * fuel,signed_mass*equivalent});
|
|
ghost(p, o.left, {density, direction * velocity, face_pressure, species, fuel,equivalent});
|
|
v.mass_flow = -mdot;
|
|
v.enthalpy = h;
|
|
if (!from)
|
|
v.air_fraction = s.air;
|
|
if (!from)
|
|
v.fuel_fraction = s.fuel;
|
|
if(!from)v.fuel_equivalent_fraction=s.fuel_equivalent;
|
|
}
|
|
void resolve_open(Open &o, float dt) {
|
|
auto &p = pipes_[o.pipe];
|
|
int i = end(p, o.left);
|
|
auto s = p.state(i);
|
|
const float c = std::sqrt(gamma * s.pressure / s.rho), gm1 = gamma - 1;
|
|
if (o.left ? s.velocity >= c : s.velocity <= -c)
|
|
throw std::runtime_error(
|
|
"supersonic inlet needs full upstream state; pressure-only boundary is insufficient");
|
|
if (o.radiating) {
|
|
const float direction = o.left ? -1.f : 1.f;
|
|
const double outward = direction*s.velocity;
|
|
GasState face = s;
|
|
if (outward < c) {
|
|
const double reference_c = c*std::pow(double(o.pressure)/s.pressure, (gamma-1)/(2*gamma));
|
|
const double incident = outward + 2*(c-reference_c)/(gamma-1);
|
|
const double radius = std::sqrt(o.area/pi);
|
|
const double ambient_c = std::sqrt(gamma*gas_constant*o.temperature);
|
|
const double reflected = o.reflection.process(incident, dt*ambient_c/radius);
|
|
const double sound = reference_c + (gamma-1)*.25*(incident+reflected);
|
|
if (!(sound > 0)) throw std::runtime_error("radiation characteristic outside subsonic model");
|
|
face.pressure = o.pressure*std::pow(sound/reference_c, 2*gamma/(gamma-1));
|
|
face.velocity = direction*.5*(incident-reflected);
|
|
face.rho = s.rho*std::pow(face.pressure/s.pressure, 1/gamma);
|
|
if (direction*face.velocity < 0) {
|
|
// Fresh ambient gas enters with ambient stagnation enthalpy.
|
|
const float temperature = o.temperature - .5f*face.velocity*face.velocity /
|
|
(gamma/(gamma-1)*gas_constant);
|
|
if (!(temperature > 0)) throw std::runtime_error("ambient inflow exceeds stagnation enthalpy");
|
|
face.rho = face.pressure/(gas_constant*temperature);
|
|
face.air = 1; face.fuel = 0;
|
|
face.fuel_equivalent=0;
|
|
}
|
|
}
|
|
ghost(p, o.left, face);
|
|
const float mass_flux = face.rho*face.velocity*o.area;
|
|
const float enthalpy = gamma/(gamma-1)*face.pressure/face.rho + .5f*face.velocity*face.velocity;
|
|
p.set_boundary_flux(o.left, {mass_flux, mass_flux*face.velocity+face.pressure*o.area,
|
|
mass_flux*enthalpy, mass_flux*face.air, mass_flux*face.fuel,mass_flux*face.fuel_equivalent});
|
|
o.last_volume_flow = direction*face.velocity*o.area;
|
|
o.last_mass_flow = direction*mass_flux;
|
|
o.last_face_pressure = face.pressure;
|
|
return;
|
|
}
|
|
const float entropy = s.pressure / std::pow(s.rho, gamma);
|
|
const float rho = std::pow(o.pressure / entropy, 1 / gamma);
|
|
const float cb = std::sqrt(gamma * o.pressure / rho);
|
|
const float u =
|
|
o.left ? s.velocity - 2 * c / gm1 + 2 * cb / gm1 : s.velocity + 2 * c / gm1 - 2 * cb / gm1;
|
|
const bool out = o.left ? s.velocity <= -c : s.velocity >= c;
|
|
const bool ambient_inflow = o.left ? u >= 0 : u <= 0;
|
|
GasState g = out ? s : GasState{rho, u, o.pressure, ambient_inflow ? 1.f : s.air,
|
|
ambient_inflow ? 0.f : s.fuel,ambient_inflow?0.f:s.fuel_equivalent};
|
|
ghost(p, o.left, g);
|
|
const float a = std::min(o.area, p.cell_area(i));
|
|
auto f = Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), a);
|
|
f.p += s.pressure * (p.cell_area(i) - a);
|
|
p.set_boundary_flux(o.left, f);
|
|
o.last_mass_flow = f.m * (o.left ? -1 : 1);
|
|
o.last_volume_flow = o.last_mass_flow / (ambient_inflow ? rho : s.rho);
|
|
o.last_face_pressure = f.p / p.cell_area(i);
|
|
}
|
|
void resolve_connection(const Connection &o) {
|
|
auto &a = pipes_[o.a];
|
|
auto &b = pipes_[o.b];
|
|
int ia = end(a, o.left_a), ib = end(b, o.left_b);
|
|
auto sa = a.state(ia), sb = b.state(ib);
|
|
if(a.solver==PipeSolver::Muscl && b.solver==PipeSolver::Muscl) {
|
|
sa=a.acoustic_edge(o.left_a); sb=b.acoustic_edge(o.left_b);
|
|
}
|
|
const float da = o.left_a ? -1.f : 1.f, db = o.left_b ? 1.f : -1.f;
|
|
sa.velocity *= da;
|
|
sb.velocity *= db;
|
|
const float area = std::min(a.cell_area(ia), b.cell_area(ib));
|
|
auto f = Pipe1D::scaled(Pipe1D::numerical_flux(sa, sb, a.solver == b.solver ? a.solver : PipeSolver::Rusanov), area);
|
|
a.set_boundary_flux(o.left_a,
|
|
{da * f.m, f.p + sa.pressure * (a.cell_area(ia) - area), da * f.e,
|
|
da * f.y, da * f.fuel,da*f.fuel_equivalent});
|
|
b.set_boundary_flux(o.left_b,
|
|
{db * f.m, f.p + sb.pressure * (b.cell_area(ib) - area), db * f.e,
|
|
db * f.y, db * f.fuel,db*f.fuel_equivalent});
|
|
}
|
|
};
|
|
|
|
} // namespace fluidsim::port
|