Athena++/Atmosphere
Planetary Atmosphere Simulator
vapor_concentration_limiter.cpp
Go to the documentation of this file.
1 // Athena++ header
2 #include "../thermodynamics/thermodynamics.hpp"
3 #include "hydro.hpp"
4 
6 {
7  MeshBlock *pmb = pmy_block;
8  Thermodynamics *pthermo = pmb->pthermo;
9  if (NVAPOR == 0) return;
10 
11  int is = pmb->is; int js = pmb->js; int ks = pmb->ks;
12  int ie = pmb->ie; int je = pmb->je; int ke = pmb->ke;
13 
14  for (int k=ks; k<=ke; ++k)
15  for (int j=js; j<=je; ++j) {
16  // fix negative vapor from top down
17  for (int i = ie; i > is; --i) {
18  Real density = 0., v1, v2, v3;
19 
20  for (int n = 0; n <= NVAPOR; ++n) density += u(n,k,j,i);
21  v1 = u(IM1,k,j,i)/density,
22  v2 = u(IM2,k,j,i)/density,
23  v3 = u(IM3,k,j,i)/density;
24 
25  // internal energy
26  Real KE = 0.5*density*(v1*v1 + v2*v2 + v3*v3);
27  Real LE = 0., rhocv = 0.;
28  for (int n = 0; n <= NVAPOR; ++n)
29  rhocv += u(n,k,j,i)*pthermo->GetCv(n);
30  Real temp = (u(IEN,k,j,i) - KE)/rhocv;
31 
32  for (int n = 1; n <= NVAPOR; ++n) {
33  Real rho = u(n,k,j,i);
34  if (rho < 0.) {
35  u(n,k,j,i-1) += rho;
36  u(IM1,k,j,i-1) += v1*rho;
37  u(IM2,k,j,i-1) += v2*rho;
38  u(IM3,k,j,i-1) += v3*rho;
39  Real en = pthermo->GetCv(n)*temp + 0.5*(v1*v1 + v2*v2 + v3*v3);
40  u(IEN,k,j,i-1) += en*rho;
41 
42  u(n,k,j,i) = 0.;
43  u(IM1,k,j,i) -= v1*rho;
44  u(IM2,k,j,i) -= v2*rho;
45  u(IM3,k,j,i) -= v3*rho;
46  u(IEN,k,j,i) -= en*rho;
47  }
48  }
49  }
50  }
51 }
double Real
Definition: athena.hpp:29
@ IM3
Definition: athena.hpp:135
@ IEN
Definition: athena.hpp:135
@ IM1
Definition: athena.hpp:135
@ IM2
Definition: athena.hpp:135
MeshBlock * pmy_block
Definition: hydro.hpp:43
AthenaArray< Real > u
Definition: hydro.hpp:46
void VaporConcentrationLimiter(AthenaArray< Real > &u)
int ie
Definition: mesh.hpp:100
int js
Definition: mesh.hpp:100
int je
Definition: mesh.hpp:100
int ke
Definition: mesh.hpp:100
int is
Definition: mesh.hpp:100
int ks
Definition: mesh.hpp:100
Thermodynamics * pthermo
Definition: mesh.hpp:135
Real GetCv(int i) const
float temp
Definition: fit_ammonia.py:6