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 "moist_adiabat_funcs.hpp"
21 
22 class MeshBlock;
23 class ParameterInput;
24 
25 // dynamic variables are ordered in the array as the following
26 // 0: dry air (non-condensible)
27 // 1+iphase*NVAPOR:1+(iphase+1)*NVAPOR
28 // iphase = 0 - moist air (condensible)
29 // iphase = 1 - primary condensible species
30 // iphase = 2..(N-1) - all other condensible species
31 
32 enum class Adiabat {reversible = 0, pseudo = 1, dry = 2, isothermal = 3};
33 
35  friend std::ostream& operator<<(std::ostream& os, Thermodynamics const& my);
36 public:
37  // members
38  static Real const Rgas;
39  static Real const kBoltz;
40  static Real const kBoltz_cgs;
43  // member functions
46 
50  Real GetRd() const {
51  return Rd_;
52  }
53 
58  Real GetCvRatio(int i) const {
59  return cv_ratios_[i];
60  }
61 
69  Real GetCv(int i) const {
70  Real cvd = Rd_/(pmy_block->peos->GetGamma() - 1.);
71  return cv_ratios_[i]*cvd;
72  }
73 
78  Real GetCpRatio(int i) const {
79  return cp_ratios_[i];
80  }
81 
89  Real GetCp(int i) const {
90  Real gamma = pmy_block->peos->GetGamma();
91  Real cpd = Rd_*gamma /(gamma - 1.);
92  return cp_ratios_[i]*cpd;
93  }
94 
103  Real GetLatent(int j, Real temp = 0.) const {
104  return latent_[j] - delta_[j]*Rd_/mu_ratios_[j]*temp;
105  }
106 
111  Real GetMassRatio(int i) const {
112  return mu_ratios_[i];
113  }
114 
125  Real GetBeta(int j) const {
126  return beta_[j];
127  }
128 
138  Real GetDelta(int j) const {
139  return delta_[j];
140  }
141 
145  void ConstructAtmosphere(Real **w, Real Ts, Real Ps,
146  Real grav, Real dzORdlnp, int len, Adiabat method, Real userp) const;
147 
148  // conversion functions
150  template<typename T1, typename T2>
151  void PrimitiveToChemical(T1 c, T2 const w) const {
152  // set molar mixing ratio
153  Real sum = 1.;
154  for (int n = 1; n <= NVAPOR; ++n) {
155  c[n] = w[n]/mu_ratios_[n];
156  sum += w[n]*(1./mu_ratios_[n] - 1.);
157  }
158  // set pressure, temperature, velocity
159  c[IPR] = w[IPR];
160  c[IDN] = w[IPR]/(w[IDN]*Rd_*sum);
161  c[IVX] = w[IVX];
162  c[IVY] = w[IVY];
163  c[IVZ] = w[IVZ];
164 
165  Real mols = c[IPR]/(c[IDN]*Rgas);
166 #pragma omp simd
167  for (int n = 1; n <= NVAPOR; ++n) {
168  c[n] *= mols/sum;
169  }
170  }
171 
173  template<typename T1, typename T2>
174  void ChemicalToPrimitive(T1 w, T2 const c) const {
175  // set mass mixing ratio
176  Real sum = 1., mols = c[IPR]/(c[IDN]*Rgas);
177  for (int n = 1; n <= NVAPOR; ++n) {
178  w[n] = c[n]/mols*mu_ratios_[n];
179  sum += c[n]/mols*(mu_ratios_[n] - 1.);
180  }
181 #pragma omp simd
182  for (int n = 1; n <= NVAPOR; ++n)
183  w[n] /= sum;
184 
185  // set pressure, density, velocity
186  w[IPR] = c[IPR];
187  w[IDN] = sum*c[IPR]/(c[IDN]*Rd_);
188  w[IVX] = c[IVX];
189  w[IVY] = c[IVY];
190  w[IVZ] = c[IVZ];
191  }
192 
194  template<typename T1, typename T2>
195  void ConservedToChemical(T1 c, T2 const u) const {
196  Real rho = 0., feps = 0., fsig = 0.;
197  for (int n = 0; n <= NVAPOR; ++n) {
198  rho += u[n];
199  c[n] = u[n]/mu_ratios_[n];
200  feps += c[n];
201  fsig += u[n]*cv_ratios_[n];
202  }
203  Real KE = 0.5*(u[IM1]*u[IM1] + u[IM2]*u[IM2] + u[IM3]*u[IM3])/rho;
204  Real gm1 = pmy_block->peos->GetGamma() - 1.;
205  c[IPR] = gm1*(u[IEN] - KE)*feps/fsig;
206  c[IDN] = c[IPR]/(feps*Rd_);
207  c[IVX] = u[IVX]/rho;
208  c[IVY] = u[IVY]/rho;
209  c[IVZ] = u[IVZ]/rho;
210 
211  Real mols = c[IPR]/(Rgas*c[IDN]);
212 #pragma omp simd
213  for (int n = 1; n <= NVAPOR; ++n)
214  c[n] *= mols/feps;
215  }
216 
218  template<typename T1, typename T2>
219  void ChemicalToConserved(T1 u, T2 const c) const {
220  Real sum = 1., mols = c[IPR]/(Rgas*c[IDN]);
221  for (int n = 1; n <= NVAPOR; ++n)
222  sum += c[n]/mols*(mu_ratios_[n] - 1.);
223  Real rho = c[IPR]*sum/(Rd_*c[IDN]);
224  Real cvd = Rd_/(pmy_block->peos->GetGamma() - 1.);
225  u[IDN] = rho;
226  u[IEN] = 0.5*rho*(c[IVX]*c[IVX] + c[IVY]*c[IVY] + c[IVZ]*c[IVZ]);
227  for (int n = 1; n <= NVAPOR; ++n) {
228  u[n] = rho*c[n]/mols*mu_ratios_[n]/sum;
229  u[IDN] -= u[n];
230  u[IEN] += u[n]*cv_ratios_[n]*cvd*c[IDN];
231  }
232  u[IEN] += u[IDN]*cvd*c[IDN];
233  u[IVX] = c[IVX]*rho;
234  u[IVY] = c[IVY]*rho;
235  u[IVZ] = c[IVZ]*rho;
236  }
237 
238  // Get thermodynamic properties
240  template<typename T>
241  Real GetGamma(T w) const {
242  Real gamma = pmy_block->peos->GetGamma();
243  Real fsig = 1., feps = 1.;
244  for (int n = 1; n <= NVAPOR; ++n) {
245  fsig += w[n]*(cv_ratios_[n] - 1.);
246  feps += w[n]*(1./mu_ratios_[n] - 1.);
247  }
248  return 1. + (gamma - 1.)*feps/fsig;
249  }
250 
251  template<typename T>
252  Real RovRd(T w) const {
253  Real feps = 1.;
254  for (int n = 1; n <= NVAPOR; ++n)
255  feps += w[n]*(1./mu_ratios_[n] - 1.);
256  return feps;
257  }
258 
260  template<typename T>
261  Real GetTemp(T w) const {
262  return w[IPR]/(w[IDN]*Rd_*RovRd(w));
263  }
264 
265  template<typename T>
266  Real GetPres(T u) const {
267  Real gm1 = pmy_block->peos->GetGamma() - 1;
268  Real rho = 0., fsig = 0., feps = 0.;
269  for (int n = 0; n <= NVAPOR; ++n) {
270  rho += u[n];
271  fsig += u[n]*cv_ratios_[n];
272  feps += u[n]/mu_ratios_[n];
273  }
274  Real KE = 0.5*(u[IM1]*u[IM1] + u[IM2]*u[IM2] + u[IM3]*u[IM3])/rho;
275  return gm1*(u[IEN] - KE)*feps/fsig;
276  }
277 
278  template<typename T>
279  Real GetChi(T w) const {
280  Real gamma = pmy_block->peos->GetGamma();
281  Real tem[1] = {GetTemp(w)};
282  update_gamma(gamma, tem);
283  Real qsig = 1., feps = 1.;
284  for (int n = 1; n <= NVAPOR; ++n) {
285  qsig += w[n]*(cp_ratios_[n] - 1.);
286  feps += w[n]*(1./mu_ratios_[n] - 1.);
287  }
288  return (gamma - 1.)/gamma*feps/qsig;
289  }
290 
291  template<typename T>
292  Real GetMeanCp(T w) const {
293  Real gamma = pmy_block->peos->GetGamma();
294  Real tem[1] = {GetTemp(w)};
295  update_gamma(gamma, tem);
296  Real qsig = 1.;
297  for (int n = 1; n <= NVAPOR; ++n)
298  qsig += w[n]*(cp_ratios_[n] - 1.);
299  return gamma/(gamma - 1.)*Rd_*qsig;
300  }
301 
302  template<typename T>
303  Real GetMeanCv(T w) const {
304  Real gamma = pmy_block->peos->GetGamma();
305  Real tem[1] = {GetTemp(w)};
306  update_gamma(gamma, tem);
307  Real qsig = 1.;
308  for (int n = 1; n <= NVAPOR; ++n)
309  qsig += w[n]*(cv_ratios_[n] - 1.);
310  return 1./(gamma - 1.)*Rd_*qsig;
311  }
312 
313  template<typename T>
315  Real feps = 1.;
316  for (int n = 1; n <= NVAPOR; ++n)
317  feps += w[n]*(mu_ratios_[n] - 1.);
318  return Rgas/Rd_*feps;
319  }
320 
325  template<typename T>
326  void SaturationSurplus(Real dv[], T v, VariableType vtype) const {
327  Real q1[NHYDRO];
328  // mass to molar mixing ratio
329  if (vtype == VariableType::prim) {
330  PrimitiveToChemical(q1, v);
331  } else if (vtype == VariableType::cons) {
332  ConservedToChemical(q1, v);
333  } else { // VariableType::chem
334  for (int n = 0; n < NHYDRO; ++n)
335  q1[n] = v[n];
336  }
337  // change molar density to molar mixing ratio
338  Real mols = q1[IPR]/(q1[IDN]*Rgas);
339  for (int n = 1; n <= NVAPOR; ++n) q1[n] /= mols;
340 
341  for (int iv = 1; iv <= NVAPOR; ++iv) {
342  int nc = q1[IDN] > t3_[iv] ? iv + NVAPOR : iv + 2*NVAPOR;
343  int ic = NHYDRO - NVAPOR + nc - 1;
344  Real rate = VaporCloudEquilibrium(q1, iv, ic, t3_[iv], p3_[iv],
345  0., beta_[nc], delta_[nc], true);
346  dv[iv] = rate/q1[iv]*v[iv];
347  }
348  }
349 
350 private:
353 
355  // Real w1_[NHYDRO+2*NVAPOR];
356 
358  // Real dw_[1+NVAPOR];
359 
360  // read from inputs
364  Real mu_ratios_[1+3*NVAPOR];
366  Real cp_ratios_[1+3*NVAPOR];
370  Real beta_[1+3*NVAPOR];
372  Real t3_[1+NVAPOR];
374  Real p3_[1+NVAPOR];
375 
376  // calculated quantities
380  Real latent_[1+3*NVAPOR];
381 
385  Real delta_[1+3*NVAPOR];
386 
388  Real cv_ratios_[1+3*NVAPOR];
389 };
390 
391 #endif
@ IPR
Definition: athena.hpp:139
@ IVZ
Definition: athena.hpp:139
@ IVY
Definition: athena.hpp:139
@ IVX
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:133
Real GetDelta(int j) const
void ChemicalToPrimitive(T1 w, T2 const c) const
Change molar mixing ratio to mass mixing ratio.
MeshBlock * pmy_block
Real GetPres(T u) const
Real GetCpRatio(int i) const
Real GetBeta(int j) const
static Real const kBoltz_cgs
Real GetChi(T w) const
static Real const kBoltz
Real mu_ratios_[1+3 *NVAPOR]
ratio of mean molecular weights
static Real const Rgas
Real GetMeanCv(T w) const
void ChemicalToConserved(T1 u, T2 const c) const
Change molar mixing ratio to density.
Real GetMeanMolecularWeight(T w) const
Real delta_[1+3 *NVAPOR]
Real cp_ratios_[1+3 *NVAPOR]
ratio of specific heat capacities at constant pressure
void PrimitiveToChemical(T1 c, T2 const w) const
Change mass mixing ratio to molar mixing ratio.
Real GetRd() const
Real RovRd(T w) const
Real GetCp(int i) const
Real GetMeanCp(T w) const
Real GetTemp(T w) const
Temperature.
Real GetGamma(T w) const
polytropic index
void SaturationSurplus(Real dv[], T v, VariableType vtype) const
Real cv_ratios_[1+3 *NVAPOR]
ratio of specific heat capacities at constant volume
Real GetCv(int i) const
Real p3_[1+NVAPOR]
triple point pressure [pa]
Real beta_[1+3 *NVAPOR]
Real t3_[1+NVAPOR]
triple point temperature [K]
void ConstructAtmosphere(Real **w, Real Ts, Real Ps, Real grav, Real dzORdlnp, int len, Adiabat method, Real userp) const
Real GetMassRatio(int i) const
Real Rd_
scratch array for storing variables
Thermodynamics(MeshBlock *pmb, ParameterInput *pin)
Real latent_[1+3 *NVAPOR]
Real GetLatent(int j, Real temp=0.) const
Real GetCvRatio(int i) const
friend std::ostream & operator<<(std::ostream &os, Thermodynamics const &my)
void ConservedToChemical(T1 c, T2 const u) const
Change density to molar mixing ratio.
void update_gamma(Real &gamma, Real const q[])
Real VaporCloudEquilibrium(Real const q[], int iv, int ic, Real t3, Real p3, Real alpha, Real beta, Real delta, bool no_cloud=false)
float temp
Definition: fit_ammonia.py:6
Adiabat