FuelType.h

/* Copyright (c) Queen's Printer for Ontario, 2020. */
/* Copyright (c) His Majesty the King in Right of Canada as represented by the Minister of Natural Resources, 2025. */
 
/* SPDX-License-Identifier: AGPL-3.0-or-later */
 
#pragma once
#include "Duff.h"
#include "FWI.h"
namespace fs
{
namespace sim
{
class SpreadInfo;
}
namespace data
{
class LogValue;
}
using sim::SpreadInfo;
using data::LogValue;
namespace fuel
{
constexpr FuelCodeSize INVALID_FUEL_CODE = 0;
// References
// Forestry Canada
// Development and Structure of the Canadian Forest Fire Behaviour Prediction System (ST-X-3)
// https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/10068.pdf
//
// Wotton, B.M., Alexander, M.E., Taylor, S.W.
// Updates and revision to the 1992 Canadian Forest Fire Behavior Prediction System (GLC-X-10)
// https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/31414.pdf
//
// Anderson, Kerry
// Incorporating Smoldering Into Fire Growth Modelling
// https://www.cfs.nrcan.gc.ca/pubwarehouse/pdfs/19950.pdf
//
// default grass fuel load (kg/m^2)
static constexpr MathSize DEFAULT_GRASS_FUEL_LOAD = 0.35;
// amount of duff to apply ffmc moisture to (cm) (1.2 cm is from Kerry's paper)
static constexpr MathSize DUFF_FFMC_DEPTH = 1.2;
[[nodiscard]] constexpr MathSize fire_intensity(const MathSize fc, const MathSize ros)
{
  return 300.0 * fc * ros;
}
class FuelType
{
public:
  [[nodiscard]] static constexpr FuelCodeSize safeCode(const FuelType* fuel)
  {
    return nullptr == fuel ? static_cast<FuelCodeSize>(INVALID_FUEL_CODE) : fuel->code();
  }
  [[nodiscard]] static constexpr const char* safeName(const FuelType* fuel)
  {
    return nullptr == fuel ? "NULL" : fuel->name();
  }
  [[nodiscard]] static constexpr MathSize criticalRos(const MathSize sfc, const MathSize csi)
  {
    return sfc > 0 ? csi / (300.0 * sfc) : 0.0;
  }
  [[nodiscard]] static constexpr bool isCrown(const MathSize csi, const MathSize sfi)
  {
    return sfi > csi;
  }
  [[nodiscard]] virtual MathSize cfl() const = 0;
  virtual ~FuelType() noexcept = default;
  constexpr FuelType(const FuelCodeSize& code,
                     const char* name,
                     const bool can_crown) noexcept
    : name_(name), can_crown_(can_crown), code_(code)
  {
  }
  FuelType(FuelType&& rhs) noexcept = delete;
  FuelType(const FuelType& rhs) noexcept = delete;
  FuelType& operator=(FuelType&& rhs) noexcept = delete;
  FuelType& operator=(const FuelType& rhs) noexcept = delete;
  [[nodiscard]] constexpr bool canCrown() const
  {
    return can_crown_;
  }
  [[nodiscard]] virtual MathSize grass_curing(const int, const wx::FwiWeather&) const
  {
    // NOTE: grass overrides this but everything else doesn't have curing
    return INVALID_CURING;
  }
 
  [[nodiscard]] virtual MathSize cbh() const = 0;
  [[nodiscard]] virtual MathSize crownFractionBurned(MathSize rss, MathSize rso) const noexcept = 0;
  [[nodiscard]] virtual MathSize probabilityPeat(MathSize mc_fraction) const noexcept = 0;
  [[nodiscard]] virtual ThresholdSize survivalProbability(
    const wx::FwiWeather& wx) const noexcept = 0;
  [[nodiscard]] virtual MathSize buiEffect(MathSize bui) const = 0;
  [[nodiscard]] virtual MathSize crownConsumption(MathSize cfb) const = 0;
  [[nodiscard]] virtual MathSize calculateRos(int nd,
                                              const wx::FwiWeather& wx,
                                              MathSize isi) const = 0;
  [[nodiscard]] virtual MathSize calculateIsf(const SpreadInfo& spread,
                                              MathSize isi) const = 0;
  [[nodiscard]] virtual MathSize surfaceFuelConsumption(
    const SpreadInfo& spread) const = 0;
  [[nodiscard]] virtual MathSize lengthToBreadth(MathSize ws) const = 0;
  [[nodiscard]] virtual MathSize finalRos(const SpreadInfo& spread,
                                          MathSize isi,
                                          MathSize cfb,
                                          MathSize rss) const = 0;
  [[nodiscard]] virtual MathSize criticalSurfaceIntensity(
    const SpreadInfo& spread) const = 0;
  [[nodiscard]] constexpr const char* name() const
  {
    return name_;
  }
  [[nodiscard]] constexpr FuelCodeSize code() const
  {
    return code_;
  }
private:
  const char* name_;
  const bool can_crown_;
  FuelCodeSize code_;
};
template <int BulkDensity, int InorganicPercent, int DuffDepth>
class FuelBase
  : public FuelType
{
public:
  ~FuelBase() override = default;
  constexpr FuelBase(const FuelCodeSize& code,
                     const char* name,
                     const bool can_crown,
                     const Duff* duff_ffmc,
                     const Duff* duff_dmc)
    : FuelType(code, name, can_crown),
      duff_ffmc_(duff_ffmc),
      duff_dmc_(duff_dmc)
  {
  }
  FuelBase(FuelBase&& rhs) noexcept = delete;
  FuelBase(const FuelBase& rhs) = delete;
  FuelBase& operator=(FuelBase&& rhs) noexcept = delete;
  FuelBase& operator=(const FuelBase& rhs) = delete;
  [[nodiscard]] MathSize crownFractionBurned(const MathSize rss,
                                             const MathSize rso) const noexcept override
  {
    // can't burn crown if it doesn't exist
    return cfl() > 0 ? max(0.0, 1.0 - exp(-0.230 * (rss - rso))) : 0.0;
  }
  [[nodiscard]] ThresholdSize probabilityPeat(const MathSize mc_fraction) const noexcept override
  {
    // Anderson table 1
    constexpr auto pb = bulkDensity();
    // Anderson table 1
    constexpr auto fi = inorganicPercent();
    constexpr auto pi = fi * pb;
    // Inorganic ratio
    constexpr auto ri = fi / (1 - fi);
    constexpr auto const_part = -19.329 + 1.7170 * ri + 23.059 * pi;
    // Anderson eq 1
    return 1 / (1 + exp(17.047 * mc_fraction / (1 - fi) + const_part));
  }
  [[nodiscard]] ThresholdSize survivalProbability(const wx::FwiWeather& wx) const noexcept
    override
  {
    // divide by 100 since we need moisture ratio
    //    IFERROR(((1 / (1 + EXP($G$43 + $I$43 *
    //            (Q$44 * $O$43 + $N$43)))) -
    //            (1 / (1 + EXP($G$43 + $I$43 * (2.5 * $O$43 + $N$43)))))
    //            / (1 / (1 + EXP($G$43 + $I$43 * $N$43))), 0)
    // HACK: use same constants for all fuels because they seem to work nicer than
    // using the ratios, but they change anyway because of the other fuel attributes
    static const auto WFfmc = 0.25;
    static const auto WDmc = 1.0;
    static const auto RatioHartford = 0.5;
    static const auto RatioFrandsen = 1.0 - RatioHartford;
    static const auto RatioAspen = 0.5;
    static const auto RatioFuel = 1.0 - RatioAspen;
    const auto mc_ffmc = wx.mcFfmc() * WFfmc + WDmc;
    static const auto McFfmcSaturated = 2.5 * WFfmc + WDmc;
    static const auto McDmc = WDmc;
    const auto prob_ffmc_peat = probabilityPeat(mc_ffmc);
    const auto prob_ffmc_peat_saturated = probabilityPeat(McFfmcSaturated);
    const auto prob_ffmc_peat_zero = probabilityPeat(McDmc);
    const auto prob_ffmc_peat_weighted = (prob_ffmc_peat - prob_ffmc_peat_saturated) / prob_ffmc_peat_zero;
    const auto prob_ffmc = duffFfmcType()->probabilityOfSurvival(mc_ffmc * 100);
    const auto prob_ffmc_saturated = duffFfmcType()->probabilityOfSurvival(
      McFfmcSaturated * 100);
    const auto prob_ffmc_zero = duffFfmcType()->probabilityOfSurvival(McDmc);
    const auto prob_ffmc_weighted = (prob_ffmc - prob_ffmc_saturated) / prob_ffmc_zero;
    const auto term_otway = exp(-3.11 + 0.12 * wx.dmc().asValue());
    const auto prob_otway = term_otway / (1 + term_otway);
    const auto mc_pct = wx.mcDmcPct() * dmcRatio() + wx.mcFfmcPct() * ffmcRatio();
    const auto prob_weight_ffmc = duffFfmcType()->probabilityOfSurvival(mc_pct);
    const auto prob_weight_ffmc_peat = probabilityPeat(mc_pct / 100);
    const auto prob_weight_dmc = duffDmcType()->probabilityOfSurvival(wx.mcDmcPct());
    const auto prob_weight_dmc_peat = probabilityPeat(wx.mcDmc());
    // chance of survival is 1 - chance of it not surviving in every fuel
    const auto tot_prob = 1 - (1 - prob_ffmc_peat_weighted) * (1 - prob_ffmc_weighted) * ((1 - prob_otway) * RatioAspen + ((1 - prob_weight_ffmc_peat) * RatioHartford + (1 - prob_weight_ffmc) * RatioFrandsen) * ((1 - prob_weight_dmc_peat) * RatioHartford + (1 - prob_weight_dmc) * RatioFrandsen) * RatioFuel);
    return tot_prob;
  }
  [[nodiscard]] static constexpr MathSize bulkDensity()
  {
    return BulkDensity / 1000.0;
  }
  [[nodiscard]] static constexpr MathSize inorganicPercent()
  {
    return InorganicPercent / 100.0;
  }
  [[nodiscard]] static constexpr MathSize duffDepth()
  {
    return DuffDepth / 10.0;
  }
  [[nodiscard]] constexpr const Duff* duffDmcType() const
  {
    return duff_dmc_;
  }
  [[nodiscard]] constexpr const Duff* duffFfmcType() const
  {
    return duff_ffmc_;
  }
  [[nodiscard]] static constexpr MathSize ffmcRatio()
  {
    return 1 - dmcRatio();
  }
  [[nodiscard]] static constexpr MathSize dmcRatio()
  {
    return (duffDepth() - DUFF_FFMC_DEPTH) / duffDepth();
  }
private:
  const Duff* duff_ffmc_;
  const Duff* duff_dmc_;
};
class InvalidFuel final
  : public FuelType
{
public:
  InvalidFuel() noexcept
    : InvalidFuel(0, nullptr)
  {
  }
  constexpr InvalidFuel(const FuelCodeSize& code, const char* name) noexcept
    : FuelType(code, name, false)
  {
  }
  ~InvalidFuel() override = default;
  InvalidFuel(const InvalidFuel& rhs) noexcept = delete;
  InvalidFuel(InvalidFuel&& rhs) noexcept = delete;
  InvalidFuel& operator=(const InvalidFuel& rhs) noexcept = delete;
  InvalidFuel& operator=(InvalidFuel&& rhs) noexcept = delete;
  [[nodiscard]] MathSize grass_curing(const int nd, const wx::FwiWeather& wx) const override;
  [[nodiscard]] MathSize cbh() const override;
  [[nodiscard]] MathSize cfl() const override;
  [[nodiscard]] MathSize buiEffect(MathSize) const override;
  [[nodiscard]] MathSize crownConsumption(MathSize) const override;
  [[nodiscard]] MathSize calculateRos(int, const wx::FwiWeather&, MathSize) const
    override;
  [[nodiscard]] MathSize calculateIsf(const SpreadInfo&, MathSize) const
    override;
  [[nodiscard]] MathSize surfaceFuelConsumption(const SpreadInfo&) const
    override;
  [[nodiscard]] MathSize lengthToBreadth(MathSize) const override;
  [[nodiscard]] MathSize finalRos(const SpreadInfo&,
                                  MathSize,
                                  MathSize,
                                  MathSize) const override;
  [[nodiscard]] MathSize criticalSurfaceIntensity(const SpreadInfo&) const
    override;
  [[nodiscard]] MathSize crownFractionBurned(MathSize,
                                             MathSize) const noexcept override;
  [[nodiscard]] ThresholdSize probabilityPeat(MathSize) const noexcept override;
  [[nodiscard]] ThresholdSize survivalProbability(const wx::FwiWeather&) const noexcept
    override;
};
}
}