Canoe
Comprehensive Atmosphere N' Ocean Engine
thermodynamics.hpp
Go to the documentation of this file.
1 #ifndef SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_
2 #define SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_
3 
4 // C/C++
5 #include <array>
6 #include <cfloat>
7 #include <iosfwd>
8 #include <map>
9 #include <memory>
10 #include <set>
11 #include <utility>
12 #include <vector>
13 
14 // external
15 #include <yaml-cpp/yaml.h>
16 
17 // athena
18 #include <athena/athena.hpp>
19 #include <athena/hydro/hydro.hpp>
20 #include <athena/mesh/mesh.hpp>
21 
22 // canoe
23 #include <air_parcel.hpp>
24 #include <configure.hpp>
25 #include <constants.hpp>
26 
27 class MeshBlock;
28 class ParameterInput;
29 
30 using IndexPair = std::pair<int, int>;
31 using IndexSet = std::vector<int>;
32 
33 using RealArray3 = std::array<Real, 3>;
34 using RealArrayX = std::vector<Real>;
35 
36 using SatVaporPresFunc1 = Real (*)(AirParcel const &, int i, int j);
37 using SatVaporPresFunc2 = Real (*)(AirParcel const &, int i, int j, int k);
38 
40 enum { MAX_REACTANT = 3 };
41 using ReactionIndx = std::array<int, MAX_REACTANT>;
42 using ReactionStoi = std::array<int, MAX_REACTANT>;
43 using ReactionInfo = std::pair<ReactionIndx, ReactionStoi>;
44 
45 Real NullSatVaporPres1(AirParcel const &, int, int);
46 Real NullSatVaporPres2(AirParcel const &, int, int, int);
47 
48 void read_thermo_property(Real var[], char const name[], int len, Real v0,
49  ParameterInput *pin);
50 Real saha_ionization_electron_density(Real T, Real num, Real ion_ev);
51 
59 inline Real SatVaporPresIdeal(Real t, Real p3, Real beta, Real delta) {
60  return p3 * exp((1. - 1. / t) * beta - delta * log(t));
61 }
62 
63 // Thermodynamic variables are ordered in the array as the following
64 // 0: dry air (non-condensible gas)
65 // 1..NVAPOR: moist air (condensible gas)
66 // NVAPOR+1..NVAPOR+NCLOUD: clouds
67 
68 // 0 is a special buffer place for cloud in equilibrium with vapor at the same
69 // temperature
70 
72  protected:
73  enum { NPHASE = NPHASE_LEGACY };
74 
78  static Thermodynamics *fromLegacyInput(ParameterInput *pin);
79  static Thermodynamics *fromYAMLInput(YAML::Node node);
80 
81  public:
82  using SVPFunc1Container = std::vector<std::vector<SatVaporPresFunc1>>;
83  using SVPFunc2Container = std::map<IndexPair, SatVaporPresFunc2>;
84 
85  enum { Size = 1 + NVAPOR + NCLOUD };
86 
87  enum class Method {
89  PseudoAdiabat = 1,
90  DryAdiabat = 2,
91  Isothermal = 3,
93  };
94 
95  static constexpr Real RefTemp = 300.;
96  static constexpr Real RefPres = 1.e5;
97 
98  // member functions
100 
102  static Thermodynamics const *GetInstance();
103 
104  static Thermodynamics const *InitFromAthenaInput(ParameterInput *pin);
105 
107  static void Destroy();
108 
111  Real GetRd() const { return Rd_; }
112 
114  Real GetGammad(AirParcel const &var) const;
115 
117  Real GetGammadRef() const { return gammad_ref_; }
118 
122  Real GetMu(int n) const { return mu_[n]; }
123 
127  Real GetInvMu(int n) const { return inv_mu_[n]; }
128 
132  Real GetMuRatio(int n) const { return mu_ratio_[n]; }
133 
137  Real GetInvMuRatio(int n) const { return inv_mu_ratio_[n]; }
138 
142  Real GetCvRatioMass(int n) const { return cv_ratio_mass_[n]; }
143 
147  Real GetCvRatioMole(int n) const { return cv_ratio_mole_[n]; }
148 
152  int GetCloudIndex(int i, int j) const { return cloud_index_set_[i][j]; }
153 
158  Real GetCvMass(AirParcel const &qfrac, int n) const {
159  Real cvd = Rd_ / (GetGammad(qfrac) - 1.);
160  return cv_ratio_mass_[n] * cvd;
161  }
162 
164  Real GetCvMassRef(int n) const {
165  Real cvd = Rd_ / (gammad_ref_ - 1.);
166  return cv_ratio_mass_[n] * cvd;
167  }
168 
171  Real GetCvMole(AirParcel const &qfrac, int n) const {
172  return GetCvMass(qfrac, n) * mu_[n];
173  }
174 
177  Real GetCpRatioMass(int n) const { return cp_ratio_mass_[n]; }
178 
181  Real GetCpRatioMole(int n) const { return cp_ratio_mole_[n]; }
182 
187  Real GetCpMass(AirParcel const &qfrac, int n) const {
188  Real gammad = GetGammad(qfrac);
189  Real cpd = Rd_ * gammad / (gammad - 1.);
190  return cp_ratio_mass_[n] * cpd;
191  }
192 
193  Real GetCpMassRef(int n) const {
194  Real cpd = Rd_ * gammad_ref_ / (gammad_ref_ - 1.);
195  return cp_ratio_mass_[n] * cpd;
196  }
197 
200  Real GetCpMole(AirParcel const &qfrac, int n) const {
201  return GetCpMass(qfrac, n) * mu_[n];
202  }
203 
206  Real GetChi(AirParcel const &qfrac) const;
207 
214  Real GetLatentEnergyMass(int n, Real temp = 0.) const {
215  return latent_energy_mass_[n] - delta_[n] * Rd_ * inv_mu_ratio_[n] * temp;
216  }
217 
222  Real GetLatentEnergyMole(int n, Real temp = 0.) const {
223  return GetLatentEnergyMass(n, temp) * mu_[n];
224  }
225 
226  Real GetLatentHeatMole(int i, std::vector<Real> const &rates,
227  Real temp) const;
228 
229  Real GetLatentHeatMass(int i, std::vector<Real> const &rates,
230  Real temp) const {
231  return GetLatentHeatMole(i, rates, temp) * inv_mu_[i];
232  }
233 
242  RealArrayX TryEquilibriumTP_VaporCloud(AirParcel const &qfrac, int ivapor,
243  Real cv_hat = 0.,
244  bool misty = false) const;
245 
255  IndexPair ij, Real cv_hat = 0.,
256  bool misty = false) const;
257 
264  void Extrapolate(AirParcel *qfrac, Real dzORdlnp, Method method,
265  Real grav = 0., Real userp = 0.) const;
266 
269  void EquilibrateTP(AirParcel *qfrac) const;
270 
274  void SaturationAdjustment(AirColumn &ac) const;
275 
278  Real RovRd(AirParcel const &qfrac) const;
279 
284  Real PotentialTemp(MeshBlock *pmb, Real p0, int k, int j, int i) const {
285  auto &w = pmb->phydro->w;
286  return GetTemp(pmb, k, j, i) *
287  pow(p0 / w(IPR, k, j, i), GetChi(pmb, k, j, i));
288  }
289 
294  Real EquivalentPotentialTemp(MeshBlock *pmb, Real p0, int v, int k, int j,
295  int i) const;
296 
301  Real GetTemp(MeshBlock const *pmb, int k, int j, int i) const {
302  auto &w = pmb->phydro->w;
303  return w(IPR, k, j, i) / (w(IDN, k, j, i) * Rd_ * RovRd(pmb, k, j, i));
304  }
305 
310  Real GetMu(MeshBlock const *pmb, int k, int j, int i) const;
311 
317  Real RovRd(MeshBlock const *pmb, int k, int j, int i) const {
318  Real feps = 1.;
319  auto &w = pmb->phydro->w;
320 #pragma omp simd reduction(+ : feps)
321  for (int n = 1; n <= NVAPOR; ++n)
322  feps += w(n, k, j, i) * (inv_mu_ratio_[n] - 1.);
323  return feps;
324  }
325 
333  Real GetCvMass(MeshBlock const *pmb, int k, int j, int i) const;
334 
342  Real GetCpMass(MeshBlock const *pmb, int k, int j, int i) const;
343 
348  Real GetChi(MeshBlock const *pmb, int k, int j, int i) const;
349 
354  Real GetGamma(MeshBlock const *pmb, int k, int j, int i) const;
355 
358  Real GetPres(MeshBlock const *pmb, int k, int j, int i) const;
359 
365  Real GetEnthalpyMass(MeshBlock const *pmb, int k, int j, int i) const;
366 
371  Real MoistStaticEnergy(MeshBlock const *pmb, Real gz, int k, int j,
372  int i) const;
373 
378  Real RelativeHumidity(MeshBlock const *pmb, int n, int k, int j, int i) const;
379 
380  Real RelativeHumidity(AirParcel const &qfrac, int n) const {
381  auto rates = TryEquilibriumTP_VaporCloud(qfrac, n, 0., true);
382  return qfrac.w[n] / (qfrac.w[n] + rates[0]);
383  }
384 
385  protected:
387  void updateTPConservingU(AirParcel *qfrac, Real rmole, Real umole) const;
388 
389  Real getInternalEnergyMole(AirParcel const &qfrac) const;
390 
391  Real getDensityMole(AirParcel const &qfrac) const;
392 
393  void setTotalEquivalentVapor(AirParcel *qfrac) const;
394 
399  Real calDlnTDlnP(AirParcel const &qfrac, Real latent[]) const;
400 
401  void rk4IntegrateLnp(AirParcel *qfrac, Real dlnp, Method method,
402  Real adlnTdlnP) const;
403  void rk4IntegrateZ(AirParcel *qfrac, Real dlnp, Method method, Real grav,
404  Real adlnTdlnP) const;
405 
407  void enrollVaporFunctions();
408 
409  protected:
411  Real Rd_;
412 
415 
417  std::array<Real, Size> mu_ratio_;
418 
420  std::array<Real, Size> inv_mu_ratio_;
421 
423  std::array<Real, Size> mu_;
424 
426  std::array<Real, Size> inv_mu_;
427 
429  std::array<Real, Size> cp_ratio_mole_;
430 
432  std::array<Real, Size> cp_ratio_mass_;
433 
435  std::array<Real, Size> cv_ratio_mole_;
436 
438  std::array<Real, Size> cv_ratio_mass_;
439 
441  std::array<Real, Size> latent_energy_mole_;
442 
446  std::array<Real, Size> latent_energy_mass_;
447 
455  std::array<Real, Size> beta_;
456 
463  std::array<Real, Size> delta_;
464 
466  std::array<Real, 1 + NVAPOR> t3_;
467 
469  std::array<Real, 1 + NVAPOR> p3_;
470 
473 
475  std::vector<IndexSet> cloud_index_set_;
476 
479 
481  std::map<IndexPair, ReactionInfo> cloud_reaction_map_;
482 
485 
487  int sa_max_iter_ = 5;
488 
490  Real sa_relax_ = 0.8;
491 
493  Real sa_ftol_ = 0.01;
494 };
495 
496 #endif // SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_
std::vector< AirParcel > AirColumn
Definition: air_parcel.hpp:121
Real *const w
Definition: air_parcel.hpp:36
Real RovRd(MeshBlock const *pmb, int k, int j, int i) const
Inverse of the mean molecular weight (no cloud)
Real GetMuRatio(int n) const
Real GetInvMu(int n) const
std::array< Real, Size > cp_ratio_mass_
ratio of specific heat capacities [J/mol] at constant pressure
static constexpr Real RefPres
RealArrayX TryEquilibriumTP_VaporCloud(AirParcel const &qfrac, int ivapor, Real cv_hat=0., bool misty=false) const
Calculate the equilibrium mole transfer by cloud reaction vapor -> cloud.
Real GetCvRatioMass(int n) const
void SaturationAdjustment(AirColumn &ac) const
std::array< Real, Size > cp_ratio_mole_
ratio of specific heat capacities [J/mol] at constant pressure
std::array< Real, Size > cv_ratio_mass_
ratio of specific heat capacities [J/mol] at constant volume
std::array< Real, Size > inv_mu_ratio_
inverse ratio of mean molecular weights
Real GetLatentEnergyMass(int n, Real temp=0.) const
Temperature dependent specific latent energy [J/kg] of condensates at constant volume .
std::array< Real, Size > beta_
Dimensionless latent heat.
Real GetCvRatioMole(int n) const
Real getDensityMole(AirParcel const &qfrac) const
Real GetCpMassRef(int n) const
static constexpr Real RefTemp
void enrollVaporFunctions()
custom functions for vapor (to be overridden by user mods)
Real PotentialTemp(MeshBlock *pmb, Real p0, int k, int j, int i) const
Calculate potential temperature from primitive variable.
Real gammad_ref_
reference polytropic index of dry air
Real MoistStaticEnergy(MeshBlock const *pmb, Real gz, int k, int j, int i) const
Moist static energy.
Real calDlnTDlnP(AirParcel const &qfrac, Real latent[]) const
Calculate moist adiabatic temperature gradient.
std::array< Real, Size > latent_energy_mass_
latent energy at constant volume [J/kg]
std::map< IndexPair, ReactionInfo > cloud_reaction_map_
reaction information map
static Thermodynamics const * GetInstance()
Return a pointer to the one and only instance of Thermodynamics.
int sa_max_iter_
maximum saturation adjustment iterations
Real GetCvMass(AirParcel const &qfrac, int n) const
int GetCloudIndex(int i, int j) const
std::vector< IndexSet > cloud_index_set_
cloud index set
Real GetRd() const
Real sa_ftol_
saturation adjustment temperature tolerance
static Thermodynamics * fromYAMLInput(YAML::Node node)
std::array< Real, Size > inv_mu_
inverse molecular weight [mol/kg]
std::array< Real, 1+NVAPOR > p3_
triple point pressure [pa]
Real GetLatentHeatMole(int i, std::vector< Real > const &rates, Real temp) const
Real RelativeHumidity(MeshBlock const *pmb, int n, int k, int j, int i) const
Relative humidity.
Real GetLatentEnergyMole(int n, Real temp=0.) const
Temperature dependent specific latent energy [J/mol] of condensates at constant volume.
Real GetMu(int n) const
void updateTPConservingU(AirParcel *qfrac, Real rmole, Real umole) const
update T/P
std::vector< std::vector< SatVaporPresFunc1 > > SVPFunc1Container
Real GetCpRatioMass(int n) const
static Thermodynamics * fromLegacyInput(ParameterInput *pin)
std::map< IndexPair, SatVaporPresFunc2 > SVPFunc2Container
SVPFunc1Container svp_func1_
saturation vapor pressure function: Vapor -> Cloud
void Extrapolate(AirParcel *qfrac, Real dzORdlnp, Method method, Real grav=0., Real userp=0.) const
std::array< Real, 1+NVAPOR > t3_
triple point temperature [K]
Real RovRd(AirParcel const &qfrac) const
std::array< Real, Size > mu_
molecular weight [kg/mol]
Real GetChi(AirParcel const &qfrac) const
std::array< Real, Size > delta_
Dimensionless differences in specific heat capacity.
Real GetCvMole(AirParcel const &qfrac, int n) const
Real GetCvMassRef(int n) const
Reference specific heat capacity [J/(kg K)] at constant volume.
static Thermodynamics * mythermo_
pointer to the single Thermodynamics instance
Real GetCpMass(AirParcel const &qfrac, int n) const
Real GetLatentHeatMass(int i, std::vector< Real > const &rates, Real temp) const
Real EquivalentPotentialTemp(MeshBlock *pmb, Real p0, int v, int k, int j, int i) const
Calculate equivalent potential temperature from primitive variable.
Real GetInvMuRatio(int n) const
Real GetCpRatioMole(int n) const
SVPFunc2Container svp_func2_
saturation vapor pressure function: Vapor + Vapor -> Cloud
Real sa_relax_
saturation adjustment relaxation parameter
std::array< Real, Size > latent_energy_mole_
latent energy at constant volume [J/mol]
Real GetTemp(MeshBlock const *pmb, int k, int j, int i) const
Calculate temperature from primitive variable.
std::array< Real, Size > cv_ratio_mole_
ratio of specific heat capacities [J/mol] at constant volume
static void Destroy()
Destroy the one and only instance of Thermodynamics.
std::array< Real, Size > mu_ratio_
ratio of mean molecular weights
Real GetGammadRef() const
reference adiabatic index of dry air [1]
Real Rd_
ideal gas constant of dry air in J/kg
static Thermodynamics const * InitFromAthenaInput(ParameterInput *pin)
void EquilibrateTP(AirParcel *qfrac) const
Real GetCpMole(AirParcel const &qfrac, int n) const
void rk4IntegrateLnp(AirParcel *qfrac, Real dlnp, Method method, Real adlnTdlnP) const
void setTotalEquivalentVapor(AirParcel *qfrac) const
Real GetEnthalpyMass(MeshBlock const *pmb, int k, int j, int i) const
specific enthalpy [J/kg] of the air parcel
Real GetPres(MeshBlock const *pmb, int k, int j, int i) const
RealArray3 TryEquilibriumTP_VaporVaporCloud(AirParcel const &air, IndexPair ij, Real cv_hat=0., bool misty=false) const
Calculate the equilibrium mole transfer by cloud reaction vapor1 + vapor2 -> cloud.
Real GetGammad(AirParcel const &var) const
adiaibatic index of dry air [1]
Real GetGamma(MeshBlock const *pmb, int k, int j, int i) const
Polytropic index.
Real RelativeHumidity(AirParcel const &qfrac, int n) const
Real getInternalEnergyMole(AirParcel const &qfrac) const
void rk4IntegrateZ(AirParcel *qfrac, Real dlnp, Method method, Real grav, Real adlnTdlnP) const
Real gz[18][8]
float temp
Definition: fit_ammonia.py:6
Real p0
Definition: robert.cpp:36
std::vector< int > IndexSet
std::array< Real, 3 > RealArray3
std::pair< int, int > IndexPair
Real NullSatVaporPres2(AirParcel const &, int, int, int)
Real(*)(AirParcel const &, int i, int j) SatVaporPresFunc1
std::pair< ReactionIndx, ReactionStoi > ReactionInfo
std::array< int, MAX_REACTANT > ReactionIndx
Real(*)(AirParcel const &, int i, int j, int k) SatVaporPresFunc2
Real SatVaporPresIdeal(Real t, Real p3, Real beta, Real delta)
std::vector< Real > RealArrayX
Real saha_ionization_electron_density(Real T, Real num, Real ion_ev)
Definition: ionization.cpp:17
Real NullSatVaporPres1(AirParcel const &, int, int)
@ MAX_REACTANT
void read_thermo_property(Real var[], char const name[], int len, Real v0, ParameterInput *pin)
std::array< int, MAX_REACTANT > ReactionStoi