Canoe
Comprehensive Atmosphere N' Ocean Engine
thermodynamics.cpp
Go to the documentation of this file.
1 // C/C++ header
2 #include <cstring>
3 #include <fstream>
4 #include <iomanip>
5 #include <iostream>
6 #include <mutex>
7 #include <sstream>
8 #include <stdexcept>
9 #include <string>
10 #include <vector> // fill
11 
12 // external
13 #include <yaml-cpp/yaml.h>
14 
15 // athena
16 #include <athena/athena.hpp>
17 #include <athena/eos/eos.hpp>
18 #include <athena/hydro/hydro.hpp>
19 #include <athena/mesh/mesh.hpp>
20 #include <athena/parameter_input.hpp>
21 
22 // application
23 #include <application/application.hpp>
24 #include <application/exceptions.hpp>
25 
26 // canoe
27 #include <air_parcel.hpp>
28 #include <configure.hpp>
29 #include <constants.hpp>
30 #include <impl.hpp>
31 #include <index_map.hpp>
32 
33 // climath
34 #include <climath/core.h>
35 
36 // snap
37 #include "molecules.hpp"
38 #include "thermodynamics.hpp"
39 
40 static std::mutex thermo_mutex;
41 
42 Real __attribute__((weak))
43 Thermodynamics::GetGammad(AirParcel const& qfrac) const {
44  return gammad_ref_;
45 }
46 
48  Application::Logger app("snap");
49  app->Log("Destroy Thermodynamics");
50 }
51 
53  mythermo_ = new Thermodynamics();
54 
55  return mythermo_;
56 }
57 
59  if (NCLOUD < (NPHASE - 1) * NVAPOR) {
60  throw RuntimeError(
61  "Thermodynamics",
62  "NCLOUD < (NPHASE-1)*NVAPOR is not supported for legacy input");
63  }
64 
65  mythermo_ = new Thermodynamics();
66 
67  // Read molecular weight ratios
68  read_thermo_property(mythermo_->mu_ratio_.data(), "eps", NPHASE, 1., pin);
69 
70  // Read cp ratios
72  pin);
73 
74  // Read beta parameter
75  read_thermo_property(mythermo_->beta_.data(), "beta", NPHASE, 0., pin);
76 
77  // Read triple point temperature
78  read_thermo_property(mythermo_->t3_.data(), "Ttriple", 1, 0., pin);
79 
80  // Read triple point pressure
81  read_thermo_property(mythermo_->p3_.data(), "Ptriple", 1, 0., pin);
82 
83  mythermo_->cloud_index_set_.resize(1 + NVAPOR);
84  mythermo_->svp_func1_.resize(1 + NVAPOR);
85 
86  for (int i = 0; i <= NVAPOR; ++i) {
87  mythermo_->cloud_index_set_[i].resize(NPHASE - 1);
88  mythermo_->svp_func1_[i].resize(NPHASE - 1);
89  for (int j = 1; j < NPHASE; ++j) {
90  mythermo_->cloud_index_set_[i][j - 1] = (j - 1) * NVAPOR + (i - 1);
91  }
92  }
93 
94  mythermo_->Rd_ = pin->GetOrAddReal("thermodynamics", "Rd", 1.);
95  mythermo_->gammad_ref_ = pin->GetReal("hydro", "gamma");
96 
97  return mythermo_;
98 }
99 
101  // RAII
102  std::unique_lock<std::mutex> lock(thermo_mutex);
103 
104  if (mythermo_ == nullptr) {
105  mythermo_ = new Thermodynamics();
106  }
107 
108  return mythermo_;
109 }
110 
112  if (mythermo_ != nullptr) {
113  throw RuntimeError("Thermodynamics", "Thermodynamics has been initialized");
114  }
115 
116  Application::Logger app("snap");
117  app->Log("Initialize Thermodynamics");
118 
119  if (pin->DoesParameterExist("thermodynamics", "thermo_file")) {
120  std::string filename = pin->GetString("thermodynamics", "thermo_file");
121  std::ifstream stream(filename);
122  if (stream.good() == false) {
123  throw RuntimeError("Thermodynamics",
124  "Cannot open thermodynamic file: " + filename);
125  }
126 
127  mythermo_ = fromYAMLInput(YAML::Load(stream));
128  } else { // legacy input
129  mythermo_ = fromLegacyInput(pin);
130  }
131 
132  // enroll vapor functions
134 
135  // alias
136  auto& Rd = mythermo_->Rd_;
137  auto& gammad = mythermo_->gammad_ref_;
138  auto& mu_ratio = mythermo_->mu_ratio_;
139  auto& inv_mu_ratio = mythermo_->inv_mu_ratio_;
140  auto& mu = mythermo_->mu_;
141  auto& inv_mu = mythermo_->inv_mu_;
142 
143  auto& cp_ratio_mole = mythermo_->cp_ratio_mole_;
144  auto& cp_ratio_mass = mythermo_->cp_ratio_mass_;
145  auto& cv_ratio_mole = mythermo_->cv_ratio_mole_;
146  auto& cv_ratio_mass = mythermo_->cv_ratio_mass_;
147 
148  auto& latent_energy_mass = mythermo_->latent_energy_mass_;
149  auto& latent_energy_mole = mythermo_->latent_energy_mole_;
150 
151  auto& beta = mythermo_->beta_;
152  auto& delta = mythermo_->delta_;
153  auto& t3 = mythermo_->t3_;
154  auto& p3 = mythermo_->p3_;
155  auto& cloud_index_set = mythermo_->cloud_index_set_;
156 
157  // molecular weight
158  mu[0] = Constants::Rgas / Rd;
159 
160  // inverse mu ratio
161  for (int n = 0; n < mu_ratio.size(); ++n) {
162  inv_mu_ratio[n] = 1. / mu_ratio[n];
163  mu[n] = mu[0] * mu_ratio[n];
164  inv_mu[n] = 1. / mu[n];
165  cp_ratio_mole[n] = cp_ratio_mass[n] * mu_ratio[n];
166  }
167 
168  // calculate latent energy = @f$\beta\frac{R_d}{\epsilon}T^r@f$
169  for (int n = 0; n <= NVAPOR; ++n) {
170  latent_energy_mass[n] = 0.;
171  latent_energy_mole[n] = 0.;
172  }
173 
174  for (int i = 1; i <= NVAPOR; ++i) {
175  for (int j = 0; j < NPHASE - 1; ++j) {
176  int n = cloud_index_set[i][j] + 1 + NVAPOR;
177  latent_energy_mass[n] = beta[n] * Rd / mu_ratio[n] * t3[i];
178  latent_energy_mole[n] = latent_energy_mass[n] * mu[n];
179  }
180  }
181 
182  // calculate delta = @f$(\sigma_j - \sigma_i)*\epsilon_i*\gamma/(\gamma - 1)@f$
183  for (int n = 0; n <= NVAPOR; ++n) delta[n] = 0.;
184  for (int i = 1; i <= NVAPOR; ++i) {
185  for (int j = 0; j < cloud_index_set[i].size(); ++j) {
186  int n = cloud_index_set[i][j] + 1 + NVAPOR;
187  delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] /
188  (1. - 1. / gammad);
189  }
190  }
191 
192  // set up extra clouds
193  auto pindex = IndexMap::GetInstance();
194 
195  for (int n = (NPHASE - 1) * NVAPOR; n < NCLOUD; ++n) {
196  Molecule mol;
197  mol.LoadThermodynamicFile(pindex->GetCloudName(n) + ".thermo");
198 
199  int j = 1 + NVAPOR + n;
200  mu[j] = mol.mu();
201  inv_mu[j] = 1. / mu[j];
202 
203  mu_ratio[j] = mol.mu() / mu[0];
204  inv_mu_ratio[j] = 1. / mu_ratio[j];
205 
206  cp_ratio_mass[j] = mol.cp() / (mol.mu() * mythermo_->GetCpMassRef(0));
207  cp_ratio_mole[j] = cp_ratio_mass[j] * mu_ratio[j];
208 
209  latent_energy_mole[j] = mol.latent();
210  latent_energy_mass[j] = mol.latent() / mol.mu();
211 
212  beta[j] = mol.latent() / (Constants::Rgas * Thermodynamics::RefTemp);
213  delta[j] = 0.;
214  }
215 
216  // finally, set up cv ratio
217  // calculate cv_ratio = @f$\sigma_i + (1. - \gamma)/\epsilon_i@f$
218  for (int n = 0; n <= NVAPOR; ++n) {
219  cv_ratio_mass[n] = gammad * cp_ratio_mass[n] + (1. - gammad) / mu_ratio[n];
220  cv_ratio_mole[n] = cv_ratio_mass[n] * mu_ratio[n];
221  }
222 
223  for (int n = 1 + NVAPOR; n < Size; ++n) {
224  cv_ratio_mass[n] = gammad * cp_ratio_mass[n];
225  cv_ratio_mole[n] = cv_ratio_mass[n] * mu_ratio[n];
226  }
227 
228  // saturation adjustment parmaeters
230  pin->GetOrAddReal("thermodynamics", "sa.max_iter", 10);
231  mythermo_->sa_relax_ = pin->GetOrAddReal("thermodynamics", "sa.relax", 0.8);
232  mythermo_->sa_ftol_ = pin->GetOrAddReal("thermodynamics", "sa.ftol", 1.e-4);
233 
234  return mythermo_;
235 }
236 
238  std::unique_lock<std::mutex> lock(thermo_mutex);
239 
240  if (Thermodynamics::mythermo_ != nullptr) {
242  Thermodynamics::mythermo_ = nullptr;
243  }
244 }
245 
246 Real Thermodynamics::GetPres(MeshBlock const* pmb, int k, int j, int i) const {
247  Real gm1 = pmb->peos->GetGamma() - 1;
248  auto& u = pmb->phydro->u;
249 
250  Real rho = 0., fsig = 0., feps = 0.;
251 #pragma omp simd reduction(+ : rho, fsig, feps)
252  for (int n = 0; n <= NVAPOR; ++n) {
253  rho += u(n, k, j, i);
254  fsig += u(n, k, j, i) * cv_ratio_mass_[n];
255  feps += u(n, k, j, i) * inv_mu_ratio_[n];
256  }
257 
258  // TODO(cli): not correct for Cubed sphere
259  Real KE =
260  0.5 *
261  (sqr(u(IM1, k, j, i)) + sqr(u(IM2, k, j, i)) + sqr(u(IM3, k, j, i))) /
262  rho;
263  return gm1 * (u(IEN, k, j, i) - KE) * feps / fsig;
264 }
265 
266 // Eq.71 in Li2019
267 Real Thermodynamics::GetChi(MeshBlock const* pmb, int k, int j, int i) const {
268  Real gammad = pmb->peos->GetGamma();
269  auto& w = pmb->phydro->w;
270 
271  Real qsig = 1., feps = 1.;
272 #pragma omp simd reduction(+ : qsig, feps)
273  for (int n = 1; n <= NVAPOR; ++n) {
274  feps += w(n, k, j, i) * (inv_mu_ratio_[n] - 1.);
275  qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.);
276  }
277 
278  return (gammad - 1.) / gammad * feps / qsig;
279 }
280 
281 // TODO(cli): check
282 Real Thermodynamics::GetChi(AirParcel const& qfrac) const {
283  Real gammad = GetGammad(qfrac);
284 
285  Real qsig = 1., feps = 1.;
286 #pragma omp simd reduction(+ : qsig)
287  for (int n = 1; n <= NVAPOR; ++n) {
288  qsig += qfrac.w[n] * (cp_ratio_mole_[n] - 1.);
289  }
290 
291 #pragma omp simd reduction(+ : qsig, feps)
292  for (int n = 0; n < NCLOUD; ++n) {
293  feps += -qfrac.c[n];
294  qsig += qfrac.c[n] * (cp_ratio_mole_[n + 1 + NVAPOR] - 1.);
295  }
296 
297  return (gammad - 1.) / gammad / qsig;
298 }
299 
300 // Eq.63 in Li2019
301 Real Thermodynamics::GetGamma(MeshBlock const* pmb, int k, int j, int i) const {
302  Real gammad = pmb->peos->GetGamma();
303 
304  auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
305 
306  Real fsig = 1., feps = 1.;
307 #pragma omp simd reduction(+ : fsig, feps)
308  for (int n = 1; n <= NVAPOR; ++n) {
309  fsig += air.w[n] * (cv_ratio_mass_[n] - 1.);
310  feps += air.w[n] * (inv_mu_ratio_[n] - 1.);
311  }
312 
313  for (int n = 0; n < NCLOUD; ++n) {
314  fsig += air.c[n] * (cv_ratio_mass_[n + 1 + NVAPOR] - 1.);
315  feps += -air.c[n];
316  }
317 
318  return 1. + (gammad - 1.) * feps / fsig;
319 }
320 
321 // Eq.94 in Li2019
322 Real Thermodynamics::RovRd(AirParcel const& qfrac) const {
323  Real fgas = 1., feps = 1.;
324 
325 #pragma omp simd reduction(+ : feps)
326  for (int n = 1; n <= NVAPOR; ++n) {
327  feps += qfrac.w[n] * (mu_ratio_[n] - 1.);
328  }
329 
330 #pragma omp simd reduction(+ : fgas, feps)
331  for (int n = 0; n < NCLOUD; ++n) {
332  fgas += -qfrac.c[n];
333  feps += qfrac.c[n] * (mu_ratio_[n] - 1.);
334  }
335 
336  return fgas / feps;
337 }
338 
339 Real Thermodynamics::MoistStaticEnergy(MeshBlock const* pmb, Real gz, int k,
340  int j, int i) const {
341  auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
342  air.ToMassConcentration();
343 
344  Real temp = GetTemp(pmb, k, j, i);
345  Real rho = 0., LE = 0., IE = 0.;
346 
347 #pragma omp simd reduction(+ : IE, rho)
348  for (int n = 0; n <= NVAPOR; ++n) {
349  IE += air.w[n] * GetCpMassRef(n) * temp;
350  rho += air.w[n];
351  }
352 
353 #pragma omp simd reduction(+ : IE, LE, rho)
354  for (int n = 0; n < NCLOUD; ++n) {
355  IE += air.c[n] * GetCpMassRef(1 + NVAPOR + n) * temp;
356  LE += -latent_energy_mass_[1 + NVAPOR + n] * air.c[n];
357  rho += air.c[n];
358  }
359 
360  return (IE + LE) / rho + gz;
361 }
362 
363 // TODO(cli): check
364 Real Thermodynamics::GetCpMass(MeshBlock const* pmb, int k, int j,
365  int i) const {
366  Real gammad = pmb->peos->GetGamma();
367  auto& w = pmb->phydro->w;
368 
369  Real qsig = 1.;
370 #pragma omp simd reduction(+ : qsig)
371  for (int n = 1; n <= NVAPOR; ++n)
372  qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.);
373  return gammad / (gammad - 1.) * Rd_ * qsig;
374 }
375 
376 // TODO(cli): check
377 Real Thermodynamics::GetCvMass(MeshBlock const* pmb, int k, int j,
378  int i) const {
379  Real gammad = pmb->peos->GetGamma();
380  auto& w = pmb->phydro->w;
381 
382  Real qsig = 1.;
383 #pragma omp simd reduction(+ : qsig)
384  for (int n = 1; n <= NVAPOR; ++n)
385  qsig += w(n, k, j, i) * (cv_ratio_mass_[n] - 1.);
386  return 1. / (gammad - 1.) * Rd_ * qsig;
387 }
388 
389 Real Thermodynamics::GetEnthalpyMass(MeshBlock const* pmb, int k, int j,
390  int i) const {
391  Real gammad = pmb->peos->GetGamma();
392  auto& w = pmb->phydro->w;
393 
394  Real qsig = 1.;
395 #pragma omp simd reduction(+ : qsig)
396  for (int n = 1; n <= NVAPOR; ++n)
397  qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.);
398  return gammad / (gammad - 1.) * Rd_ * qsig * w(IDN, k, j, i);
399 }
400 
401 Real Thermodynamics::GetMu(MeshBlock const* pmb, int k, int j, int i) const {
402  auto& w = pmb->phydro->w;
403 
404  Real feps = 1.;
405 #pragma omp simd reduction(+ : feps)
406  for (int n = 1; n <= NVAPOR; ++n) feps += w(n, k, j, i) * (mu_ratio_[n] - 1.);
407  return mu_[0] * feps;
408 }
409 
410 Real Thermodynamics::RelativeHumidity(MeshBlock const* pmb, int n, int k, int j,
411  int i) const {
412  auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
413  air.ToMoleFraction();
414  return RelativeHumidity(air, n);
415 }
416 
417 void Thermodynamics::Extrapolate(AirParcel* qfrac, Real dzORdlnp, Method method,
418  Real grav, Real userp) const {
419  qfrac->ToMoleFraction();
420 
421  // equilibrate vapor with clouds
422  for (int i = 1; i <= NVAPOR; ++i) {
423  auto rates = TryEquilibriumTP_VaporCloud(*qfrac, i);
424 
425  // vapor condensation rate
426  qfrac->w[i] += rates[0];
427 
428  // cloud concentration rates
429  for (int j = 1; j < rates.size(); ++j)
430  qfrac->c[cloud_index_set_[i][j - 1]] += rates[j];
431  }
432 
433  // RK4 integration
434 #ifdef HYDROSTATIC
435  rk4IntegrateLnp(qfrac, dzORdlnp, method, userp);
436 #else
437  rk4IntegrateZ(qfrac, dzORdlnp, method, grav, userp);
438 #endif
439 }
440 
441 Real Thermodynamics::GetLatentHeatMole(int i, std::vector<Real> const& rates,
442  Real temp) const {
443  if (std::abs(rates[0]) < 1.E-8) return 0.;
444 
445  Real heat = 0.;
446  for (int j = 1; j < rates.size(); ++j) {
447  int n = cloud_index_set_[i][j - 1] + 1 + NVAPOR;
448  heat += rates[j] * GetLatentEnergyMole(n, temp);
449  }
450 
451  return heat / std::abs(rates[0]);
452 }
453 
454 // Eq.4.5.11 in Emanuel (1994)
455 Real Thermodynamics::EquivalentPotentialTemp(MeshBlock* pmb, Real p0, int v,
456  int k, int j, int i) const {
457 #if (NVAPOR > 0)
458  AirParcel&& air_mass = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
459 
460  AirParcel air_mole = air_mass;
461  air_mole.ToMoleFraction();
462 
463  // get dry air mixing ratio
464  Real xg = 1., qd = 1.;
465 #pragma omp simd reduction(+ : xg, qd)
466  for (int n = 0; n < NCLOUD; ++n) {
467  xg += -air_mole.c[n];
468  qd += -air_mass.c[n];
469  }
470 
471  Real xd = xg;
472 #pragma omp simd reduction(+ : xd, qd)
473  for (int n = 1; n <= NVAPOR; ++n) {
474  xd += -air_mole.w[n];
475  qd += -air_mass.w[n];
476  }
477 
478  Real temp = air_mole.w[IDN];
479  Real pres = air_mole.w[IPR];
480 
481  Real qv = air_mass.w[v];
482  Real qc = air_mass.c[cloud_index_set_[v][0]];
483  Real qt = qv + qc;
484 
485  Real Rd = Rd_;
486  Real Rv = Rd_ / mu_ratio_[v];
487 
488  Real cpd = GetCpMassRef(0);
489  Real cl = GetCpMassRef(cloud_index_set_[v][0] + 1 + NVAPOR);
490 
491  std::vector<Real> rates{-0.01, 0.01};
492  Real lv = GetLatentHeatMass(v, rates, temp);
493 
494  Real rh = RelativeHumidity(pmb, v, k, j, i);
495  Real pd = xd / xg * pres;
496  Real cpt = cpd * qd + cl * qt;
497 
498  return temp * pow(p0 / pd, Rd * qd / cpt) * pow(rh, -Rv * qv / cpt) *
499  exp(lv * qv / (cpt * temp));
500 #else
501  return PotentialTemp(pmb, p0, k, j, i);
502 #endif
503 }
504 
506  Real umole) const {
507  Real cvd = Constants::Rgas / (GetGammad(*qfrac) - 1.);
508  Real fsig = 1., qgas = 1.;
509 
510  for (int i = 1; i <= NVAPOR; ++i) {
511  // vapor
512  fsig += (cv_ratio_mole_[i] - 1.) * qfrac->w[i];
513 
514  // clouds
515  for (auto j : cloud_index_set_[i]) {
516  int n = j + 1 + NVAPOR;
517  Real qc = qfrac->c[j];
518 
519  fsig += (cv_ratio_mole_[n] - 1.) * qc;
520  umole += latent_energy_mole_[n] * qc;
521  }
522  }
523 
524  // clouds
525 #pragma omp simd reduction(+ : qgas)
526  for (int n = 0; n < NCLOUD; ++n) qgas += -qfrac->c[n];
527 
528  qfrac->w[IDN] = umole / (cvd * fsig);
529  qfrac->w[IPR] = rmole * qgas * Constants::Rgas * qfrac->w[IDN];
530 }
531 
533  Real cvd = Constants::Rgas / (GetGammad(qfrac) - 1.);
534  Real fsig = 1., LE = 0.;
535 
536  for (int i = 1; i <= NVAPOR; ++i) {
537  // vapor
538  fsig += (cv_ratio_mole_[i] - 1.) * qfrac.w[i];
539 
540  // clouds
541  for (auto j : cloud_index_set_[i]) {
542  int n = j + 1 + NVAPOR;
543  Real qc = qfrac.c[j];
544 
545  fsig += (cv_ratio_mole_[n] - 1.) * qc;
546  LE += latent_energy_mole_[n] * qc;
547  }
548  }
549 
550  return cvd * qfrac.w[IDN] * fsig - LE;
551 }
552 
553 Real Thermodynamics::getDensityMole(AirParcel const& qfrac) const {
554  Real qgas = 1.;
555 #pragma omp simd reduction(+ : qgas)
556  for (int n = 0; n < NCLOUD; ++n) qgas += -qfrac.c[n];
557  return qfrac.w[IPR] / (Constants::Rgas * qfrac.w[IDN] * qgas);
558 }
559 
561  // vpaor <=> cloud
562  for (int i = 1; i <= NVAPOR; ++i)
563  for (auto& j : cloud_index_set_[i]) {
564  qfrac->w[i] += qfrac->c[j];
565  qfrac->c[j] = 0.;
566  }
567 
568  // vapor + vapor <=> cloud
569  for (auto const& [ij, info] : cloud_reaction_map_) {
570  auto& indx = info.first;
571  auto& stoi = info.second;
572 
573  qfrac->w[indx[0]] += stoi[0] / stoi[2] * qfrac->c[indx[2]];
574  qfrac->c[indx[2]] = 0.;
575  qfrac->w[indx[1]] += stoi[1] / stoi[2] * qfrac->c[indx[2]];
576  qfrac->c[indx[2]] = 0.;
577  }
578 }
579 
AirParcel & ToMoleFraction()
Definition: air_parcel.cpp:104
Real *const w
Definition: air_parcel.hpp:36
Real *const c
cloud data
Definition: air_parcel.hpp:39
static IndexMap const * GetInstance()
Definition: index_map.cpp:29
virtual double cp(double T) const
Definition: molecules.cpp:57
void LoadThermodynamicFile(std::string chemfile)
Definition: molecules.cpp:36
double mu() const
Definition: molecules.hpp:107
virtual double latent(double T) const
Definition: molecules.cpp:104
std::array< Real, Size > cp_ratio_mass_
ratio of specific heat capacities [J/mol] at constant pressure
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.
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
std::array< Real, Size > beta_
Dimensionless latent heat.
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.
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
std::vector< IndexSet > cloud_index_set_
cloud index set
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
static Thermodynamics * fromLegacyInput(ParameterInput *pin)
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.
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 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 Rd_
ideal gas constant of dry air in J/kg
static Thermodynamics const * InitFromAthenaInput(ParameterInput *pin)
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
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 getInternalEnergyMole(AirParcel const &qfrac) const
void rk4IntegrateZ(AirParcel *qfrac, Real dlnp, Method method, Real grav, Real adlnTdlnP) const
double sqr(double x)
Definition: core.h:14
Real gz[18][8]
AirParcel gather_from_primitive(MeshBlock const *pmb, int k, int j, int i)
Definition: air_parcel.cpp:507
double const Rgas
Definition: constants.hpp:5
float temp
Definition: fit_ammonia.py:6
void read_thermo_property(Real var[], char const name[], int len, Real v0, ParameterInput *pin)
Real p0
Definition: robert.cpp:36
Real Rd
Definition: straka.cpp:66
Real __attribute__((weak)) Thermodynamics
static std::mutex thermo_mutex