本节继续学习个例二，Robert problem。模拟了一个热气团的上升和充分发展的Kelvin–Helmholtz不稳定性，不考虑额外的Forcing。
问题描述
本示例展示了一个热空气泡从近地表上升的情景模拟。这个测试案例由André Robert，1993提出。本个例将其从2D扩展到3D模拟。模拟空间的尺寸为 $1000m \times 1000m \times 2000m$ ，分辨率在每个维度上均匀为 $25m$。 起始时刻，气泡中心最初位于 $x_0=y_0=500m$ ，$z_0=260m$ 处，其位温异常（Anomaly of Potential Temperature）呈高斯分布，表示为：
$ \Delta \theta =5 \; \text{K}, \qquad\quad\quad r \leq a \\ \Delta \theta =5 e^{-(r-a)^2/s^2} \; \text{K} , \quad r > a $
其中 $r$ 为到气泡中心点的距离，参数 $a=50 m，s=100 m$ 。背景大气是等熵的(绝热、无摩擦)，地表气压 $p_0=1 bar$，地表温度 $T_s =303.15 K$ 。随后释放气泡，模拟900秒内的演变过程。
源码分析
本个例的<head>和UOV的设置同Straka问题一样，不再赘述。
// 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参数吗？
// 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类进行初始化，遍历设置每个格点的初始条件：
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{gh}{c_p}$
- 计算压强：$p=p_0(\frac{T}{T_s})^{c_p/R_d}$
- 计算温度异常量 ：$ \Delta T = \Delta\theta (\frac{p}{p_0})^{Rd/c_p}$
- 计算密度：$ \rho=\frac{P}{Rd T}$
运行和结果
项目构建和编译的过程和Straka问题一样，这里不重复说了。这里有一个不同点，为了提高对3D模拟中众多网格的运算效率，使用了8核并行运算。
$ mpiexec -n 8 ./robert.release -i robert3d.inp
这里为什么 np=8呢，是因为我们在输入文件中对mesh为维数和单核MeshBlock维数进行了设计：
# /robert3d.inp
<mesh>
nx1 = 80
nx2 = 40
nx3 = 40
<meshblock>
nx1 = 40
nx2 = 20
nx3 = 20
也就是总共需要处理的网格数是 $80 \times 40 \times 40$， 是单核处理个数 $40 \times 20 \times 20$ 的八倍。 如果单核运行，会出现无法同步的问题，导致报错没有设定边界条件。
温度场、位温场的的演变：
剖面位置$x3=500m$ 处
问题
- 如果改变P0，Ts的初始化位置，会分别发生什么问题？
- 如果加上Straka问题中的Fusion作为Forcing，模拟结果会如何？
- 改变dT的大小，结果变化如何？
- 湍流场对Scheme的依赖性如何？