Canoe Example#3：Bryan Moist Uprising Bubble

#include <athena/bvals/bvals.hpp>  // defines BoundaryBase, BoundaryValues classes used for setting BCs on all data
// application
#include <application/application.hpp> // app注册
#include <application/exceptions.hpp>  // 报错信息
// 定义了vapor,cloud,chemitry,tracer,particle等介质的索引和函数，如GetVaporId();
#include <index_map.hpp>  

// climath
#include <climath/core.h> //基础的数学函数和单位转换
#include <climath/interpolation.h> // 插值工具包，包括定位，线性内插，样条等方法；
#include <climath/root.hpp> 
// special includes
#include "bryan_vapor_functions.hpp"//包含了本个例使用的湿空气物理，计算饱和水汽压

int iH2O, iH2Oc;  // indices for water vapor and cloud
Real p0, grav;  //  global var declaration on surface pressure and the gravity

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  AllocateUserOutputVariables(7);
  SetUserOutputVariableName(0, "temp");
  SetUserOutputVariableName(1, "theta");
  SetUserOutputVariableName(2, "theta_v");
  SetUserOutputVariableName(3, "mse");
  SetUserOutputVariableName(4, "theta_e");
  SetUserOutputVariableName(5, "rh");
  SetUserOutputVariableName(6, "qtol");
}

void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
  auto pthermo = Thermodynamics::GetInstance();

  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j)
      for (int i = is; i <= ie; ++i) {
        user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
        user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, p0, k, j, i);
        // theta_v
        user_out_var(2, k, j, i) =
            user_out_var(1, j, i) * pthermo->RovRd(this, k, j, i); // Theta/Rd
        // mse
        user_out_var(3, k, j, i) =
            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
        // theta_e
        user_out_var(4, k, j, i) =
            pthermo->EquivalentPotentialTemp(this, p0, iH2O, k, j, i);
        user_out_var(5, k, j, i) =
            pthermo->RelativeHumidity(this, iH2O, k, j, i);
        // total mixing ratio
        auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, i);
        user_out_var(6, k, j, i) = air.w[iH2O] + air.c[iH2Oc];
      }
}

  Real *const w; //! hydro data
  Real *const c; //! cloud data
  Real *const q; //! chemistry data
  Real *const x; //! tracer data
  Real *const t; //! turbulence data
  Real const *d; //! particle data

  AirParcel &ToMassFraction();
  AirParcel &ToMassConcentration();
  AirParcel &ToMoleFraction();
  AirParcel &ToMoleConcentration();

void Mesh::InitUserMeshData(ParameterInput *pin) {
  auto pindex = IndexMap::GetInstance();

  grav = -pin->GetReal("hydro", "grav_acc1");
  p0 = pin->GetReal("problem", "p0");

  // index
  iH2O = pindex->GetVaporId("H2O");
  iH2Oc = pindex->GetCloudId("H2O(c)");
}

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  auto pthermo = Thermodynamics::GetInstance();

  Real Ps = p0;
  Real Ts = pin->GetReal("problem", "Ts");  // 289.85
  Real xc = pin->GetReal("problem", "xc");  // 10,000 m
  Real zc = pin->GetReal("problem", "zc");  // 2,000 m
  Real xr = pin->GetReal("problem", "xr");  // 2,000
  Real zr = pin->GetReal("problem", "zr");  // 2,000
  Real dT = pin->GetReal("problem", "dT");  // 1E5 pa
  Real qt = pin->GetReal("problem", "qt");  // 0.0196

  AirParcel air(AirParcel::Type::MassFrac); // 定义一个MassFrac为单位的空气块
  air.w[iH2O] = qt; // 设定比湿
  air.c[iH2Oc] = 0.; 

  air.ToMoleFraction();  // 转换为摩尔分数
  qt = air.w[iH2O];

  // construct a reversible adiabat
  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j) {
      air.w[iH2O] = qt;
      air.w[IPR] = Ps;
      air.w[IDN] = Ts;
      air.c[iH2Oc] = 0.;

      //这里先使用地表气压和地表温度插值出第一层（is）和地表(0)中间高度处的气温和气压，水汽
      pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
                           Thermodynamics::Method::ReversibleAdiabat, grav);
      // 循环空气柱，自下而上构建出一条湿绝热廓线
      for (int i = is; i <= ie; ++i) {
        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
        pthermo->Extrapolate(&air, pcoord->dx1f(i),
                             Thermodynamics::Method::ReversibleAdiabat, grav);
      }
    }

void Thermodynamics::Extrapolate	(	AirParcel * 	qfrac,
Real 	dzORdlnp,
Method 	method,
Real 	grav = 0.,
Real 	userp = 0. 
)		const
Construct an 1d atmosphere

Parameters
[in,out]	qfrac	mole fraction representation of air parcel
[in]	dzORdlnp	vertical grid spacing
[in]	method	choose from [reversible, pseudo, dry, isothermal]
[in]	grav	gravitational acceleration
[in]	userp	user parameter to adjust the temperature gradient

#include <algorithm>
void  Thermodynamics::rk4IntegrateLnp(AirParcel *air, Real dlnp, Method method,
                                      Real adlnTdlnP)

 Real Thermodynamics::calDlnTDlnP(AirParcel const& qfrac, Real latent[]) const {
   // calculate gammad
   Real gammad = GetGammad(qfrac);
  
   Real q_gas = 1.;
   for (int n = 0; n < NCLOUD; ++n) q_gas -= qfrac.c[n];
  
   Real f_sig = 1.;
   // vapor
   for (int n = 1; n <= NVAPOR; ++n)
     f_sig += qfrac.w[n] * (cp_ratio_mole_[n] - 1.);
   // cloud
   for (int n = 0; n < NCLOUD; ++n)
     f_sig += qfrac.c[n] * (cp_ratio_mole_[1 + NVAPOR + n] - 1.);
   Real cphat_ov_r = gammad / (gammad - 1.) * f_sig / q_gas;
  
   // vapor
   Real xd = q_gas;
   for (int n = 1; n <= NVAPOR; ++n) xd -= qfrac.w[n];
  
   Real c1 = 0., c2 = 0., c3 = 0.;
   for (int iv = 1; iv <= NVAPOR; ++iv) {
     c1 += qfrac.w[iv] / xd * latent[iv];
     c2 += qfrac.w[iv] / xd * latent[iv] * latent[iv];
     c3 += qfrac.w[iv] / xd;
   }
  
   return (1. + c1) / (cphat_ov_r + (c2 + c1 * c1) / (1. + c3));
 }

  // add temperature anomaly
  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j)
      for (int i = is; i <= ie; ++i) {
        Real x1 = pcoord->x1v(i);
        Real x2 = pcoord->x2v(j);
        Real L = sqrt(sqr((x2 - xc) / xr) + sqr((x1 - zc) / zr));

        if (L < 1.) {
          auto &&air = AirParcelHelper::gather_from_conserved(this, k, j, i);
          air.ToMoleFraction();
          Real temp_v = air.w[IDN] * pthermo->RovRd(air);   /// 计算虚温
          temp_v *= 1. + dT * sqr(cos(M_PI * L / 2.)) / 300.;   /// 叠加保护水汽的虚温

          Real temp;
          int err = root(air.w[IDN], air.w[IDN] + dT, 1.E-8, &temp,
                         [&pthermo, &air, temp_v](Real temp) {
                           air.w[IDN] = temp;

                           auto rates =
                               pthermo->TryEquilibriumTP_VaporCloud(air, iH2O);
                           air.w[iH2O] += rates[0];
                           air.c[iH2Oc] += rates[1];

                           return temp * pthermo->RovRd(air) - temp_v;
                         });

          if (err) throw RuntimeError("pgen", "TVSolver doesn't converge");

          air.w[IDN] = temp;
          auto rates = pthermo->TryEquilibriumTP_VaporCloud(air, iH2O);
          air.w[iH2O] += rates[0];
          air.c[iH2Oc] += rates[1];

          AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
        }
      }
}