Athena++/Atmosphere
Planetary Atmosphere Simulator
adiabatic_hydro.cpp
Go to the documentation of this file.
1 //========================================================================================
2 // Athena++ astrophysical MHD code
3 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
4 // Licensed under the 3-clause BSD License, see LICENSE file for details
5 //========================================================================================
7 // \brief implements functions in class EquationOfState for adiabatic hydrodynamics`
8 
9 // C headers
10 
11 // C++ headers
12 #include <cmath> // sqrt()
13 #include <sstream>
14 #include <stdexcept>
15 
16 // Athena++ headers
17 #include "../athena.hpp"
18 #include "../athena_arrays.hpp"
19 #include "../field/field.hpp"
20 #include "../hydro/hydro.hpp"
21 #include "../mesh/mesh.hpp"
22 #include "../parameter_input.hpp"
23 #include "../thermodynamics/thermodynamics.hpp"
24 #include "../globals.hpp"
25 #include "../debugger/debugger.hpp"
26 #include "eos.hpp"
27 
28 // EquationOfState constructor
29 
31  pmy_block_(pmb),
32  gamma_{pin->GetReal("hydro", "gamma")},
33  density_floor_{pin->GetOrAddReal("hydro", "dfloor", std::sqrt(1024*float_min))},
34  pressure_floor_{pin->GetOrAddReal("hydro", "pfloor", std::sqrt(1024*float_min))},
35  scalar_floor_{pin->GetOrAddReal("hydro", "sfloor", std::sqrt(1024*float_min))} {}
36 
37 //----------------------------------------------------------------------------------------
38 // \!fn void EquationOfState::ConservedToPrimitive(AthenaArray<Real> &cons,
39 // const AthenaArray<Real> &prim_old, const FaceField &b,
40 // AthenaArray<Real> &prim, AthenaArray<Real> &bcc, Coordinates *pco,
41 // int il, int iu, int jl, int ju, int kl, int ku)
42 // \brief Converts conserved into primitive variables in adiabatic hydro.
43 
45  AthenaArray<Real> &cons, const AthenaArray<Real> &prim_old, const FaceField &b,
47  Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku) {
48  Real gm1 = GetGamma() - 1.0;
49  std::stringstream msg;
50  Thermodynamics *pthermo = pmy_block_->pthermo;
51 
52  for (int k=kl; k<=ku; ++k) {
53  for (int j=jl; j<=ju; ++j) {
54 #pragma omp simd
55  for (int i=il; i<=iu; ++i) {
56  Real& u_d = cons(IDN,k,j,i);
57  Real& u_m1 = cons(IM1,k,j,i);
58  Real& u_m2 = cons(IM2,k,j,i);
59  Real& u_m3 = cons(IM3,k,j,i);
60  Real& u_e = cons(IEN,k,j,i);
61 
62  Real& w_d = prim(IDN,k,j,i);
63  Real& w_vx = prim(IVX,k,j,i);
64  Real& w_vy = prim(IVY,k,j,i);
65  Real& w_vz = prim(IVZ,k,j,i);
66  Real& w_p = prim(IPR,k,j,i);
67 
68  // apply density floor, without changing momentum or energy
69  u_d = (u_d > density_floor_) ? u_d : density_floor_;
70 
71  Real density = 0.;
72  for (int n = 0; n <= NVAPOR; ++n)
73  density += cons(n,k,j,i);
74  w_d = density;
75  Real di = 1./density;
76 
77  // mass mixing ratio
78  for (int n = 1; n <= NVAPOR; ++n)
79  prim(n,k,j,i) = cons(n,k,j,i)*di;
80 
81  #ifdef DEBUG
82  if (std::isnan(w_d) || (w_d < density_floor_)) { // IDN may be NAN
83  msg << "### FATAL ERROR in function ConservedToPrimitive"
84  << std::endl << "Density reaches lowest value: " << w_d
85  << std::endl << "At position ("
86  << k << "," << j << "," << i << ") in rank " << Globals::my_rank <<
87  std::endl;
88  for (int ii = std::max(i-3, il); ii <= std::min(i+3, iu); ++ii) {
89  msg << "i = " << ii << " ";
90  for (int jj = std::max(j-3, jl); jj <= std::min(j+3, ju); ++jj)
91  msg << cons(IDN,k,jj,ii) << " ";
92  msg << std::endl;
93  }
94  ATHENA_ERROR(msg);
95  }
96  #endif
97 
98  //Real di = 1.0/u_d;
99  w_vx = u_m1*di;
100  w_vy = u_m2*di;
101  w_vz = u_m3*di;
102 
103  // internal energy
104  Real KE = 0.5*di*(u_m1*u_m1 + u_m2*u_m2 + u_m3*u_m3);
105  Real fsig = 1., feps = 1.;
106  // vapors
107  for (int n = 1; n <= NVAPOR; ++n) {
108  fsig += prim(n,k,j,i)*(pthermo->GetCvRatio(n) - 1.);
109  feps += prim(n,k,j,i)*(1./pthermo->GetMassRatio(n) - 1.);
110  }
111  w_p = gm1*(u_e - KE)*feps/fsig;
112 
113  // apply pressure floor, correct total energy
114  u_e = (w_p > pressure_floor_) ? u_e : ((pressure_floor_/gm1)*fsig/feps + KE);
115  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
116 
117  #ifdef DEBUG
118  if (std::isnan(w_p) || (w_p < pressure_floor_)) {
119  msg << "### FATAL ERROR in function ConservedToPrimitive"
120  << std::endl << "Pressure reaches lowest value: " << w_p
121  << std::endl << "At position ("
122  << k << "," << j << "," << i << ") in rank " << Globals::my_rank <<
123  std::endl;
124  ATHENA_ERROR(msg);
125  }
126  #endif
127  }
128  }
129  }
130 
131 #if DEBUG_LEVEL > 0
132  Debugger *pdbg = pmy_block_->pdebug;
133  pdbg = pdbg->StartTracking("EquationOfStates::ConservedToPrimitive");
134  pdbg->Track3D("rho", IsPositive, prim, IDN);
135  pdbg->Track3D("pres", IsPositive, prim, IPR);
136 #endif
137 
138  return;
139 }
140 
141 //----------------------------------------------------------------------------------------
142 // \!fn void EquationOfState::PrimitiveToConserved(const AthenaArray<Real> &prim,
143 // const AthenaArray<Real> &bc, AthenaArray<Real> &cons, Coordinates *pco,
144 // int il, int iu, int jl, int ju, int kl, int ku);
145 // \brief Converts primitive variables into conservative variables
146 
148  const AthenaArray<Real> &prim, const AthenaArray<Real> &bc,
149  AthenaArray<Real> &cons, Coordinates *pco,
150  int il, int iu, int jl, int ju, int kl, int ku) {
151  Real igm1 = 1.0/(GetGamma() - 1.0);
152  Thermodynamics *pthermo = pmy_block_->pthermo;
153 
154  // Force outer-loop vectorization
155 #pragma omp simd
156  for (int k=kl; k<=ku; ++k) {
157  for (int j=jl; j<=ju; ++j) {
158  //#pragma omp simd
159 #pragma novector
160  for (int i=il; i<=iu; ++i) {
161  Real& u_d = cons(IDN,k,j,i);
162  Real& u_m1 = cons(IM1,k,j,i);
163  Real& u_m2 = cons(IM2,k,j,i);
164  Real& u_m3 = cons(IM3,k,j,i);
165  Real& u_e = cons(IEN,k,j,i);
166 
167  const Real& w_d = prim(IDN,k,j,i);
168  const Real& w_vx = prim(IVX,k,j,i);
169  const Real& w_vy = prim(IVY,k,j,i);
170  const Real& w_vz = prim(IVZ,k,j,i);
171  const Real& w_p = prim(IPR,k,j,i);
172 
173  // density
174  u_d = w_d;
175  for (int n = 1; n <= NVAPOR; ++n) {
176  cons(n,k,j,i) = prim(n,k,j,i)*w_d;
177  cons(IDN,k,j,i) -= cons(n,k,j,i);
178  }
179 
180  // momentum
181  u_m1 = w_vx*w_d;
182  u_m2 = w_vy*w_d;
183  u_m3 = w_vz*w_d;
184 
185  // total energy
186  Real KE = 0.5*w_d*(w_vx*w_vx + w_vy*w_vy + w_vz*w_vz);
187  Real fsig = 1., feps = 1.;
188  // vapors
189  for (int n = 1; n <= NVAPOR; ++n) {
190  fsig += prim(n,k,j,i)*(pthermo->GetCvRatio(n) - 1.);
191  feps += prim(n,k,j,i)*(1./pthermo->GetMassRatio(n) - 1.);
192  }
193  u_e = igm1*w_p*fsig/feps + KE ;
194  }
195  }
196  }
197 
198  return;
199 }
200 
201 //----------------------------------------------------------------------------------------
202 // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
203 // \brief returns adiabatic sound speed given vector of primitive variables
204 Real EquationOfState::SoundSpeed(const Real prim[NHYDRO]) {
205  Thermodynamics *pthermo = pmy_block_->pthermo;
206 
207  Real fsig = 1., feps = 1.;
208  for (int n = 1; n <= NVAPOR; ++n) {
209  fsig += prim[n]*(pthermo->GetCvRatio(n) - 1.);
210  feps += prim[n]*(1./pthermo->GetMassRatio(n) - 1.);
211  }
212 
213  return std::sqrt((1. + (gamma_ - 1)*feps/fsig)*prim[IPR]/prim[IDN]);
214 }
215 
216 //----------------------------------------------------------------------------------------
217 // \!fn void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real> &prim, int k, int j,
218 // =int i)
219 // \brief Apply density and pressure floors to reconstructed L/R cell interface states
221  Real& w_d = prim(IDN,i);
222  Real& w_p = prim(IPR,i);
223 
224  // apply (prim) density floor
225  w_d = (w_d > density_floor_) ? w_d : density_floor_;
226  // apply pressure floor
227  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
228 
229  return;
230 }
231 
232 //----------------------------------------------------------------------------------------
233 // \!fn void EquationOfState::ApplyPrimitiveConservedFloors(AthenaArray<Real> &prim,
234 // AthenaArray<Real> &cons, FaceField &b, int k, int j, int i) {
235 // \brief Apply pressure (prim) floor and correct energy (cons) (typically after W(U))
238  int k, int j, int i) {
239  Real gm1 = GetGamma() - 1.0;
240  Real& w_d = prim(IDN,k,j,i);
241  Real& w_p = prim(IPR,k,j,i);
242 
243  Real& u_d = cons(IDN,k,j,i);
244  Real& u_e = cons(IEN,k,j,i);
245  // apply (prim) density floor, without changing momentum or energy
246  w_d = (w_d > density_floor_) ? w_d : density_floor_;
247  // ensure cons density matches
248  u_d = w_d;
249 
250  Real e_k = 0.5*w_d*(SQR(prim(IVX,k,j,i)) + SQR(prim(IVY,k,j,i))
251  + SQR(prim(IVZ,k,j,i)));
252  // apply pressure floor, correct total energy
253  u_e = (w_p > pressure_floor_) ?
254  u_e : ((pressure_floor_/gm1) + e_k);
255  w_p = (w_p > pressure_floor_) ?
256  w_p : pressure_floor_;
257 
258  return;
259 }
@ 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
Real SoundSpeed(const Real prim[(NHYDRO)])
Real gamma_
Definition: eos.hpp:178
void PrimitiveToConserved(const AthenaArray< Real > &prim, const AthenaArray< Real > &bc, AthenaArray< Real > &cons, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
void ApplyPrimitiveFloors(AthenaArray< Real > &prim, int k, int j, int i)
Real pressure_floor_
Definition: eos.hpp:179
EquationOfState(MeshBlock *pmb, ParameterInput *pin)
Real GetGamma() const
Definition: eos.hpp:171
MeshBlock * pmy_block_
Definition: eos.hpp:177
void ApplyPrimitiveConservedFloors(AthenaArray< Real > &prim, AthenaArray< Real > &cons, AthenaArray< Real > &bcc, int k, int j, int i)
void ConservedToPrimitive(AthenaArray< Real > &cons, const AthenaArray< Real > &prim_old, const FaceField &b, AthenaArray< Real > &prim, AthenaArray< Real > &bcc, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
Real density_floor_
Definition: eos.hpp:179
Debugger * pdebug
Definition: mesh.hpp:137
Thermodynamics * pthermo
Definition: mesh.hpp:132
Real GetMassRatio(int i) const
Real GetCvRatio(int i) const
double min(double x1, double x2, double x3)
Definition: core.h:15
double max(double x1, double x2, double x3)
Definition: core.h:18
#define SQR(x)