Athena++/Atmosphere
Planetary Atmosphere Simulator
thermodynamics.hpp
Go to the documentation of this file.
1 
9 #ifndef THERMODYNAMICS_HPP
10 #define THERMODYNAMICS_HPP
11 
12 // C/C++ headers
13 #include <cfloat>
14 #include <iosfwd>
15 
16 // Athena headers
17 #include "../mesh/mesh.hpp"
18 #include "../eos/eos.hpp"
19 #include "../athena.hpp"
20 #include "thermodynamic_funcs.hpp"
22 #include "vapors/water_vapors.hpp"
23 
24 class MeshBlock;
25 class ParameterInput;
26 
27 // dynamic variables are ordered in the array as the following
28 // 0: dry air (non-condensible)
29 // 1+iphase*NVAPOR:1+(iphase+1)*NVAPOR
30 // iphase = 0 - moist air (condensible)
31 // iphase = 1 - primary condensible species
32 // iphase = 2..(N-1) - all other condensible species
33 
34 enum class Adiabat {reversible = 0, pseudo = 1, dry = 2, isothermal = 3};
35 
37  friend std::ostream& operator<<(std::ostream& os, Thermodynamics const& my);
38 public:
39  // members
42  // member functions
45 
49  Real GetRd() const {
50  return Rd_;
51  }
52 
57  Real GetCvRatio(int i) const {
58  return cv_ratios_[i];
59  }
60 
68  Real GetCv(int i) const {
69  Real cvd = Rd_/(pmy_block->peos->GetGamma() - 1.);
70  return cv_ratios_[i]*cvd;
71  }
72 
77  Real GetCpRatio(int i) const {
78  return cp_ratios_[i];
79  }
80 
88  Real GetCp(int i) const {
89  Real gamma = pmy_block->peos->GetGamma();
90  Real cpd = Rd_*gamma /(gamma - 1.);
91  return cp_ratios_[i]*cpd;
92  }
93 
102  Real GetLatent(int j, Real temp = 0.) const {
103  return latent_[j] - delta_[j]*Rd_/mu_ratios_[j]*temp;
104  }
105 
110  Real GetMassRatio(int i) const {
111  return mu_ratios_[i];
112  }
113 
124  Real GetBeta(int j) const {
125  return beta_[j];
126  }
127 
137  Real GetDelta(int j) const {
138  return delta_[j];
139  }
140 
146  void ConstructAtmosphere(Real **w, Real Ts, Real Ps,
147  Real grav, Real dz, int len, Adiabat method, Real dTdz = 0.) const;
148 
149  // uhat is the molar interal energy defined as:
150  // u^\hat = (q_d^\hat c_d^\hat T
151  // + \sum_i q_i^\hat c_i^\hat T
152  // + \sum_{i,j} q_{ij}^\hat c_{ij}^\hat T
153  // - \sum_{i,j} q_{ij}^\hat \mu_{ij}^\hat)/R_d^\hat
154  void UpdateTPConservingU(Real q[], Real rho, Real uhat) const;
155 
159  //void SaturationAdjustment(AthenaArray<Real> &u, AthenaArray<Real> &c) const;
160 
172  int il, int iu, int jl, int ju, int kl, int ku) const;
173 
174  // Thermodynamic variables to conserved variables
175  // density of dry air, momentum and total energy are not updated
177  int il, int iu, int jl, int ju, int kl, int ku) const;
178 
179  //void PolytropicIndex(AthenaArray<Real> &gamma, AthenaArray<Real> &w,
180  // int kl, int ku, int jl, int ju, int il, int iu) const;
181 
183  template<typename T>
184  Real GetGamma(T w) {
185  Real gamma = pmy_block->peos->GetGamma();
186  Real fsig = 1., feps = 1.;
187  for (int n = 1; n <= NVAPOR; ++n) {
188  fsig += w[n]*(cv_ratios_[n] - 1.);
189  feps += w[n]*(1./mu_ratios_[n] - 1.);
190  }
191  return 1. + (gamma - 1.)*feps/fsig;
192  }
193 
194  // template functions
195  template<typename T>
196  void ChemicalToPrimitive(T w, Real const q[]) {
198  for (int n = 0; n < NHYDRO; ++n) w[n] = w1_[n];
199  }
200 
201  template<typename T>
202  void PrimitiveToChemical(Real q[], T const w) {
203  for (int n = 0; n < NHYDRO; ++n) w1_[n] = w[n];
205  }
206 
207  template<typename T>
208  Real RovRd(T w) {
209  Real feps = 1.;
210  for (int n = 1; n <= NVAPOR; ++n)
211  feps += w[n]*(1./mu_ratios_[n] - 1.);
212  return feps;
213  }
214 
216  template<typename T>
217  Real GetTemp(T w) {
218  return w[IPR]/(w[IDN]*Rd_*RovRd(w));
219  }
220 
221  template<typename T>
222  Real GetPres(T u) {
223  Real gm1 = pmy_block->peos->GetGamma() - 1;
224  Real fsig = 0., feps = 0., rho = u[IDN];
225  for (int n = 0; n <= NVAPOR; ++n) {
226  rho += u[n];
227  fsig += u[n]*cv_ratios_[n];
228  feps += u[n]/mu_ratios_[n];
229  }
230  Real KE = 0.5*(u[IM1]*u[IM1] + u[IM2]*u[IM2] + u[IM3]*u[IM3])/rho;
231  return gm1*(u[IEN] - KE)*feps/fsig;
232  }
233 
234  template<typename T>
235  Real GetChi(T w) {
236  for (int n = 0; n < NHYDRO; ++n) w1_[n] = w[n];
237  Real gamma = pmy_block->peos->GetGamma();
238  Real tem[1] = {GetTemp(w)};
239  update_gamma(gamma, tem);
240  Real qsig = 1., feps = 1.;
241  for (int n = 1; n <= NVAPOR; ++n) {
242  qsig += w[n]*(cp_ratios_[n] - 1.);
243  feps += w[n]*(1./mu_ratios_[n] - 1.);
244  }
245  return (gamma - 1.)/gamma*feps/qsig;
246  }
247 
248  template<typename T>
249  Real GetCp(T w) {
250  Real gamma = pmy_block->peos->GetGamma();
251  Real tem[1] = {GetTemp(w)};
252  update_gamma(gamma, tem);
253  Real qsig = 1.;
254  for (int n = 1; n <= NVAPOR; ++n)
255  qsig += w[n]*(cp_ratios_[n] - 1.);
256  return gamma/(gamma - 1.)*Rd_*qsig;
257  }
258 
259  template<typename T>
260  Real GetCv(T w) {
261  Real gamma = pmy_block->peos->GetGamma();
262  Real tem[1] = {GetTemp(w)};
263  update_gamma(gamma, tem);
264  Real qsig = 1.;
265  for (int n = 1; n <= NVAPOR; ++n)
266  qsig += w[n]*(cv_ratios_[n] - 1.);
267  return 1./(gamma - 1.)*Rd_*qsig;
268  }
269 
271  template<typename T>
272  Real GetTheta(T w, Real p0) {
273  Real chi = GetChi(w);
274  Real temp = GetTemp(w);
275  return temp*pow(p0/w[IPR], chi);
276  }
277 
282  template<typename T>
284  // mass to molar mixing ratio
285  if (vtype == VariableType::prim) {
286  Real sum = 1.;
287  for (int n = 1; n <= NVAPOR; ++n) {
288  w1_[n] = v[n]/mu_ratios_[n];
289  sum += v[n]*(1./mu_ratios_[n] - 1.);
290  }
291  for (int n = 1; n <= NVAPOR; ++n)
292  w1_[n] /= sum;
293 
294  w1_[IDN] = GetTemp(v);
295  w1_[IPR] = v[IPR];
296  } else if (vtype == VariableType::cons) {
297  Real sum = 0., feps = 0.;
298  for (int n = 0; n <= NVAPOR; ++n) {
299  w1_[n] = v[n]/mu_ratios_[n];
300  sum += w1_[n];
301  }
302  for (int n = 0; n <= NVAPOR; ++n)
303  w1_[n] /= sum;
304  w1_[IPR] = GetPres(v);
305  w1_[IDN] = w1_[IPR]/(sum*Rd_);
306  } else { // vtype == VariableType::chem
307  for (int n = 0; n < NHYDRO+2*NVAPOR; ++n)
308  w1_[n] = v[n];
309  }
310 
311  for (int iv = 1; iv <= NVAPOR; ++iv) {
312  int nc = w1_[IDN] > t3_[iv] ? iv + NVAPOR : iv + 2*NVAPOR;
313  int ic = NHYDRO - NVAPOR + nc - 1;
314  Real rate = VaporCloudEquilibrium(w1_, iv, ic, t3_[iv], p3_[iv],
315  0., beta_[nc], delta_[nc], true);
316  dw[iv] = rate/w1_[iv]*v[iv];
317  }
318  }
319 
321  template<typename T>
322  Real GetRelativeHumidity(T w, int iv) {
324  return w[iv]/(w[iv] - dw_[iv]);
325  }
326 
328  //template<typename T>
329  //Real GetThetaE(T prim, Real p0);
330 
331 private:
334 
336  Real w1_[NHYDRO+2*NVAPOR];
337 
339  Real dw_[1+NVAPOR];
340 
341  // read from inputs
345  Real mu_ratios_[1+3*NVAPOR];
347  Real cp_ratios_[1+3*NVAPOR];
351  Real beta_[1+3*NVAPOR];
353  Real t3_[1+NVAPOR];
355  Real p3_[1+NVAPOR];
356 
357  // calculated quantities
361  Real latent_[1+3*NVAPOR];
362 
366  Real delta_[1+3*NVAPOR];
367 
369  Real cv_ratios_[1+3*NVAPOR];
370 };
371 
372 //#if (NVAPOR > 0)
373 // #include "get_theta_e.hpp"
374 //#endif
375 
376 #endif
@ IPR
Definition: athena.hpp:139
double Real
Definition: athena.hpp:29
@ IM3
Definition: athena.hpp:135
@ IEN
Definition: athena.hpp:135
@ IM1
Definition: athena.hpp:135
@ IDN
Definition: athena.hpp:135
@ IM2
Definition: athena.hpp:135
VariableType
Definition: athena.hpp:168
Real GetGamma() const
Definition: eos.hpp:171
EquationOfState * peos
Definition: mesh.hpp:130
Real GetDelta(int j) const
MeshBlock * pmy_block
Real dw_[1+NVAPOR]
scratch array for storing variables
Real GetTemp(T w)
Temperature.
Real GetCpRatio(int i) const
void UpdateTPConservingU(Real q[], Real rho, Real uhat) const
void ThermoToConserved(AthenaArray< Real > &u, AthenaArray< Real > const &q, int il, int iu, int jl, int ju, int kl, int ku) const
Real GetBeta(int j) const
Real mu_ratios_[1+3 *NVAPOR]
ratio of mean molecular weights
void ConservedToThermo(AthenaArray< Real > &q, AthenaArray< Real > const &u, int il, int iu, int jl, int ju, int kl, int ku) const
void SaturationSurplus(Real dw[], T v, VariableType vtype=VariableType::prim)
Real delta_[1+3 *NVAPOR]
Real cp_ratios_[1+3 *NVAPOR]
ratio of specific heat capacities at constant pressure
Real GetRd() const
Real GetGamma(T w)
polytropic index
Real GetCp(int i) const
Real cv_ratios_[1+3 *NVAPOR]
ratio of specific heat capacities at constant volume
void ConstructAtmosphere(Real **w, Real Ts, Real Ps, Real grav, Real dz, int len, Adiabat method, Real dTdz=0.) const
Real GetTheta(T w, Real p0)
Potential temperature.
void PrimitiveToChemical(Real q[], T const w)
Real GetCv(int i) const
Real GetRelativeHumidity(T w, int iv)
Relative humidity.
Real p3_[1+NVAPOR]
triple point pressure [pa]
Real beta_[1+3 *NVAPOR]
Real t3_[1+NVAPOR]
triple point temperature [K]
void ChemicalToPrimitive(T w, Real const q[])
Real GetMassRatio(int i) const
Real Rd_
ideal gas constant of dry air in J/kg
Thermodynamics(MeshBlock *pmb, ParameterInput *pin)
Real latent_[1+3 *NVAPOR]
Real GetLatent(int j, Real temp=0.) const
Real ftol_
Equivalent potential temperature.
Real w1_[NHYDRO+2 *NVAPOR]
scratch array for storing variables
Real GetCvRatio(int i) const
friend std::ostream & operator<<(std::ostream &os, Thermodynamics const &my)
void update_gamma(Real &gamma, Real rcp[], Real const prim[])
Real temp[400]
Real VaporCloudEquilibrium(Real const q[], int iv, int ic, Real t3, Real p3, Real alpha, Real beta, Real delta, bool no_cloud=false)
void mass_to_molar(Real q[], Real const w[], Real const eps[], Real Rd)
void molar_to_mass(Real w[], Real const q[], Real const eps[], Real Rd)
Adiabat