Canoe
Comprehensive Atmosphere N' Ocean Engine
ideal_moist_hydro.cpp
Go to the documentation of this file.
1 // C/C++
2 #include <cmath> // sqrt()
3 #include <sstream>
4 #include <stdexcept>
5 
6 // application
7 #include <application/exceptions.hpp>
8 
9 // athena
10 #include <athena/athena.hpp>
11 #include <athena/athena_arrays.hpp>
12 #include <athena/eos/eos.hpp>
13 #include <athena/field/field.hpp>
14 #include <athena/globals.hpp>
15 #include <athena/hydro/hydro.hpp>
16 #include <athena/mesh/mesh.hpp>
17 #include <athena/parameter_input.hpp>
18 
19 // canoe
20 #include <configure.hpp>
21 #include <impl.hpp>
22 
23 // snap
24 #include "../thermodynamics/thermodynamics.hpp"
25 #include "eos_helper.hpp"
26 
27 #ifdef ENABLE_GLOG
28 #include <glog/logging.h>
29 #endif
30 
31 std::string str_grid_ij(AthenaArray<Real> const& var, int n, int k, int j,
32  int i, int il, int iu, int jl, int ju) {
33  std::stringstream msg;
34  for (int ii = std::max(i - 3, il); ii <= std::min(i + 3, iu); ++ii) {
35  msg << "i = " << ii << " ";
36  for (int jj = std::max(j - 3, jl); jj <= std::min(j + 3, ju); ++jj)
37  msg << var(n, k, jj, ii) << " ";
38  msg << std::endl;
39  }
40 
41  return msg.str();
42 }
43 
44 // EquationOfState constructor
45 
46 EquationOfState::EquationOfState(MeshBlock* pmb, ParameterInput* pin)
47  : pmy_block_(pmb),
48  gamma_{pin->GetReal("hydro", "gamma")},
49  density_floor_{
50  pin->GetOrAddReal("hydro", "dfloor", std::sqrt(1024 * float_min))},
51  pressure_floor_{
52  pin->GetOrAddReal("hydro", "pfloor", std::sqrt(1024 * float_min))},
53  scalar_floor_{
54  pin->GetOrAddReal("hydro", "sfloor", std::sqrt(1024 * float_min))} {}
55 
56 //----------------------------------------------------------------------------------------
57 // \!fn void EquationOfState::ConservedToPrimitive(AthenaArray<Real> &cons,
58 // const AthenaArray<Real> &prim_old, const FaceField &b,
59 // AthenaArray<Real> &prim, AthenaArray<Real> &bcc, Coordinates *pco,
60 // int il, int iu, int jl, int ju, int kl, int ku)
61 // \brief Converts conserved into primitive variables in adiabatic hydro.
62 
63 void EquationOfState::ConservedToPrimitive(
64  AthenaArray<Real>& cons, const AthenaArray<Real>& prim_old,
65  const FaceField& b, AthenaArray<Real>& prim, AthenaArray<Real>& bcc,
66  Coordinates* pco, int il, int iu, int jl, int ju, int kl, int ku) {
67  Real gm1 = GetGamma() - 1.0;
68  std::stringstream msg;
69  auto pthermo = Thermodynamics::GetInstance();
70 
71  apply_vapor_limiter(&cons, pmy_block_);
72 
73  for (int k = kl; k <= ku; ++k) {
74  for (int j = jl; j <= ju; ++j) {
75  for (int i = il; i <= iu; ++i) {
76  Real& u_d = cons(IDN, k, j, i);
77  Real& u_m1 = cons(IM1, k, j, i);
78  Real& u_m2 = cons(IM2, k, j, i);
79  Real& u_m3 = cons(IM3, k, j, i);
80  Real& u_e = cons(IEN, k, j, i);
81 
82  if (u_e < 0.) {
83  throw RuntimeError("ideal_moist_hydro", "negative energy");
84  }
85 
86  Real& w_d = prim(IDN, k, j, i);
87  Real& w_vx = prim(IVX, k, j, i);
88  Real& w_vy = prim(IVY, k, j, i);
89  Real& w_vz = prim(IVZ, k, j, i);
90  Real& w_p = prim(IPR, k, j, i);
91 
92  // apply density floor, without changing momentum or energy
93  u_d = (u_d > density_floor_) ? u_d : density_floor_;
94 
95  Real density = 0.;
96  for (int n = 0; n <= NVAPOR; ++n) density += cons(n, k, j, i);
97  w_d = density;
98  Real di = 1. / density;
99 
100  // mass mixing ratio
101  for (int n = 1; n <= NVAPOR; ++n)
102  prim(n, k, j, i) = cons(n, k, j, i) * di;
103 
104 #ifdef ENABLE_GLOG
105  LOG_IF(ERROR, std::isnan(w_d) || (w_d < density_floor_))
106  << "rank = " << Globals::my_rank << ", (k,j,i) = "
107  << "(" << k << "," << j << "," << i << ")"
108  << ", w_d = " << w_d << std::endl;
109 
110  LOG_IF(INFO, std::isnan(w_d) || (w_d < density_floor_))
111  << str_grid_ij(cons, IDN, k, j, i, il, iu, jl, ju);
112 #endif
113 
114  Real KE, fsig = 1., feps = 1.;
115  // vapors
116  for (int n = 1; n <= NVAPOR; ++n) {
117  fsig += prim(n, k, j, i) * (pthermo->GetCvRatioMass(n) - 1.);
118  feps += prim(n, k, j, i) * (pthermo->GetInvMuRatio(n) - 1.);
119  }
120 
121  int decay_factor = 1;
122  do {
123  u_m1 /= decay_factor;
124  u_m2 /= decay_factor;
125  u_m3 /= decay_factor;
126 
127  // Real di = 1.0/u_d;
128  w_vx = u_m1 * di;
129  w_vy = u_m2 * di;
130  w_vz = u_m3 * di;
131 
132  // internal energy
133  KE = 0.5 * di * (u_m1 * u_m1 + u_m2 * u_m2 + u_m3 * u_m3);
134  w_p = gm1 * (u_e - KE) * feps / fsig;
135 
136 #ifdef ENABLE_GLOG
137  LOG_IF(ERROR, w_p < 0.)
138  << "rank = " << Globals::my_rank << ", (k,j,i) = "
139  << "(" << k << "," << j << "," << i << ")"
140  << ", w_p = " << w_p << std::endl;
141 #endif
142  decay_factor = 2;
143  } while (w_p < 0.);
144 
145  // apply pressure floor, correct total energy
146  u_e = (w_p > pressure_floor_)
147  ? u_e
148  : ((pressure_floor_ / gm1) * fsig / feps + KE);
149  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
150  }
151  }
152  }
153 
154  // #if DEBUG_LEVEL > 0
155  // Debugger *pdbg = pmy_block_->pdebug;
156  // pdbg = pdbg->StartTracking("EquationOfStates::ConservedToPrimitive");
157  // pdbg->Track3D("rho", IsPositive, prim, IDN);
158  // pdbg->Track3D("pres", IsPositive, prim, IPR);
159  // #endif
160 
161  return;
162 }
163 
164 //----------------------------------------------------------------------------------------
165 // \!fn void EquationOfState::PrimitiveToConserved(const AthenaArray<Real>
166 // &prim,
167 // const AthenaArray<Real> &bc, AthenaArray<Real> &cons, Coordinates
168 // *pco, int il, int iu, int jl, int ju, int kl, int ku);
169 // \brief Converts primitive variables into conservative variables
170 
171 void EquationOfState::PrimitiveToConserved(const AthenaArray<Real>& prim,
172  const AthenaArray<Real>& bc,
173  AthenaArray<Real>& cons,
174  Coordinates* pco, int il, int iu,
175  int jl, int ju, int kl, int ku) {
176  Real igm1 = 1.0 / (GetGamma() - 1.0);
177  auto pthermo = Thermodynamics::GetInstance();
178 
179  for (int k = kl; k <= ku; ++k) {
180  for (int j = jl; j <= ju; ++j) {
181  for (int i = il; i <= iu; ++i) {
182  Real& u_d = cons(IDN, k, j, i);
183  Real& u_m1 = cons(IM1, k, j, i);
184  Real& u_m2 = cons(IM2, k, j, i);
185  Real& u_m3 = cons(IM3, k, j, i);
186  Real& u_e = cons(IEN, k, j, i);
187 
188  const Real& w_d = prim(IDN, k, j, i);
189  const Real& w_vx = prim(IVX, k, j, i);
190  const Real& w_vy = prim(IVY, k, j, i);
191  const Real& w_vz = prim(IVZ, k, j, i);
192  const Real& w_p = prim(IPR, k, j, i);
193 
194  // density
195  u_d = w_d;
196  for (int n = 1; n <= NVAPOR; ++n) {
197  cons(n, k, j, i) = prim(n, k, j, i) * w_d;
198  cons(IDN, k, j, i) -= cons(n, k, j, i);
199  }
200 
201  // momentum
202  u_m1 = w_vx * w_d;
203  u_m2 = w_vy * w_d;
204  u_m3 = w_vz * w_d;
205 
206  // total energy
207  Real KE = 0.5 * w_d * (w_vx * w_vx + w_vy * w_vy + w_vz * w_vz);
208  Real fsig = 1., feps = 1.;
209  // vapors
210  for (int n = 1; n <= NVAPOR; ++n) {
211  fsig += prim(n, k, j, i) * (pthermo->GetCvRatioMass(n) - 1.);
212  feps += prim(n, k, j, i) * (pthermo->GetInvMuRatio(n) - 1.);
213  }
214  u_e = igm1 * w_p * fsig / feps + KE;
215  }
216  }
217  }
218 
219  return;
220 }
221 
222 //----------------------------------------------------------------------------------------
223 // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
224 // \brief returns adiabatic sound speed given vector of primitive variables
225 Real EquationOfState::SoundSpeed(const Real prim[NHYDRO]) {
226  auto pthermo = Thermodynamics::GetInstance();
227 
228  Real fsig = 1., feps = 1.;
229  for (int n = 1; n <= NVAPOR; ++n) {
230  fsig += prim[n] * (pthermo->GetCvRatioMass(n) - 1.);
231  feps += prim[n] * (pthermo->GetInvMuRatio(n) - 1.);
232  }
233 
234  return std::sqrt((1. + (gamma_ - 1) * feps / fsig) * prim[IPR] / prim[IDN]);
235 }
236 
237 //----------------------------------------------------------------------------------------
238 // \!fn void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real> &prim, int
239 // k, int j,
240 // =int i)
241 // \brief Apply density and pressure floors to reconstructed L/R cell interface
242 // states
243 void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real>& prim, int k,
244  int j, int i) {
245  Real& w_d = prim(IDN, i);
246  Real& w_p = prim(IPR, i);
247 
248  // apply (prim) density floor
249  w_d = (w_d > density_floor_) ? w_d : density_floor_;
250  // apply pressure floor
251  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
252 
253  return;
254 }
255 
256 //----------------------------------------------------------------------------------------
257 // \!fn void EquationOfState::ApplyPrimitiveConservedFloors(AthenaArray<Real>
258 // &prim,
259 // AthenaArray<Real> &cons, FaceField &b, int k, int j, int i) {
260 // \brief Apply pressure (prim) floor and correct energy (cons) (typically after
261 // W(U))
262 void EquationOfState::ApplyPrimitiveConservedFloors(AthenaArray<Real>& prim,
263  AthenaArray<Real>& cons,
264  AthenaArray<Real>& bcc,
265  int k, int j, int i) {
266  Real gm1 = GetGamma() - 1.0;
267  Real& w_d = prim(IDN, k, j, i);
268  Real& w_p = prim(IPR, k, j, i);
269 
270  Real& u_d = cons(IDN, k, j, i);
271  Real& u_e = cons(IEN, k, j, i);
272  // apply (prim) density floor, without changing momentum or energy
273  w_d = (w_d > density_floor_) ? w_d : density_floor_;
274  // ensure cons density matches
275  u_d = w_d;
276 
277  Real e_k = 0.5 * w_d *
278  (SQR(prim(IVX, k, j, i)) + SQR(prim(IVY, k, j, i)) +
279  SQR(prim(IVZ, k, j, i)));
280  // apply pressure floor, correct total energy
281  u_e = (w_p > pressure_floor_) ? u_e : ((pressure_floor_ / gm1) + e_k);
282  w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
283 
284  return;
285 }
apply_vapor_limiter
void apply_vapor_limiter(AthenaArray< Real > *pu, MeshBlock *pmb)
Definition: apply_vapor_limiter.cpp:11
Thermodynamics::GetInstance
static Thermodynamics const * GetInstance()
Return a pointer to the one and only instance of Thermodynamics.
Definition: thermodynamics.cpp:100
min
double min(double x1, double x2, double x3)
Definition: core.h:16
max
double max(double x1, double x2, double x3)
Definition: core.h:19
eos_helper.hpp
str_grid_ij
std::string str_grid_ij(AthenaArray< Real > const &var, int n, int k, int j, int i, int il, int iu, int jl, int ju)
Definition: ideal_moist_hydro.cpp:31
fit_ammonia.b
b
Definition: fit_ammonia.py:8
make_plots.j
j
Definition: make_plots.py:26