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 
53 
54  for (int k=kl; k<=ku; ++k) {
55  for (int j=jl; j<=ju; ++j) {
56 #pragma omp simd
57  for (int i=il; i<=iu; ++i) {
58  Real& u_d = cons(IDN,k,j,i);
59  Real& u_m1 = cons(IM1,k,j,i);
60  Real& u_m2 = cons(IM2,k,j,i);
61  Real& u_m3 = cons(IM3,k,j,i);
62  Real& u_e = cons(IEN,k,j,i);
63 
64  Real& w_d = prim(IDN,k,j,i);
65  Real& w_vx = prim(IVX,k,j,i);
66  Real& w_vy = prim(IVY,k,j,i);
67  Real& w_vz = prim(IVZ,k,j,i);
68  Real& w_p = prim(IPR,k,j,i);
69 
70  // apply density floor, without changing momentum or energy
71  u_d = (u_d > density_floor_) ? u_d : density_floor_;
72 
73  Real density = 0.;
74  for (int n = 0; n <= NVAPOR; ++n)
75  density += cons(n,k,j,i);
76  w_d = density;
77  Real di = 1./density;
78 
79  // mass mixing ratio
80  for (int n = 1; n <= NVAPOR; ++n)
81  prim(n,k,j,i) = cons(n,k,j,i)*di;
82 
83 #if DEBUG_LEVEL > 3
84  if (std::isnan(w_d) || (w_d < density_floor_)) { // IDN may be NAN
85  msg << "### FATAL ERROR in function ConservedToPrimitive"
86  << std::endl << "Density reaches lowest value: " << w_d
87  << std::endl << "At position ("
88  << k << "," << j << "," << i << ") in rank " << Globals::my_rank <<
89  std::endl;
90  for (int ii = std::max(i-3, il); ii <= std::min(i+3, iu); ++ii) {
91  msg << "i = " << ii << " ";
92  for (int jj = std::max(j-3, jl); jj <= std::min(j+3, ju); ++jj)
93  msg << cons(IDN,k,jj,ii) << " ";
94  msg << std::endl;
95  }
96  ATHENA_ERROR(msg);
97  }
98 #endif
99 
100  //Real di = 1.0/u_d;
101  w_vx = u_m1*di;
102  w_vy = u_m2*di;
103  w_vz = u_m3*di;
104 
105  // internal energy
106  Real KE = 0.5*di*(u_m1*u_m1 + u_m2*u_m2 + u_m3*u_m3);
107  Real fsig = 1., feps = 1.;
108  // vapors
109  for (int n = 1; n <= NVAPOR; ++n) {
110  fsig += prim(n,k,j,i)*(pthermo->GetCvRatio(n) - 1.);
111  feps += prim(n,k,j,i)*(1./pthermo->GetMassRatio(n) - 1.);
112  }
113  w_p = gm1*(u_e - KE)*feps/fsig;
114 
115  // apply pressure floor, correct total energy
116  u_e = (w_p > pressure_floor_) ? u_e : ((pressure_floor_/gm1)*fsig/feps + KE);
117  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
118 
119 #if DEBUG_LEVEL > 3
120  if (std::isnan(w_p) || (w_p < pressure_floor_)) {
121  msg << "### FATAL ERROR in function ConservedToPrimitive"
122  << std::endl << "Pressure reaches lowest value: " << w_p
123  << std::endl << "At position ("
124  << k << "," << j << "," << i << ") in rank " << Globals::my_rank <<
125  std::endl;
126  ATHENA_ERROR(msg);
127  }
128 #endif
129  }
130  }
131  }
132 
133 //#if DEBUG_LEVEL > 0
134 // Debugger *pdbg = pmy_block_->pdebug;
135 // pdbg = pdbg->StartTracking("EquationOfStates::ConservedToPrimitive");
136 // pdbg->Track3D("rho", IsPositive, prim, IDN);
137 // pdbg->Track3D("pres", IsPositive, prim, IPR);
138 //#endif
139 
140  return;
141 }
142 
143 //----------------------------------------------------------------------------------------
144 // \!fn void EquationOfState::PrimitiveToConserved(const AthenaArray<Real> &prim,
145 // const AthenaArray<Real> &bc, AthenaArray<Real> &cons, Coordinates *pco,
146 // int il, int iu, int jl, int ju, int kl, int ku);
147 // \brief Converts primitive variables into conservative variables
148 
150  const AthenaArray<Real> &prim, const AthenaArray<Real> &bc,
151  AthenaArray<Real> &cons, Coordinates *pco,
152  int il, int iu, int jl, int ju, int kl, int ku) {
153  Real igm1 = 1.0/(GetGamma() - 1.0);
154  Thermodynamics *pthermo = pmy_block_->pthermo;
155 
156  // Force outer-loop vectorization
157 #pragma omp simd
158  for (int k=kl; k<=ku; ++k) {
159  for (int j=jl; j<=ju; ++j) {
160  //#pragma omp simd
161 #pragma novector
162  for (int i=il; i<=iu; ++i) {
163  Real& u_d = cons(IDN,k,j,i);
164  Real& u_m1 = cons(IM1,k,j,i);
165  Real& u_m2 = cons(IM2,k,j,i);
166  Real& u_m3 = cons(IM3,k,j,i);
167  Real& u_e = cons(IEN,k,j,i);
168 
169  const Real& w_d = prim(IDN,k,j,i);
170  const Real& w_vx = prim(IVX,k,j,i);
171  const Real& w_vy = prim(IVY,k,j,i);
172  const Real& w_vz = prim(IVZ,k,j,i);
173  const Real& w_p = prim(IPR,k,j,i);
174 
175  // density
176  u_d = w_d;
177  for (int n = 1; n <= NVAPOR; ++n) {
178  cons(n,k,j,i) = prim(n,k,j,i)*w_d;
179  cons(IDN,k,j,i) -= cons(n,k,j,i);
180  }
181 
182  // momentum
183  u_m1 = w_vx*w_d;
184  u_m2 = w_vy*w_d;
185  u_m3 = w_vz*w_d;
186 
187  // total energy
188  Real KE = 0.5*w_d*(w_vx*w_vx + w_vy*w_vy + w_vz*w_vz);
189  Real fsig = 1., feps = 1.;
190  // vapors
191  for (int n = 1; n <= NVAPOR; ++n) {
192  fsig += prim(n,k,j,i)*(pthermo->GetCvRatio(n) - 1.);
193  feps += prim(n,k,j,i)*(1./pthermo->GetMassRatio(n) - 1.);
194  }
195  u_e = igm1*w_p*fsig/feps + KE ;
196  }
197  }
198  }
199 
200  return;
201 }
202 
203 //----------------------------------------------------------------------------------------
204 // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
205 // \brief returns adiabatic sound speed given vector of primitive variables
206 Real EquationOfState::SoundSpeed(const Real prim[NHYDRO]) {
207  Thermodynamics *pthermo = pmy_block_->pthermo;
208 
209  Real fsig = 1., feps = 1.;
210  for (int n = 1; n <= NVAPOR; ++n) {
211  fsig += prim[n]*(pthermo->GetCvRatio(n) - 1.);
212  feps += prim[n]*(1./pthermo->GetMassRatio(n) - 1.);
213  }
214 
215  return std::sqrt((1. + (gamma_ - 1)*feps/fsig)*prim[IPR]/prim[IDN]);
216 }
217 
218 //----------------------------------------------------------------------------------------
219 // \!fn void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real> &prim, int k, int j,
220 // =int i)
221 // \brief Apply density and pressure floors to reconstructed L/R cell interface states
223  Real& w_d = prim(IDN,i);
224  Real& w_p = prim(IPR,i);
225 
226  // apply (prim) density floor
227  w_d = (w_d > density_floor_) ? w_d : density_floor_;
228  // apply pressure floor
229  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
230 
231  return;
232 }
233 
234 //----------------------------------------------------------------------------------------
235 // \!fn void EquationOfState::ApplyPrimitiveConservedFloors(AthenaArray<Real> &prim,
236 // AthenaArray<Real> &cons, FaceField &b, int k, int j, int i) {
237 // \brief Apply pressure (prim) floor and correct energy (cons) (typically after W(U))
240  int k, int j, int i) {
241  Real gm1 = GetGamma() - 1.0;
242  Real& w_d = prim(IDN,k,j,i);
243  Real& w_p = prim(IPR,k,j,i);
244 
245  Real& u_d = cons(IDN,k,j,i);
246  Real& u_e = cons(IEN,k,j,i);
247  // apply (prim) density floor, without changing momentum or energy
248  w_d = (w_d > density_floor_) ? w_d : density_floor_;
249  // ensure cons density matches
250  u_d = w_d;
251 
252  Real e_k = 0.5*w_d*(SQR(prim(IVX,k,j,i)) + SQR(prim(IVY,k,j,i))
253  + SQR(prim(IVZ,k,j,i)));
254  // apply pressure floor, correct total energy
255  u_e = (w_p > pressure_floor_) ?
256  u_e : ((pressure_floor_/gm1) + e_k);
257  w_p = (w_p > pressure_floor_) ?
258  w_p : pressure_floor_;
259 
260  return;
261 }
@ 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
#define SQR(x)
Definition: cdisort.h:482
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
void VaporConcentrationLimiter(AthenaArray< Real > &u)
Hydro * phydro
Definition: mesh.hpp:128
Thermodynamics * pthermo
Definition: mesh.hpp:135
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
int my_rank
Definition: globals.cpp:23