Canoe Example#2：Robert Uprisng Bubble

It has been 501 days since the last update, the content of the article may be outdated.

本节继续学习个例二，Robert problem。模拟了一个热气团的上升和充分发展的Kelvin–Helmholtz不稳定性，不考虑额外的Forcing。

问题描述

本示例展示了一个热空气泡从近地表上升的情景模拟。这个测试案例由André Robert，1993提出。本个例将其从2D扩展到3D模拟。模拟空间的尺寸为 $1000 m \times 1000 m \times 2000 m$ ，分辨率在每个维度上均匀为 $25 m$ 。
起始时刻，气泡中心最初位于 $x_{0} = y_{0} = 500 m$ ， $z_{0} = 260 m$ 处，其位温异常（Anomaly of Potential Temperature）呈高斯分布，表示为：
$\begin{array}{rcl} Δ θ & = & {\begin{cases} 5; K, r \leq a \\ 5 e^{- (r - a)^{2} / s^{2}}; K, r > a \end{cases} \end{array}$
其中 $r$ 为到气泡中心点的距离，参数 $a = 50 m ， s = 100 m$ 。背景大气是等熵的(绝热、无摩擦)，地表气压 $p_{0} = 1 b a r$ ，地表温度 $T_{s} = 303.15 K$ 。随后释放气泡，模拟900秒内的演变过程。

源码分析

本个例的和UOV的设置同Straka问题一样，不再赘述。

cpp

// We only need one global variables here, the surface pressure
Real p0;

// Same as that in @ref straka, make outputs of temperature and potential
// temperature.
void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  AllocateUserOutputVariables(2);
  SetUserOutputVariableName(0, "temp");
  SetUserOutputVariableName(1, "theta");
}

// Set temperature and potential temperature.
void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
  auto pthermo = Thermodynamics::GetInstance();

  Real gamma = peos->GetGamma();
  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);
      }
}

对Mesh类进行初始化，本个例不考虑额外的forcing，只有模式内部的重力场作用。
这里为什么只对Mesh初始化了 $P_{0}$ 这一个初始值呢？如果在这里先初始化Ts会怎样？必须要作为global参数吗？

cpp

// Initialize surface pressure from input file.
void Mesh::InitUserMeshData(ParameterInput *pin) {
  p0 = pin->GetReal("problem", "p0");
}
// We do not need forcings other than gravity in this problem,
// so we go directly to the initial condition.

对MeshBlock类进行初始化，遍历设置每个格点的初始条件：

cpp

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  // Similar to @ref straka, read variables in the input file
  Real gamma = pin->GetReal("hydro", "gamma");
  Real grav = -phydro->hsrc.GetG1();
  Real Ts = pin->GetReal("problem", "Ts");
  Real Rd = pin->GetReal("thermodynamics", "Rd");
  Real cp = gamma / (gamma - 1.) * Rd;

  Real xc = pin->GetReal("problem", "xc");
  Real yc = pin->GetReal("problem", "yc");
  Real zc = pin->GetReal("problem", "zc");
  Real s = pin->GetReal("problem", "s");
  Real a = pin->GetReal("problem", "a");
  Real dT = pin->GetReal("problem", "dT");

  // Whether to do a uniform bubble or not.
  bool uniform_bubble =
      pin->GetOrAddBoolean("problem", "uniform_bubble", false);

  // Loop over the grids and set initial condition
  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 x3 = pcoord->x3v(k);
        Real r = sqrt((x3 - yc) * (x3 - yc) + (x2 - xc) * (x2 - xc) +
                      (x1 - zc) * (x1 - zc));
        Real temp = Ts - grav * x1 / cp;  //  Adiabatic temperature gradient.
        phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd);  // $p=p_0(\frac{T}{T_s})^{c_p/R_d}$
        if (r <= a)
          temp += dT * pow(phydro->w(IPR, k, j, i) / p0, Rd / cp);
        else if (!uniform_bubble)
          temp += dT * exp(-(r - a) * (r - a) / (s * s)) *
                  pow(phydro->w(IPR, k, j, i) / p0, Rd / cp);
        phydro->w(IDN, k, j, i) = phydro->w(IPR, k, j, i) / (Rd * temp); // $ \pho=\frac{P}{Rd T}$
        phydro->w(IVX, k, j, i) = phydro->w(IVY, k, j, i) = 0.;
      }

  // Change primitive variables to conserved variables
  peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
                             js, je, ks, ke);
}

这里涉及到的物理有：

计算温度： $T = T_{s} - \frac{g h}{c_{p}}$
计算压强： $p = p_{0} (\frac{T}{T_{s}})^{c_{p} / R_{d}}$
计算温度异常量： $Δ T = Δ θ (\frac{p}{p_{0}})^{R d / c_{p}}$
计算密度： $ρ = \frac{P}{R d T}$

运行和结果

项目构建和编译的过程和Straka问题一样，这里不重复说了。这里有一个不同点，为了提高对3D模拟中众多网格的运算效率，使用了8核并行运算。

bash

$ mpiexec -n 8 ./robert.release -i robert3d.inp

这里为什么 np=8呢，是因为我们在输入文件中对mesh为维数和单核MeshBlock维数进行了设计：

yaml

# /robert3d.inp
<mesh>
nx1        = 80          
nx2        = 40          
nx3        = 40     
<meshblock>
nx1        = 40
nx2        = 20
nx3        = 20