Canoe
Comprehensive Atmosphere N' Ocean Engine
Example #2 : The robert Problem

Table of contents

Introduction

The second example shows the the rising of a warm bubble, which to some extent, resembles the mushroom cloud of an atomic bomb. This test case was proposed by [1]. The original setting was is 2D and we extended it to 3D. The size of the domain is 1 km \(\times\) 1 km \(\times\) 2.5km. The model resolution is uniformlly 5 m in each dimension. The bubble is initially centered at \(x_0 = y_0 = 500 \text{m}, z_0 = 260 \text{m}\) with a Gaussian-shaped potential temperature anomaly as:

\begin{eqnarray*} \Delta \theta &=& 0.5\;\text{K}, \qquad\quad\quad r \leq a \\ &=& 0.5e^{-(r-a)^2/s^2}\;\text{K}, \quad r > a \end{eqnarray*}

where \(r\) is the distance to the center, \(a=50\) m and \(s=100\) m. The background atmosphere is isentropic with surface pressure being 1 bar and surface temperature 303.15 K. The bubble is then released and evolve to 1200 s.

The commented program

/* -------------------------------------------------------------------------------------
* SNAP Example Program
*
* Contributer:
* Cheng Li, University of Michigan
*
* Year: 2021
* Contact: chengcli@umich.edu
* Reference: Robert et al., 1992
* -------------------------------------------------------------------------------------
*/

Include files

These input files are just like those in the Example #1 : The straka Problem, so additional comments are not required.

#include <athena/athena.hpp>
#include <athena/athena_arrays.hpp>
#include <athena/coordinates/coordinates.hpp>
#include <athena/eos/eos.hpp>
#include <athena/field/field.hpp>
#include <athena/hydro/hydro.hpp>
#include <athena/mesh/mesh.hpp>
#include <athena/parameter_input.hpp>
#include <athena/stride_iterator.hpp>

canoe

#include <impl.hpp>

snap

#include <snap/thermodynamics/thermodynamics.hpp>

Preamble

We only need one global variables here, the surface pressure

Real p0;
p0
Real p0
Definition: robert.cpp:36

Same as that in Example #1 : The straka Problem, 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);
}
}
Thermodynamics::GetInstance
static Thermodynamics const * GetInstance()
Return a pointer to the one and only instance of Thermodynamics.
Definition: thermodynamics.cpp:100
make_plots.j
j
Definition: make_plots.py:26

Initialize surface pressure from input file.

void Mesh::InitUserMeshData(ParameterInput *pin) {
p0 = pin->GetReal("problem", "p0");
}

Initial condition

We do not need forcings other than gravity in this problem, so we go directly to the initial condition.

void MeshBlock::ProblemGenerator(ParameterInput *pin) {

Similar to Example #1 : The straka Problem, 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");
fit_ammonia.a
a
Definition: fit_ammonia.py:8
make_plots.Ts
int Ts
Definition: make_plots.py:17
Rd
Real Rd
Definition: straka.cpp:66
cp
Real cp
Definition: straka.cpp:66

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;
phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd);
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);
phydro->w(IVX, k, j, i) = phydro->w(IVY, k, j, i) = 0.;
}
fit_ammonia.r
r
Definition: fit_ammonia.py:8
fit_ammonia.temp
float temp
Definition: fit_ammonia.py:6
make_plots.x1
int x1
Definition: make_plots.py:10
make_plots.x2
int x2
Definition: make_plots.py:11

Change primitive variables to conserved variables

peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
js, je, ks, ke);
}

Results

Here we show the animation of the rising bubble

Animation of a rising air bubble

The temperature anomaly of the bubble is colored in blue/red with blue being colder and red being warmer. Initially, the warm part of the bubble is hiden at its center and you only see the colder outer part. Multiple convective overturnings develop seen during the rising motion to deliver the warm gas upward and outward. Unlike the 2D solution shown in [1], the 3D evolution of the bubble features a vortex ring surounding the core of the rising center and the head of the rising bubble is minimum compared to the ring.

The plain program

/* -------------------------------------------------------------------------------------
* SNAP Example Program
*
* Contributer:
* Cheng Li, University of Michigan
*
* Year: 2021
* Contact: chengcli@umich.edu
* Reference: Robert et al., 1992
* -------------------------------------------------------------------------------------
*/
#include <athena/athena.hpp>
#include <athena/athena_arrays.hpp>
#include <athena/coordinates/coordinates.hpp>
#include <athena/eos/eos.hpp>
#include <athena/field/field.hpp>
#include <athena/hydro/hydro.hpp>
#include <athena/mesh/mesh.hpp>
#include <athena/parameter_input.hpp>
#include <athena/stride_iterator.hpp>
#include <impl.hpp>
#include <snap/thermodynamics/thermodynamics.hpp>
Real p0;
void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
AllocateUserOutputVariables(2);
SetUserOutputVariableName(0, "temp");
SetUserOutputVariableName(1, "theta");
}
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);
}
}
void Mesh::InitUserMeshData(ParameterInput *pin) {
p0 = pin->GetReal("problem", "p0");
}
void MeshBlock::ProblemGenerator(ParameterInput *pin) {
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");
bool uniform_bubble =
pin->GetOrAddBoolean("problem", "uniform_bubble", false);
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;
phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd);
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);
phydro->w(IVX, k, j, i) = phydro->w(IVY, k, j, i) = 0.;
}
peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
js, je, ks, ke);
}