fastsim_core/vehicle/powertrain/

fuel_converter.rs

use super::utils::ScalingMethods;
use super::*;
use crate::prelude::*;
use crate::utils::interp::InterpolatorMutMethods;
use std::f64::consts::PI;

// TODO: think about how to incorporate life modeling for Fuel Cells and other tech

#[serde_api]
#[derive(Deserialize, Serialize, Debug, Clone, PartialEq, StateMethods, SetCumulative)]
/// Struct for modeling [FuelConverter] (e.g. engine, fuel cell.) thermal plant
#[non_exhaustive]
#[serde(deny_unknown_fields)]
#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
pub struct FuelConverter {
    /// [Self] Thermal plant, including thermal management controls
    #[serde(default, skip_serializing_if = "FuelConverterThermalOption::is_none")]
    #[has_state]
    pub thrml: FuelConverterThermalOption,
    /// [Self] mass
    #[serde(default)]
    pub(in super::super) mass: Option<si::Mass>,
    /// FuelConverter specific power
    pub(in super::super) specific_pwr: Option<si::SpecificPower>,
    /// max rated brake output power
    pub pwr_out_max: si::Power,
    /// starting/baseline transient power limit
    #[serde(default)]
    pub pwr_out_max_init: si::Power,
    // TODO: consider a ramp down rate, which may be needed for fuel cells
    /// lag time for ramp up
    pub pwr_ramp_lag: si::Time,
    /// interpolator for calculating [Self] efficiency as a function of output power
    pub eff_interp_from_pwr_out: InterpolatorEnumOwned<f64>,
    /// power at which peak efficiency occurs
    #[serde(skip)]
    pub(crate) pwr_for_peak_eff: si::Power,
    /// idle fuel power to overcome internal friction (not including aux load) \[W\]
    pub pwr_idle_fuel: si::Power,
    /// time step interval between saves. 1 is a good option. If None, no saving occurs.
    pub save_interval: Option<usize>,
    /// struct for tracking current state
    #[serde(default)]
    pub state: FuelConverterState,
    /// Custom vector of [Self::state]
    #[serde(
        default,
        skip_serializing_if = "FuelConverterStateHistoryVec::is_empty"
    )]
    pub history: FuelConverterStateHistoryVec,
}

#[pyo3_api]
impl FuelConverter {
    // optional, custom, struct-specific pymethods
    #[getter("eff_max")]
    fn get_eff_max_py(&self) -> PyResult<f64> {
        Ok(*self.get_eff_max()?)
    }

    #[setter("__eff_max")]
    fn set_eff_max_py(&mut self, eff_max: f64) -> PyResult<()> {
        Ok(self.set_eff_max(eff_max, None)?)
    }

    #[getter("eff_min")]
    fn get_eff_min_py(&self) -> PyResult<f64> {
        Ok(*self.get_eff_min()?)
    }

    #[setter("__eff_min")]
    fn set_eff_min_py(&mut self, eff_min: f64) -> PyResult<()> {
        Ok(self.set_eff_min(eff_min, None)?)
    }

    #[setter("__eff_range")]
    fn set_eff_range_py(&mut self, eff_range: f64) -> PyResult<()> {
        self.set_eff_range(eff_range)?;
        Ok(())
    }

    // TODO: handle `side_effects` and uncomment
    // #[setter("__mass_kg")]
    // fn set_mass_py(&mut self, mass_kg: Option<f64>) -> anyhow::Result<()> {
    //     self.set_mass(mass_kg.map(|m| m * uc::KG))?;
    //     Ok(())
    // }

    #[getter("mass_kg")]
    fn get_mass_py(&self) -> PyResult<Option<f64>> {
        Ok(self.mass()?.map(|m| m.get::<si::kilogram>()))
    }

    #[getter]
    fn get_specific_pwr_kw_per_kg(&self) -> Option<f64> {
        self.specific_pwr
            .map(|x| x.get::<si::kilowatt_per_kilogram>())
    }
}

impl SerdeAPI for FuelConverter {}
impl Init for FuelConverter {
    fn init(&mut self) -> Result<(), Error> {
        let _ = self
            .mass()
            .map_err(|err| Error::InitError(format_dbg!(err)))?;
        self.thrml.init()?;
        self.state
            .init()
            .map_err(|err| Error::InitError(format_dbg!(err)))?;
        let eff_max = self
            .get_eff_max()
            .map_err(|err| Error::InitError(format_dbg!(err)))?;
        self.pwr_for_peak_eff = match &self.eff_interp_from_pwr_out {
            InterpolatorEnum::Interp1D(interp) => *interp.data.grid[0]
                .get(
                    interp
                        .data
                        .values
                        .iter()
                        .position(|eff| eff == eff_max)
                        .ok_or_else(|| Error::InitError(format_dbg!()))?,
                )
                .ok_or_else(|| Error::InitError(format_dbg!()))?,
            _ => {
                return Err(Error::InitError(format_dbg!(
                    "Only 1-D interpolators are supported"
                )))
            }
        } * self.pwr_out_max;
        Ok(())
    }
}
impl HistoryMethods for FuelConverter {
    fn save_interval(&self) -> anyhow::Result<Option<usize>> {
        Ok(self.save_interval)
    }
    fn set_save_interval(&mut self, save_interval: Option<usize>) -> anyhow::Result<()> {
        self.save_interval = save_interval;
        self.thrml.set_save_interval(save_interval)?;
        Ok(())
    }
    fn clear(&mut self) {
        self.history.clear();
        self.thrml.clear();
    }
}

impl Mass for FuelConverter {
    fn mass(&self) -> anyhow::Result<Option<si::Mass>> {
        let derived_mass = self
            .derived_mass()
            .with_context(|| anyhow!(format_dbg!()))?;
        if let (Some(derived_mass), Some(set_mass)) = (derived_mass, self.mass) {
            ensure!(
                utils::almost_eq_uom(&set_mass, &derived_mass, None),
                format!(
                    "{}",
                    format_dbg!(utils::almost_eq_uom(&set_mass, &derived_mass, None)),
                )
            );
        }
        Ok(self.mass)
    }

    fn set_mass(
        &mut self,
        new_mass: Option<si::Mass>,
        side_effect: MassSideEffect,
    ) -> anyhow::Result<()> {
        let derived_mass = self
            .derived_mass()
            .with_context(|| anyhow!(format_dbg!()))?;
        if let (Some(derived_mass), Some(new_mass)) = (derived_mass, new_mass) {
            if derived_mass != new_mass {
                match side_effect {
                    MassSideEffect::Extensive => {
                        self.pwr_out_max = self.specific_pwr.ok_or_else(|| {
                            anyhow!(
                                "{}\nExpected `self.specific_pwr` to be `Some`.",
                                format_dbg!()
                            )
                        })? * new_mass;
                    }
                    MassSideEffect::Intensive => {
                        self.specific_pwr = Some(self.pwr_out_max / new_mass);
                    }
                    MassSideEffect::None => {
                        self.specific_pwr = None;
                    }
                }
            }
        } else if new_mass.is_none() {
            self.specific_pwr = None;
        }
        self.mass = new_mass;
        Ok(())
    }

    fn derived_mass(&self) -> anyhow::Result<Option<si::Mass>> {
        Ok(self
            .specific_pwr
            .map(|specific_pwr| self.pwr_out_max / specific_pwr))
    }

    fn expunge_mass_fields(&mut self) {
        self.mass = None;
        self.specific_pwr = None;
    }
}

// non-py methods
impl FuelConverter {
    /// Sets maximum possible total power [FuelConverter]
    /// can produce.
    /// # Arguments
    /// - `dt`: simulation time step size
    pub fn set_curr_pwr_out_max(&mut self, dt: si::Time) -> anyhow::Result<()> {
        if self.pwr_out_max_init == si::Power::ZERO {
            // TODO: think about how to initialize power
            self.pwr_out_max_init = self.pwr_out_max / 10.
        };
        let pwr_out_max = (*self.state.pwr_prop.get_stale(|| format_dbg!())?
            + *self.state.pwr_aux.get_stale(|| format_dbg!())?
            + self.pwr_out_max / self.pwr_ramp_lag * dt)
            .min(self.pwr_out_max)
            .max(self.pwr_out_max_init);
        self.state
            .pwr_out_max
            .update(pwr_out_max, || format_dbg!())?;
        Ok(())
    }

    /// Sets maximum possible propulsion-related power [FuelConverter]
    /// can produce, accounting for any aux-related power required.
    /// # Arguments
    /// - `pwr_aux`: aux-related power required from this component
    pub fn set_curr_pwr_prop_max(&mut self, pwr_aux: si::Power) -> anyhow::Result<()> {
        ensure!(
            pwr_aux >= si::Power::ZERO,
            format!(
                "{}\n`pwr_aux` must be >= 0",
                format_dbg!(pwr_aux >= si::Power::ZERO),
            )
        );
        self.state.pwr_aux.update(pwr_aux, || format_dbg!())?;
        self.state.pwr_prop_max.update(
            *self.state.pwr_out_max.get_fresh(|| format_dbg!())? - pwr_aux,
            || format_dbg!(),
        )?;
        Ok(())
    }

    /// Solves for this powertrain system/component efficiency and sets/returns power output values.
    /// # Arguments
    /// - `pwr_out_req`: tractive power output required to achieve presribed speed
    /// - `fc_on`: whether component is actively running
    /// - `dt`: simulation time step size
    pub fn solve(
        &mut self,
        pwr_out_req: si::Power,
        fc_on: bool,
        dt: si::Time,
    ) -> anyhow::Result<()> {
        self.state.fc_on.update(fc_on, || format_dbg!())?;
        if fc_on {
            self.state.time_on.increment(dt, || format_dbg!())?;
        } else {
            self.state
                .time_on
                .update(si::Time::ZERO, || format_dbg!())?;
        }
        // NOTE: think about the possibility of engine braking, not urgent
        ensure!(
            pwr_out_req >= si::Power::ZERO,
            format!(
                "{}\n`pwr_out_req` must be >= 0",
                format_dbg!(pwr_out_req >= si::Power::ZERO),
            )
        );
        // if the engine is not on, `pwr_out_req` should be 0.0
        ensure!(
            fc_on || (pwr_out_req == si::Power::ZERO && *self.state.pwr_aux.get_fresh(|| format_dbg!())? == si::Power::ZERO),
            format!(
                "{}\nEngine is off but pwr_out_req + pwr_aux is non-zero\n`pwr_out_req`: {} kW\n`self.state.pwr_aux`: {} kW",
                format_dbg!(
                    fc_on
                        || (pwr_out_req == si::Power::ZERO
                            && *self.state.pwr_aux.get_fresh(|| format_dbg!())? == si::Power::ZERO)
                ),
               pwr_out_req.get::<si::kilowatt>(),
               self.state.pwr_aux.get_fresh(|| format_dbg!())?.get::<si::kilowatt>()
            )
        );
        self.state.pwr_prop.update(pwr_out_req, || format_dbg!())?;
        self.state.eff.update(
            if fc_on {
                uc::R
                    * self
                        .eff_interp_from_pwr_out
                        .interpolate(&[((pwr_out_req
                            + *self.state.pwr_aux.get_fresh(|| format_dbg!())?)
                            / self.pwr_out_max)
                            .get::<si::ratio>()])
                        .with_context(|| {
                            anyhow!(
                                "{}\n failed to calculate {}",
                                format_dbg!(),
                                stringify!(self.state.eff)
                            )
                        })?
            } else {
                si::Ratio::ZERO
            } * match self.thrml.temp_eff_coeff() {
                Some(tec) => *tec.get_fresh(|| format_dbg!())?,
                None => 1.0 * uc::R,
            },
            || format_dbg!(),
        )?;
        ensure!(
            (*self.state.eff.get_fresh(|| format_dbg!())? >= 0.0 * uc::R
                && *self.state.eff.get_fresh(|| format_dbg!())? <= 1.0 * uc::R),
            format!(
                "fc efficiency ({}) must be either between 0 and 1",
                self.state
                    .eff
                    .get_fresh(|| format_dbg!())?
                    .get::<si::ratio>()
            )
        );

        self.state.pwr_fuel.update(
            if *self.state.fc_on.get_fresh(|| format_dbg!())? {
                ((pwr_out_req + *self.state.pwr_aux.get_fresh(|| format_dbg!())?)
                    / *self.state.eff.get_fresh(|| format_dbg!())?)
                .max(self.pwr_idle_fuel)
            } else {
                si::Power::ZERO
            },
            || format_dbg!(),
        )?;
        self.state.pwr_loss.update(
            *self.state.pwr_fuel.get_fresh(|| format_dbg!())?
                - *self.state.pwr_prop.get_fresh(|| format_dbg!())?,
            || format_dbg!(),
        )?;

        // TODO: put this in `SetCumulative::set_custom_cumulative`
        // ensure!(
        //     self.state.energy_loss.get::<si::joule>() >= 0.0,
        //     format!(
        //         "{}\nEnergy loss must be non-negative",
        //         format_dbg!(self.state.energy_loss.get::<si::joule>() >= 0.0)
        //     )
        // );
        Ok(())
    }

    pub fn solve_thermal(
        &mut self,
        te_amb: si::Temperature,
        pwr_thrml_fc_to_cab: Option<si::Power>,
        veh_state: &mut VehicleState,
        dt: si::Time,
    ) -> anyhow::Result<()> {
        let veh_speed = *veh_state.speed_ach.get_stale(|| format_dbg!())?;
        self.thrml
            .solve_thermal(&self.state, te_amb, pwr_thrml_fc_to_cab, veh_speed, dt)
            .with_context(|| format_dbg!())
    }

    /// If thermal model is appropriately configured, returns current lumped [Self] temperature
    pub fn temperature(&self) -> Option<&TrackedState<si::Temperature>> {
        match &self.thrml {
            FuelConverterThermalOption::FuelConverterThermal(fct) => Some(&fct.state.temperature),
            FuelConverterThermalOption::None => None,
        }
    }

    /// Returns max value of [Self::eff_interp_from_pwr_out]
    pub fn get_eff_max(&self) -> anyhow::Result<&f64> {
        self.eff_interp_from_pwr_out.max()
    }

    /// Returns min value of [Self::eff_interp_from_pwr_out]
    pub fn get_eff_min(&self) -> anyhow::Result<&f64> {
        self.eff_interp_from_pwr_out.min()
    }

    /// Scales eff_interp_fwd and eff_interp_bwd by ratio of new `eff_max` per
    /// current calculated max (Note: this may change eff_min)
    pub fn set_eff_max(
        &mut self,
        eff_max: f64,
        scaling: Option<ScalingMethods>,
    ) -> anyhow::Result<()> {
        if (0.0..=1.0).contains(&eff_max) {
            self.eff_interp_from_pwr_out.set_max(eff_max, scaling)?;
        } else {
            return Err(anyhow!(
                "`eff_max` ({:.3}) must be between 0.0 and 1.0",
                eff_max,
            ));
        }
        // to update any dependent fields
        self.init().map_err(|err| anyhow!("{:?}", err))?;
        Ok(())
    }

    /// Scales eff_interp_fwd and eff_interp_bwd by ratio of new `eff_min` per
    /// current calculated min (Note: this may change eff_max)
    pub fn set_eff_min(
        &mut self,
        eff_min: f64,
        scaling: Option<ScalingMethods>,
    ) -> anyhow::Result<()> {
        self.eff_interp_from_pwr_out.set_min(eff_min, scaling)
    }

    /// Scales values of `eff_interp_fwd.f_x` and `eff_interp_bwd.f_x` without
    /// changing max such that max - min is equal to new range.  Will change max
    /// if needed to ensure no values are less than zero.
    pub fn set_eff_range(&mut self, eff_range: f64) -> anyhow::Result<()> {
        if eff_range <= 1.0 && eff_range >= 0. {
            self.eff_interp_from_pwr_out.set_range(eff_range)
        } else {
            Err(anyhow!(format!(
                "`eff_range` ({:.3}) must be between 0.0 and 1.0",
                eff_range,
            )))
        }
    }

    pub fn fc_thrml_state_mut(&mut self) -> Option<&mut FuelConverterThermalState> {
        match &mut self.thrml {
            FuelConverterThermalOption::FuelConverterThermal(fct) => Some(&mut fct.state),
            FuelConverterThermalOption::None => None,
        }
    }
}

#[serde_api]
#[derive(
    Clone,
    Debug,
    Default,
    Deserialize,
    Serialize,
    PartialEq,
    HistoryVec,
    StateMethods,
    SetCumulative,
)]
#[non_exhaustive]
#[serde(default)]
#[serde(deny_unknown_fields)]
#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
pub struct FuelConverterState {
    /// time step index
    pub i: TrackedState<usize>,
    /// max total output power fc can produce at current time
    pub pwr_out_max: TrackedState<si::Power>,
    /// max propulsion power fc can produce at current time
    pub pwr_prop_max: TrackedState<si::Power>,
    /// efficiency evaluated at current demand
    pub eff: TrackedState<si::Ratio>,
    /// instantaneous power going to drivetrain, not including aux
    pub pwr_prop: TrackedState<si::Power>,
    /// integral of [Self::pwr_prop]
    pub energy_prop: TrackedState<si::Energy>,
    /// power going to auxiliaries
    pub pwr_aux: TrackedState<si::Power>,
    /// Integral of [Self::pwr_aux]
    pub energy_aux: TrackedState<si::Energy>,
    /// instantaneous fuel power flow
    pub pwr_fuel: TrackedState<si::Power>,
    /// Integral of [Self::pwr_fuel]
    pub energy_fuel: TrackedState<si::Energy>,
    /// loss power, including idle
    pub pwr_loss: TrackedState<si::Power>,
    /// Integral of [Self::pwr_loss]
    pub energy_loss: TrackedState<si::Energy>,
    /// If true, engine is on, and if false, off (no idle)
    pub fc_on: TrackedState<bool>,
    /// Time the engine has been on
    pub time_on: TrackedState<si::Time>,
}

#[pyo3_api]
impl FuelConverterState {}
impl SerdeAPI for FuelConverterState {}
impl Init for FuelConverterState {}

/// Options for handling [FuelConverter] thermal model
#[derive(
    Clone, Default, Debug, Serialize, Deserialize, PartialEq, IsVariant, derive_more::From, TryInto,
)]
pub enum FuelConverterThermalOption {
    /// Basic thermal plant for [FuelConverter]
    FuelConverterThermal(Box<FuelConverterThermal>),
    /// no thermal plant for [FuelConverter]
    #[default]
    None,
}

impl StateMethods for FuelConverterThermalOption {}

impl SaveState for FuelConverterThermalOption {
    fn save_state<F: Fn() -> String>(&mut self, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => fct.save_state(loc)?,
            Self::None => {}
        }
        Ok(())
    }
}
impl TrackedStateMethods for FuelConverterThermalOption {
    fn check_and_reset<F: Fn() -> String>(&mut self, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => {
                fct.check_and_reset(|| format!("{}\n{}", loc(), format_dbg!()))?
            }
            Self::None => {}
        }
        Ok(())
    }

    fn mark_fresh<F: Fn() -> String>(&mut self, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => {
                fct.mark_fresh(|| format!("{}\n{}", loc(), format_dbg!()))?
            }
            Self::None => {}
        }
        Ok(())
    }
}
impl Step for FuelConverterThermalOption {
    fn step<F: Fn() -> String>(&mut self, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => fct.step(|| format!("{}\n{}", loc(), format_dbg!())),
            Self::None => Ok(()),
        }
    }

    fn reset_step<F: Fn() -> String>(&mut self, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => {
                fct.reset_step(|| format!("{}\n{}", loc(), format_dbg!()))
            }
            Self::None => Ok(()),
        }
    }
}
impl Init for FuelConverterThermalOption {
    fn init(&mut self) -> Result<(), Error> {
        match self {
            Self::FuelConverterThermal(fct) => fct.init()?,
            Self::None => {}
        }
        Ok(())
    }
}
impl SerdeAPI for FuelConverterThermalOption {}
impl SetCumulative for FuelConverterThermalOption {
    fn set_cumulative<F: Fn() -> String>(&mut self, dt: si::Time, loc: F) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => {
                fct.set_cumulative(dt, || format!("{}\n{}", loc(), format_dbg!()))?
            }
            Self::None => {}
        }
        Ok(())
    }
}
impl HistoryMethods for FuelConverterThermalOption {
    fn save_interval(&self) -> anyhow::Result<Option<usize>> {
        match self {
            FuelConverterThermalOption::FuelConverterThermal(fct) => fct.save_interval(),
            FuelConverterThermalOption::None => Ok(None),
        }
    }
    fn set_save_interval(&mut self, save_interval: Option<usize>) -> anyhow::Result<()> {
        match self {
            FuelConverterThermalOption::FuelConverterThermal(fct) => {
                fct.set_save_interval(save_interval)
            }
            FuelConverterThermalOption::None => Ok(()),
        }
    }
    fn clear(&mut self) {
        match self {
            FuelConverterThermalOption::FuelConverterThermal(fct) => {
                fct.clear();
            }
            FuelConverterThermalOption::None => {}
        }
    }
}
impl FuelConverterThermalOption {
    /// Solve change in temperature and other thermal effects
    /// # Arguments
    /// - `fc_state`: [FuelConverter] state
    /// - `te_amb`: ambient temperature
    /// - `pwr_thrml_fc_to_cab`: heat demand from [Vehicle::hvac] system -- zero if `None` is passed
    /// - `veh_speed`: current achieved speed
    fn solve_thermal(
        &mut self,
        fc_state: &FuelConverterState,
        te_amb: si::Temperature,
        pwr_thrml_fc_to_cab: Option<si::Power>,
        veh_speed: si::Velocity,
        dt: si::Time,
    ) -> anyhow::Result<()> {
        match self {
            Self::FuelConverterThermal(fct) => fct
                .solve(
                    fc_state,
                    te_amb,
                    pwr_thrml_fc_to_cab.unwrap_or_default(),
                    veh_speed,
                    dt,
                )
                .with_context(|| format_dbg!())?,
            Self::None => {
                ensure!(
                    pwr_thrml_fc_to_cab.is_none(),
                    format_dbg!(
                        "`FuelConverterThermal needs to be configured to provide heat demand`"
                    )
                );
            }
        }
        Ok(())
    }

    /// If appropriately configured, returns temperature-dependent efficiency coefficient
    fn temp_eff_coeff(&self) -> Option<&TrackedState<si::Ratio>> {
        match self {
            Self::FuelConverterThermal(fct) => Some(&fct.state.eff_coeff),
            Self::None => None,
        }
    }
}

#[serde_api]
#[derive(Deserialize, Serialize, Debug, Clone, PartialEq, StateMethods)]
#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
#[non_exhaustive]
#[serde(deny_unknown_fields)]
/// Struct for modeling Fuel Converter (e.g. engine, fuel cell.)
pub struct FuelConverterThermal {
    /// [FuelConverter] thermal capacitance
    pub heat_capacitance: si::HeatCapacity,
    /// parameter for engine characteristic length for heat transfer calcs
    pub length_for_convection: si::Length,
    /// parameter for heat transfer coeff from [FuelConverter] to ambient during vehicle stop
    pub htc_to_amb_stop: si::HeatTransferCoeff,

    /// Heat transfer coefficient between adiabatic flame temperature and [FuelConverterThermal] temperature
    pub conductance_from_comb: si::ThermalConductance,
    /// Max coefficient for fraction of combustion heat that goes to [FuelConverter]
    /// (engine) thermal mass. Remainder goes to environment (e.g. via tailpipe).
    pub max_frac_from_comb: si::Ratio,
    /// parameter for temperature at which thermostat starts to open
    pub tstat_te_sto: Option<si::Temperature>,
    /// temperature delta over which thermostat is partially open
    pub tstat_te_delta: Option<si::TemperatureInterval>,
    #[serde(default = "tstat_interp_default")]
    pub tstat_interp: Interp1DOwned<f64, strategy::Linear>,
    /// Radiator effectiveness -- ratio of active heat rejection from
    /// radiator to passive heat rejection, always greater than 1
    pub radiator_effectiveness: si::Ratio,
    /// Model for [FuelConverter] dependence on efficiency
    pub fc_eff_model: FCTempEffModel,
    /// struct for tracking current state
    #[serde(default)]
    pub state: FuelConverterThermalState,
    /// Custom vector of [Self::state]
    #[serde(
        default,
        skip_serializing_if = "FuelConverterThermalStateHistoryVec::is_empty"
    )]
    pub history: FuelConverterThermalStateHistoryVec,
    pub save_interval: Option<usize>,
}

#[pyo3_api]
impl FuelConverterThermal {
    #[staticmethod]
    #[pyo3(name = "default")]
    fn default_py() -> Self {
        Default::default()
    }
}

impl HistoryMethods for FuelConverterThermal {
    fn save_interval(&self) -> anyhow::Result<Option<usize>> {
        Ok(self.save_interval)
    }
    fn set_save_interval(&mut self, save_interval: Option<usize>) -> anyhow::Result<()> {
        self.save_interval = save_interval;
        Ok(())
    }
    fn clear(&mut self) {
        self.history.clear();
    }
}

/// Dummy interpolator that will be overridden in [FuelConverterThermal::init]
fn tstat_interp_default() -> Interp1DOwned<f64, strategy::Linear> {
    Interp1D::new(
        array![85.0, 90.0],
        array![0.0, 1.0],
        strategy::Linear,
        Extrapolate::Clamp,
    )
    .unwrap()
}

lazy_static! {
    /// gasoline stoichiometric air-fuel ratio https://en.wikipedia.org/wiki/Air%E2%80%93fuel_ratio
    pub static ref AFR_STOICH_GASOLINE: si::Ratio = uc::R * 14.7;
    /// gasoline density in https://inchem.org/documents/icsc/icsc/eics1400.htm
    /// This is reasonably constant with respect to temperature and pressure
    pub static ref GASOLINE_DENSITY: si::MassDensity = 0.75 * uc::KG / uc::L;
    /// TODO: find a source for this value
    pub static ref GASOLINE_LHV: si::SpecificEnergy = 33.7 * uc::KWH / uc::GALLON / *GASOLINE_DENSITY;
    pub static ref TE_ADIABATIC_STD: si::Temperature = Air::get_te_from_u(
            Air::get_specific_energy(*TE_STD_AIR).with_context(|| format_dbg!()).unwrap()
                + (Octane::get_specific_energy(*TE_STD_AIR).with_context(|| format_dbg!()).unwrap()
                    + *GASOLINE_LHV)
                    / *AFR_STOICH_GASOLINE,
        )
        .with_context(|| format_dbg!()).unwrap_or_else(|_| panic!("{}\nFailed to calculate adiabatic flame temp for gasoline", format_dbg!()));
}

impl FuelConverterThermal {
    /// Solve change in temperature and other thermal effects
    /// # Arguments
    /// - `fc_state`: [FuelConverter] state
    /// - `te_amb`: ambient temperature
    /// - `pwr_thrml_fc_to_cab`: heat demand from [Vehicle::hvac] system
    /// - `veh_speed`: current achieved speed
    /// - `dt`: simulation time step size
    fn solve(
        &mut self,
        fc_state: &FuelConverterState,
        te_amb: si::Temperature,
        pwr_thrml_fc_to_cab: si::Power,
        veh_speed: si::Velocity,
        dt: si::Time,
    ) -> anyhow::Result<()> {
        self.state
            .pwr_thrml_fc_to_cab
            .update(pwr_thrml_fc_to_cab, || format_dbg!())?;
        // film temperature for external convection calculations
        let te_air_film: si::Temperature = 0.5
            * (self
                .state
                .temperature
                .get_stale(|| format_dbg!())?
                .get::<si::kelvin_abs>()
                + te_amb.get::<si::kelvin_abs>())
            * uc::KELVIN;
        // Reynolds number = density * speed * diameter / dynamic viscosity
        // NOTE: might be good to pipe in elevation
        let fc_air_film_re =
            Air::get_density(Some(te_air_film), None) * veh_speed * self.length_for_convection
                / Air::get_dyn_visc(te_air_film).with_context(|| format_dbg!())?;

        // calculate heat transfer coeff. from engine to ambient [W / (m ** 2 * K)]
        self.state.htc_to_amb.update(
            if veh_speed < 1.0 * uc::MPS {
                // if stopped, scale based on thermostat opening and constant convection
                self.state.tstat_open_frac.update(
                    self.tstat_interp
                        .interpolate(&[self
                            .state
                            .temperature
                            .get_stale(|| format_dbg!())?
                            .get::<si::degree_celsius>()])
                        .with_context(|| format_dbg!())?,
                    || format_dbg!(),
                )?;
                (uc::R
                    + *self.state.tstat_open_frac.get_fresh(|| format_dbg!())?
                        * self.radiator_effectiveness)
                    * self.htc_to_amb_stop
            } else {
                // Calculate heat transfer coefficient for sphere,
                // from Incropera's Intro to Heat Transfer, 5th Ed., eq. 7.44
                let sphere_conv_params = get_sphere_conv_params(fc_air_film_re.get::<si::ratio>());
                let htc_to_amb_sphere: si::HeatTransferCoeff = sphere_conv_params.0
                    * fc_air_film_re.get::<si::ratio>().powf(sphere_conv_params.1)
                    * Air::get_pr(te_air_film)
                        .with_context(|| format_dbg!())?
                        .get::<si::ratio>()
                        .powf(1.0 / 3.0)
                    * Air::get_therm_cond(te_air_film).with_context(|| format_dbg!())?
                    / self.length_for_convection;
                // if stopped, scale based on thermostat opening and constant convection
                self.state.tstat_open_frac.update(
                    self.tstat_interp
                        .interpolate(&[self
                            .state
                            .temperature
                            .get_stale(|| format_dbg!())?
                            .get::<si::degree_celsius>()])
                        .with_context(|| format_dbg!())?,
                    || format_dbg!(),
                )?;
                *self.state.tstat_open_frac.get_fresh(|| format_dbg!())? * htc_to_amb_sphere
            },
            || format_dbg!(),
        )?;

        self.state.pwr_thrml_to_amb.update(
            *self.state.htc_to_amb.get_fresh(|| format_dbg!())?
                * PI
                * self.length_for_convection.powi(typenum::P2::new())
                / 4.0
                * (self
                    .state
                    .temperature
                    .get_stale(|| format_dbg!())?
                    .get::<si::degree_celsius>()
                    - te_amb.get::<si::degree_celsius>())
                * uc::KELVIN_INT,
            || format_dbg!(),
        )?;

        // let heat_to_amb = ;
        // assumes fuel/air mixture is entering combustion chamber at block temperature
        // assumes stoichiometric combustion
        self.state.te_adiabatic.update(
            Air::get_te_from_u(
                Air::get_specific_energy(*self.state.temperature.get_stale(|| format_dbg!())?)
                    .with_context(|| format_dbg!())?
                    + (Octane::get_specific_energy(*self.state.temperature.get_stale(|| format_dbg!())?)
                    .with_context(|| format_dbg!())?
                    // TODO: make config. for other fuels -- e.g. with enum for specific fuels and/or fuel properties
                    + *GASOLINE_LHV)
                        / *AFR_STOICH_GASOLINE,
            )
            .with_context(|| format_dbg!())?,
            || format_dbg!(),
        )?;
        // heat that will go both to the block and out the exhaust port
        self.state.pwr_fuel_as_heat.update(
            *fc_state.pwr_fuel.get_stale(|| format_dbg!())?
                - (*fc_state.pwr_prop.get_stale(|| format_dbg!())?
                    + *fc_state.pwr_aux.get_stale(|| format_dbg!())?),
            || format_dbg!(),
        )?;
        self.state.pwr_thrml_to_tm.update(
            (self.conductance_from_comb
                * (self
                    .state
                    .te_adiabatic
                    .get_fresh(|| format_dbg!())?
                    .get::<si::degree_celsius>()
                    - self
                        .state
                        .temperature
                        .get_stale(|| format_dbg!())?
                        .get::<si::degree_celsius>())
                * uc::KELVIN_INT)
                .min(
                    self.max_frac_from_comb
                        * *self.state.pwr_fuel_as_heat.get_fresh(|| format_dbg!())?,
                ),
            || format_dbg!(),
        )?;
        let delta_temp: si::TemperatureInterval =
            ((*self.state.pwr_thrml_to_tm.get_fresh(|| format_dbg!())?
                - *self.state.pwr_thrml_fc_to_cab.get_fresh(|| format_dbg!())?
                - *self.state.pwr_thrml_to_amb.get_fresh(|| format_dbg!())?)
                * dt)
                / self.heat_capacitance;
        // Interestingly, it seems to be ok to add a `TemperatureInterval` to a `Temperature` here
        self.state.temperature.update(
            *self.state.temperature.get_stale(|| format_dbg!())? + delta_temp,
            || format_dbg!(),
        )?;

        self.state.eff_coeff.update(
            match self.fc_eff_model {
                FCTempEffModel::Linear(FCTempEffModelLinear {
                    offset,
                    slope_per_kelvin: slope,
                    minimum,
                }) => minimum.max(
                    {
                        let calc_unbound: si::Ratio = offset
                            + slope * uc::R / uc::KELVIN
                                * *self.state.temperature.get_fresh(|| format_dbg!())?;
                        calc_unbound
                    }
                    .min(1.0 * uc::R),
                ),
                FCTempEffModel::Exponential(FCTempEffModelExponential {
                    offset,
                    lag,
                    minimum,
                }) => {
                    let dte: si::TemperatureInterval = (self
                        .state
                        .temperature
                        .get_fresh(|| format_dbg!())?
                        .get::<si::kelvin_abs>()
                        - offset.get::<si::kelvin_abs>())
                        * uc::KELVIN_INT;
                    ((1.0 - f64::exp((-dte / lag).get::<si::ratio>())) * uc::R).max(minimum)
                }
            },
            || format_dbg!(),
        )?;
        Ok(())
    }
}
impl SerdeAPI for FuelConverterThermal {}
impl SetCumulative for FuelConverterThermal {
    fn set_cumulative<F: Fn() -> String>(&mut self, dt: si::Time, loc: F) -> anyhow::Result<()> {
        self.state
            .set_cumulative(dt, || format!("{}\n{}", loc(), format_dbg!()))
    }
}
impl Init for FuelConverterThermal {
    fn init(&mut self) -> Result<(), Error> {
        self.tstat_te_sto = self
            .tstat_te_sto
            .or(Some((85. + uc::CELSIUS_TO_KELVIN) * uc::KELVIN));
        self.tstat_te_delta = self.tstat_te_delta.or(Some(5. * uc::KELVIN_INT));
        self.tstat_interp = Interp1D::new(
            array![
                self.tstat_te_sto.unwrap().get::<si::degree_celsius>(),
                self.tstat_te_sto.unwrap().get::<si::degree_celsius>()
                    + self.tstat_te_delta.unwrap().get::<si::kelvin>(),
            ],
            array![0.0, 1.0],
            strategy::Linear,
            Extrapolate::Clamp,
        )
        .map_err(|err| {
            Error::InitError(format!(
                "{}\n{}\n{}",
                err,
                format_dbg!(self.tstat_te_sto),
                format_dbg!(self.tstat_te_delta)
            ))
        })?;
        Ok(())
    }
}
impl Default for FuelConverterThermal {
    fn default() -> Self {
        let mut fct = Self {
            heat_capacitance: Default::default(),
            length_for_convection: Default::default(),
            htc_to_amb_stop: Default::default(),
            conductance_from_comb: Default::default(),
            max_frac_from_comb: Default::default(),
            tstat_te_sto: None,
            tstat_te_delta: None,
            tstat_interp: tstat_interp_default(),
            radiator_effectiveness: Default::default(),
            fc_eff_model: Default::default(),
            state: Default::default(),
            history: Default::default(),
            save_interval: Some(1),
        };
        fct.init().unwrap();
        fct
    }
}

#[serde_api]
#[derive(
    Clone, Debug, Deserialize, Serialize, PartialEq, HistoryVec, StateMethods, SetCumulative,
)]
#[serde(default)]
#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
#[serde(deny_unknown_fields)]
pub struct FuelConverterThermalState {
    /// time step index
    pub i: TrackedState<usize>,
    /// Adiabatic flame temperature assuming complete (i.e. all fuel is consumed
    /// if fuel lean or stoich or all air is consumed if fuel rich) combustion
    pub te_adiabatic: TrackedState<si::Temperature>,
    /// Current engine thermal mass temperature (lumped engine block and coolant)
    pub temperature: TrackedState<si::Temperature>,
    /// thermostat open fraction (1 = fully open, 0 = fully closed)
    pub tstat_open_frac: TrackedState<f64>,
    /// Current heat transfer coefficient from [FuelConverter] to ambient
    pub htc_to_amb: TrackedState<si::HeatTransferCoeff>,
    /// Current heat transfer power to ambient
    pub pwr_thrml_to_amb: TrackedState<si::Power>,
    /// Cumulative heat transfer energy to ambient
    pub energy_thrml_to_amb: TrackedState<si::Energy>,
    /// Efficency coefficient, used to modify [FuelConverter] effciency based on temperature
    pub eff_coeff: TrackedState<si::Ratio>,
    /// Thermal power flowing from fuel converter to cabin
    pub pwr_thrml_fc_to_cab: TrackedState<si::Power>,
    /// Cumulative thermal energy flowing from fuel converter to cabin
    pub energy_thrml_fc_to_cab: TrackedState<si::Energy>,
    /// Fuel power that is not converted to mechanical work
    pub pwr_fuel_as_heat: TrackedState<si::Power>,
    /// Cumulative fuel energy that is not converted to mechanical work
    pub energy_fuel_as_heat: TrackedState<si::Energy>,
    /// Thermal power flowing from combustion to [FuelConverter] thermal mass
    pub pwr_thrml_to_tm: TrackedState<si::Power>,
    /// Cumulative thermal energy flowing from combustion to [FuelConverter] thermal mass
    pub energy_thrml_to_tm: TrackedState<si::Energy>,
}
#[pyo3_api]
impl FuelConverterThermalState {}

impl Init for FuelConverterThermalState {}
impl SerdeAPI for FuelConverterThermalState {}
impl Default for FuelConverterThermalState {
    fn default() -> Self {
        Self {
            i: Default::default(),
            te_adiabatic: TrackedState::new(*TE_ADIABATIC_STD),
            temperature: TrackedState::new(*TE_STD_AIR),
            tstat_open_frac: Default::default(),
            htc_to_amb: Default::default(),
            eff_coeff: TrackedState::new(uc::R),
            pwr_thrml_fc_to_cab: Default::default(),
            energy_thrml_fc_to_cab: Default::default(),
            pwr_thrml_to_amb: Default::default(),
            energy_thrml_to_amb: Default::default(),
            pwr_fuel_as_heat: Default::default(),
            energy_fuel_as_heat: Default::default(),
            pwr_thrml_to_tm: Default::default(),
            energy_thrml_to_tm: Default::default(),
        }
    }
}

/// Model variants for how FC efficiency depends on temperature
#[derive(
    Debug, Clone, Deserialize, Serialize, PartialEq, IsVariant, derive_more::From, TryInto,
)]
pub enum FCTempEffModel {
    /// Linear temperature dependence
    Linear(FCTempEffModelLinear),
    /// Exponential temperature dependence
    Exponential(FCTempEffModelExponential),
}

impl Default for FCTempEffModel {
    fn default() -> Self {
        FCTempEffModel::Exponential(FCTempEffModelExponential::default())
    }
}

#[derive(Debug, Clone, Deserialize, Serialize, PartialEq)]
#[serde(deny_unknown_fields)]
pub struct FCTempEffModelLinear {
    pub offset: si::Ratio,
    /// Change in efficiency factor per change in temperature /[K/]
    pub slope_per_kelvin: f64,
    pub minimum: si::Ratio,
}

impl Default for FCTempEffModelLinear {
    fn default() -> Self {
        Self {
            offset: 0.0 * uc::R,
            slope_per_kelvin: 25.0,
            minimum: 0.2 * uc::R,
        }
    }
}

#[derive(Debug, Clone, Deserialize, Serialize, PartialEq)]
#[serde(deny_unknown_fields)]
pub struct FCTempEffModelExponential {
    /// temperature at which `fc_eta_temp_coeff` begins to grow
    pub offset: si::Temperature,
    /// exponential lag parameter [K^-1]
    pub lag: si::TemperatureInterval,
    /// minimum value that `fc_eta_temp_coeff` can take
    pub minimum: si::Ratio,
}

impl Default for FCTempEffModelExponential {
    fn default() -> Self {
        Self {
            // TODO: update after reasonable calibration
            offset: 0.0 * uc::KELVIN,
            lag: 25.0 * uc::KELVIN_INT,
            minimum: 0.2 * uc::R,
        }
    }
}