Canoe
Comprehensive Atmosphere N' Ocean Engine
shallow_xy_hydro.cpp
Go to the documentation of this file.
1 // \brief implements functions in class EquationOfState for shallow water
3 // hydrodynamics
4 
5 // C/C++
6 #include <cmath> // sqrt()
7 
8 // athena
9 #include <athena/athena.hpp>
10 #include <athena/eos/eos.hpp>
11 #include <athena/field/field.hpp>
12 #include <athena/hydro/hydro.hpp>
13 #include <athena/mesh/mesh.hpp>
14 #include <athena/parameter_input.hpp>
15 
16 // EquationOfState constructor
17 //
18 EquationOfState::EquationOfState(MeshBlock *pmb, ParameterInput *pin)
19  : pmy_block_(pmb),
20  density_floor_{
21  pin->GetOrAddReal("hydro", "dfloor", std::sqrt(1024 * float_min))},
22  pressure_floor_{
23  pin->GetOrAddReal("hydro", "pfloor", std::sqrt(1024 * float_min))} {}
24 
25 //----------------------------------------------------------------------------------------
26 // \!fn void EquationOfState::ConservedToPrimitive(AthenaArray<Real> &cons,
27 // const AthenaArray<Real> &prim_old, const FaceField &b,
28 // AthenaArray<Real> &prim, AthenaArray<Real> &bcc, Coordinates *pco,
29 // int il, int iu, int jl, int ju, int kl, int ku)
30 // \brief Converts conserved into primitive variables in shallow water.
31 
32 void EquationOfState::ConservedToPrimitive(
33  AthenaArray<Real> &cons, const AthenaArray<Real> &prim_old,
34  const FaceField &b, AthenaArray<Real> &prim, AthenaArray<Real> &bcc,
35  Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku) {
36  for (int k = kl; k <= ku; ++k)
37  for (int j = jl; j <= ju; ++j)
38  for (int i = il; i <= iu; ++i) {
39  Real &ud = cons(IDN, k, j, i);
40  Real &wd = prim(IDN, k, j, i);
41 
42  // apply density floor, without changing momentum or energy
43  ud = (ud > density_floor_) ? ud : density_floor_;
44  wd = ud;
45 
46  Real di = 1. / ud;
47  prim(IVX, k, j, i) = cons(IM1, k, j, i) * di;
48  prim(IVY, k, j, i) = cons(IM2, k, j, i) * di;
49  prim(IVZ, k, j, i) = 0.;
50  }
51 }
52 
53 //----------------------------------------------------------------------------------------
54 // \!fn void EquationOfState::PrimitiveToConserved(const AthenaArray<Real>
55 // &prim,
56 // const AthenaArray<Real> &bc, AthenaArray<Real> &cons, Coordinates
57 // *pco, int il, int iu, int jl, int ju, int kl, int ku);
58 // \brief Converts primitive variables into conservative variables
59 
60 void EquationOfState::PrimitiveToConserved(const AthenaArray<Real> &prim,
61  const AthenaArray<Real> &bc,
62  AthenaArray<Real> &cons,
63  Coordinates *pco, int il, int iu,
64  int jl, int ju, int kl, int ku) {
65  for (int k = kl; k <= ku; ++k)
66  for (int j = jl; j <= ju; ++j)
67  for (int i = il; i <= iu; ++i) {
68  const Real &wd = prim(IDN, k, j, i);
69  cons(IDN, k, j, i) = wd;
70  cons(IM1, k, j, i) = prim(IVX, k, j, i) * wd;
71  cons(IM2, k, j, i) = prim(IVY, k, j, i) * wd;
72  cons(IM3, k, j, i) = 0.;
73  }
74 }
75 
76 //----------------------------------------------------------------------------------------
77 // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
78 // \brief returns shallow water gravity wave speed given vector of primitive
79 // variables
80 
81 Real EquationOfState::SoundSpeed(const Real prim[NHYDRO]) {
82  return std::sqrt(prim[IDN]);
83 }
84 
85 //---------------------------------------------------------------------------------------
86 // \!fn void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real> &prim,
87 // int k, int j, int i)
88 // \brief Apply density and pressure floors to reconstructed L/R cell interface
89 // states
90 
91 void EquationOfState::ApplyPrimitiveFloors(AthenaArray<Real> &prim, int k,
92  int j, int i) {
93  Real &w_d = prim(IDN, i);
94 
95  // apply density floor
96  w_d = (w_d > density_floor_) ? w_d : density_floor_;
97 
98  return;
99 }
100 
101 void EquationOfState::ApplyPrimitiveConservedFloors(AthenaArray<Real> &prim,
102  AthenaArray<Real> &cons,
103  AthenaArray<Real> &bcc,
104  int k, int j, int i) {
105  Real &w_d = prim(IDN, k, j, i);
106  Real &u_d = cons(IDN, k, j, i);
107 
108  // apply density floor
109  w_d = (w_d > density_floor_) ? w_d : density_floor_;
110  u_d = w_d;
111 }