thermodynamics.hpp

Go to the documentation of this file.

#ifndef SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_
#define SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_
 
// C/C++
#include <array>
#include <cfloat>
#include <iosfwd>
#include <map>
#include <memory>
#include <set>
#include <utility>
#include <vector>
 
// external
#include <yaml-cpp/yaml.h>
 
// athena
#include <athena/athena.hpp>
#include <athena/hydro/hydro.hpp>
#include <athena/mesh/mesh.hpp>
 
// canoe
#include <air_parcel.hpp>
#include <configure.hpp>
#include <constants.hpp>
 
class MeshBlock;
class ParameterInput;
 
using IndexPair = std::pair<int, int>;
using IndexSet = std::vector<int>;
 
using RealArray3 = std::array<Real, 3>;
using RealArrayX = std::vector<Real>;
 
using SatVaporPresFunc1 = Real (*)(AirParcel const &, int i, int j);
using SatVaporPresFunc2 = Real (*)(AirParcel const &, int i, int j, int k);
 
enum { MAX_REACTANT = 3 };
using ReactionIndx = std::array<int, MAX_REACTANT>;
using ReactionStoi = std::array<int, MAX_REACTANT>;
using ReactionInfo = std::pair<ReactionIndx, ReactionStoi>;
 
Real NullSatVaporPres1(AirParcel const &, int, int);
Real NullSatVaporPres2(AirParcel const &, int, int, int);
 
void read_thermo_property(Real var[], char const name[], int len, Real v0,
                          ParameterInput *pin);
Real saha_ionization_electron_density(Real T, Real num, Real ion_ev);
 
inline Real SatVaporPresIdeal(Real t, Real p3, Real beta, Real delta) {
  return p3 * exp((1. - 1. / t) * beta - delta * log(t));
}
 
// Thermodynamic variables are ordered in the array as the following
// 0: dry air (non-condensible gas)
// 1..NVAPOR: moist air (condensible gas)
// NVAPOR+1..NVAPOR+NCLOUD: clouds
 
// 0 is a special buffer place for cloud in equilibrium with vapor at the same
// temperature
 
class Thermodynamics {
 protected:
  enum { NPHASE = NPHASE_LEGACY };
 
  Thermodynamics() {}
  static Thermodynamics *fromLegacyInput(ParameterInput *pin);
  static Thermodynamics *fromYAMLInput(YAML::Node node);
 
 public:
  using SVPFunc1Container = std::vector<std::vector<SatVaporPresFunc1>>;
  using SVPFunc2Container = std::map<IndexPair, SatVaporPresFunc2>;
 
  enum { Size = 1 + NVAPOR + NCLOUD };
 
  enum class Method {
    ReversibleAdiabat = 0,
    PseudoAdiabat = 1,
    DryAdiabat = 2,
    Isothermal = 3,
    NeutralStability = 4
  };
 
  static constexpr Real RefTemp = 300.;
  static constexpr Real RefPres = 1.e5;
 
  // member functions
  ~Thermodynamics();
 
  static Thermodynamics const *GetInstance();
 
  static Thermodynamics const *InitFromAthenaInput(ParameterInput *pin);
 
  static void Destroy();
 
  Real GetRd() const { return Rd_; }
 
  Real GetGammad(AirParcel const &var) const;
 
  Real GetGammadRef() const { return gammad_ref_; }
 
  Real GetMu(int n) const { return mu_[n]; }
 
  Real GetInvMu(int n) const { return inv_mu_[n]; }
 
  Real GetMuRatio(int n) const { return mu_ratio_[n]; }
 
  Real GetInvMuRatio(int n) const { return inv_mu_ratio_[n]; }
 
  Real GetCvRatioMass(int n) const { return cv_ratio_mass_[n]; }
 
  Real GetCvRatioMole(int n) const { return cv_ratio_mole_[n]; }
 
  int GetCloudIndex(int i, int j) const { return cloud_index_set_[i][j]; }
 
  Real GetCvMass(AirParcel const &qfrac, int n) const {
    Real cvd = Rd_ / (GetGammad(qfrac) - 1.);
    return cv_ratio_mass_[n] * cvd;
  }
 
  Real GetCvMassRef(int n) const {
    Real cvd = Rd_ / (gammad_ref_ - 1.);
    return cv_ratio_mass_[n] * cvd;
  }
 
  Real GetCvMole(AirParcel const &qfrac, int n) const {
    return GetCvMass(qfrac, n) * mu_[n];
  }
 
  Real GetCpRatioMass(int n) const { return cp_ratio_mass_[n]; }
 
  Real GetCpRatioMole(int n) const { return cp_ratio_mole_[n]; }
 
  Real GetCpMass(AirParcel const &qfrac, int n) const {
    Real gammad = GetGammad(qfrac);
    Real cpd = Rd_ * gammad / (gammad - 1.);
    return cp_ratio_mass_[n] * cpd;
  }
 
  Real GetCpMassRef(int n) const {
    Real cpd = Rd_ * gammad_ref_ / (gammad_ref_ - 1.);
    return cp_ratio_mass_[n] * cpd;
  }
 
  Real GetCpMole(AirParcel const &qfrac, int n) const {
    return GetCpMass(qfrac, n) * mu_[n];
  }
 
  Real GetChi(AirParcel const &qfrac) const;
 
  Real GetLatentEnergyMass(int n, Real temp = 0.) const {
    return latent_energy_mass_[n] - delta_[n] * Rd_ * inv_mu_ratio_[n] * temp;
  }
 
  Real GetLatentEnergyMole(int n, Real temp = 0.) const {
    return GetLatentEnergyMass(n, temp) * mu_[n];
  }
 
  Real GetLatentHeatMole(int i, std::vector<Real> const &rates,
                         Real temp) const;
 
  Real GetLatentHeatMass(int i, std::vector<Real> const &rates,
                         Real temp) const {
    return GetLatentHeatMole(i, rates, temp) * inv_mu_[i];
  }
 
  RealArrayX TryEquilibriumTP_VaporCloud(AirParcel const &qfrac, int ivapor,
                                         Real cv_hat = 0.,
                                         bool misty = false) const;
 
  RealArray3 TryEquilibriumTP_VaporVaporCloud(AirParcel const &air,
                                              IndexPair ij, Real cv_hat = 0.,
                                              bool misty = false) const;
 
  void Extrapolate(AirParcel *qfrac, Real dzORdlnp, Method method,
                   Real grav = 0., Real userp = 0.) const;
 
  void EquilibrateTP(AirParcel *qfrac) const;
 
  void SaturationAdjustment(AirColumn &ac) const;
 
  Real RovRd(AirParcel const &qfrac) const;
 
  Real PotentialTemp(MeshBlock *pmb, Real p0, int k, int j, int i) const {
    auto &w = pmb->phydro->w;
    return GetTemp(pmb, k, j, i) *
           pow(p0 / w(IPR, k, j, i), GetChi(pmb, k, j, i));
  }
 
  Real EquivalentPotentialTemp(MeshBlock *pmb, Real p0, int v, int k, int j,
                               int i) const;
 
  Real GetTemp(MeshBlock const *pmb, int k, int j, int i) const {
    auto &w = pmb->phydro->w;
    return w(IPR, k, j, i) / (w(IDN, k, j, i) * Rd_ * RovRd(pmb, k, j, i));
  }
 
  Real GetMu(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real RovRd(MeshBlock const *pmb, int k, int j, int i) const {
    Real feps = 1.;
    auto &w = pmb->phydro->w;
#pragma omp simd reduction(+ : feps)
    for (int n = 1; n <= NVAPOR; ++n)
      feps += w(n, k, j, i) * (inv_mu_ratio_[n] - 1.);
    return feps;
  }
 
  Real GetCvMass(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real GetCpMass(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real GetChi(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real GetGamma(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real GetPres(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real GetEnthalpyMass(MeshBlock const *pmb, int k, int j, int i) const;
 
  Real MoistStaticEnergy(MeshBlock const *pmb, Real gz, int k, int j,
                         int i) const;
 
  Real RelativeHumidity(MeshBlock const *pmb, int n, int k, int j, int i) const;
 
  Real RelativeHumidity(AirParcel const &qfrac, int n) const {
    auto rates = TryEquilibriumTP_VaporCloud(qfrac, n, 0., true);
    return qfrac.w[n] / (qfrac.w[n] + rates[0]);
  }
 
 protected:
  void updateTPConservingU(AirParcel *qfrac, Real rmole, Real umole) const;
 
  Real getInternalEnergyMole(AirParcel const &qfrac) const;
 
  Real getDensityMole(AirParcel const &qfrac) const;
 
  void setTotalEquivalentVapor(AirParcel *qfrac) const;
 
  Real calDlnTDlnP(AirParcel const &qfrac, Real latent[]) const;
 
  void rk4IntegrateLnp(AirParcel *qfrac, Real dlnp, Method method,
                       Real adlnTdlnP) const;
  void rk4IntegrateZ(AirParcel *qfrac, Real dlnp, Method method, Real grav,
                     Real adlnTdlnP) const;
 
  void enrollVaporFunctions();
 
 protected:
  Real Rd_;
 
  Real gammad_ref_;
 
  std::array<Real, Size> mu_ratio_;
 
  std::array<Real, Size> inv_mu_ratio_;
 
  std::array<Real, Size> mu_;
 
  std::array<Real, Size> inv_mu_;
 
  std::array<Real, Size> cp_ratio_mole_;
 
  std::array<Real, Size> cp_ratio_mass_;
 
  std::array<Real, Size> cv_ratio_mole_;
 
  std::array<Real, Size> cv_ratio_mass_;
 
  std::array<Real, Size> latent_energy_mole_;
 
  std::array<Real, Size> latent_energy_mass_;
 
  std::array<Real, Size> beta_;
 
  std::array<Real, Size> delta_;
 
  std::array<Real, 1 + NVAPOR> t3_;
 
  std::array<Real, 1 + NVAPOR> p3_;
 
  SVPFunc1Container svp_func1_;
 
  std::vector<IndexSet> cloud_index_set_;
 
  SVPFunc2Container svp_func2_;
 
  std::map<IndexPair, ReactionInfo> cloud_reaction_map_;
 
  static Thermodynamics *mythermo_;
 
  int sa_max_iter_ = 5;
 
  Real sa_relax_ = 0.8;
 
  Real sa_ftol_ = 0.01;
};
 
#endif  // SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_