Athena++/Atmosphere Planetary Atmosphere Simulator
Example #1 : The 2d.straka Problem

# Introduction

This example shows the simulation of a sinking air bubble in two dimensions. It is proposed by [2] for benchmarking the numerical methods. The temperature anomaly of the cold air bubble is described by:

\begin{equation*} \Delta T = -15\;\text{K} \frac{\cos(\pi L)+1}{2}, \quad L < 1, \end{equation*}

where $$L$$ is the normalized distance to the center of the cold air bubble at $$(x_c,z_c)$$.

\begin{equation*} L = \Big(\big(\frac{x - x_c}{x_r}\big)^2+\big(\frac{z - z_c}{z_r}\big)^2\Big)^{1/2}. \end{equation*}

Particularly, $$x_c=0$$ km, $$z_c=3$$ km, $$x_r=4$$ km and $$z_r=2$$ km. The domain is two-dimensional and is 6.4 km tall and 25.6 km wide. Either the horizontal or the vertical resolution is 100 m. The background atmosphere is isentropic, where the temperature is 300 K at 1 bar surface pressure. Initially, the bubble is aloft in the air, hydrostatically balanced. Convective instability causes the bubble to sink while trigging horizontal and vertical velocity fields. The resulting density current propagates out at the surface developing Kelvin-Helmholtz instability. The purpose of the program is to resolve the instability and calculate the evolution of the density current.

### Patch the code

Since this example is the first of the series. We would spend a few paragraphs on how to configure, compile and run the code. The first step of running the code is to patch the code. The Athena++/Atmosphere model is developed based on the Athena++ project. We have modified many source codes and have added support for atmospheric simulation. In order to keep track of our modifications and keep it seperate from the main stream development of the Athena++ code, all source files that are in conflict with the original Athena++ code are placed in a dedicated folder drum. Then we patch the code using a python script:

./patch.py

After this script was executed, a file named patch_files will appear containing replacement files and supplementary files. Then the program will configure according to the rules specified by the patch_files.

### Configure the program

Next, configure the program using the following command.

./configure.py --prob=straka --nghost=3 -netcdf

straka is the name of the application. Internally, the configure script searches for file named straka.cpp under the directory src/pgen. It is the problem generator file that will be complied in the next section. 3 is the number of grids in the ghost cells outside of the computational domain. A higher order numerical scheme requires a larger number of ghost cells. We are using the Weighted Essentially Non-Oscillatory (WENO) method, which requires 3 ghost celss. Lastly, we enable the NetCDF output by the option -netcdf. If successful, a messge will appear saying:

Your Athena++ distribution has now been configured with the following options:
Problem generator:          straka
Coordinate system:          cartesian
X1 Grid ratio:              1.0
Ammonia vapor id:           -1
Water vapor id:             -1
Riemann solver:             hllc
Chemistry:                  chemistry
Forcing Jacobian:           JacobianGravityCoriolis
Magnetic fields:            OFF
Number of vapors:           0
Number of phases:           1
Number of scalars:          0
Special relativity:         OFF
General relativity:         OFF
Frame transformations:      OFF
Self-Gravity:               OFF
Super-Time-Stepping:        OFF
Shearing Box BCs:           OFF
Debug level:                0
Code coverage flags:        OFF
Floating-point precision:   double
Number of ghost cells:      3
MPI parallelism:            OFF
OpenMP parallelism:         OFF
FFT:                        OFF
HDF5 output:                OFF
NETCDF output:              ON
PNETCDF output:             OFF
Compiler:                   g++
Compilation command:        g++  -O3 -std=c++11


### Compile the program

Then we compile the program by:

make

It is possible to speed up the compilation using multiple core available to the computer. To do it, we use:

make -j4

in which the option -j4 suggests that using 4 cores to compile the program simultaneously. If the compilation is successfull, an executable file named "straka.ex" will appear in the bin folder.

### Run the program

Change directory to the running directory

cd example/2d.straka

It is recommended to make a soft link of the executable file to the running directory:

ln -s ../../bin/straka.ex ./

so that everything can be perfrom under the running directory. Finally, execute the program by

./straka.ex -i straka.inp > log.straka &

The straka.inp is the input file and the log file for the run is log.straka. The & symbol at the end lets the program to run in the background so that you can do other things. When the program finishes, many .nc files will appear. They record the dynamic fields at specified output frequency in the input file. The next step is to combine them into a single one. This functionality is provided by the script combine.py located in the root directory.

../../combine.py -o test

Now you will have a single output named straka-test-main.nc. Use the following command to check out the fields stored in the netcdf file :

ncdump -h straka-test-main.nc
int main(int argc, char *argv[])
Definition: main.cpp:63

The output is:

netcdf output {
dimensions:
time = UNLIMITED ; // (4 currently)
x1 = 64 ;
x1f = 65 ;
x2 = 256 ;
x2f = 257 ;
x3 = 1 ;
variables:
float time(time) ;
time:axis = "T" ;
time:units = "s" ;
float x1(x1) ;
x1:axis = "Z" ;
x1:units = "m" ;
float x1f(x1f) ;
float x2(x2) ;
x2:axis = "Y" ;
x2:units = "m" ;
float x2f(x2f) ;
float x3(x3) ;
x3:axis = "X" ;
x3:units = "m" ;
float rho(time, x1, x2, x3) ;
rho:units = "kg/m^3" ;
rho:long_name = "density" ;
float press(time, x1, x2, x3) ;
press:units = "pa" ;
press:long_name = "pressure" ;
float vel1(time, x1, x2, x3) ;
vel1:units = "m/s" ;
vel1:long_name = "velocity" ;
float vel2(time, x1, x2, x3) ;
vel2:units = "m/s" ;
vel2:long_name = "velocity" ;
float vel3(time, x1, x2, x3) ;
vel3:units = "m/s" ;
vel3:long_name = "velocity" ;
float temp(time, x1, x2, x3) ;
temp:units = "K" ;
temp:long_name = "temperature" ;
float theta(time, x1, x2, x3) ;
theta:units = "K" ;
theta:long_name = "potential temperature" ;
}


The variables are self-descriptive with units and long name, which is the advantage of using NetCDF as the output format!

# The commented program

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

### Include files

First we include some Athena++ headers files. They contain declarations of various components of the solver system. Athena++ is able to run both single precision (float) and double precision (double) applications. A macro Real is defined to indicate the precision and it is define in the following header file.

#include "../athena.hpp"

Then the basic multi-dimension array that stores fluid dynamic data is declared in the AthenaArray class.

#include "../athena_arrays.hpp"

The model executes according to the parameters specified in an input file. The input file and the parameters within are managed by the ParameterInput class.

#include "../parameter_input.hpp"

Coordinate related information are stored in the Coordinates class.

#include "../coordinates/coordinates.hpp"

Everying regarding the equation of state is treated by the EquationOfState class and declared in the following file.

#include "../eos/eos.hpp"

Evolution of hydrodynamic fields like pressure, density and velocities are managed by the Hydro class.

#include "../hydro/hydro.hpp"

Athena++ can alo simulate magnetohydrodynamics (MHD). However, since this example is a hydro-only application, the MHD file is only included as a place holder.

#include "../field/field.hpp"

Dynamic fields are evolving on a mesh and this is the file working with meshes and the partition of meshes among processors

#include "../mesh/mesh.hpp"

Finally, the Thermodynamics class works with thermodynamic aspects of the problem such as the temperature, potential temperature, condensation of vapor, etc.

#include "../thermodynamics/thermodynamic_funcs.hpp"

Functions in the math library are protected by a specific namespace because math functions are usually concise in the names, such as min, sqr. It is better to protect them in a namespace to avoid conflicts.

namespace math {
#include "../math/core.h"
}

Now, we get to the real program. First, we define several global variables specific to this appliciation For example, here K is the kinematic visocisty, p0 is the surface pressure, cp is the specific heat capacity, and Rd is the ideal gas constant of dry air. The Prandtl number is 1.

Real K, p0, cp, Rd;
double Real
Definition: athena.hpp:29

### User-defined output variables

The hydrodynamic solver evolves density, pressure and velocities with time, meaning that the (potential) temperature is a diagnostic quantity for output only. This block of code allocates the memory for the outputs of temperature and potential temperature.

{
SetUserOutputVariableName(0, "temp", "temperature", "K");
SetUserOutputVariableName(1, "theta", "potential temperature", "K");
}
void AllocateUserOutputVariables(int n)
Definition: meshblock.cpp:504
void InitUserMeshBlockData(ParameterInput *pin)
void SetUserOutputVariableName(int n, const char *name, const char *long_name="", const char *units="")
Definition: meshblock.cpp:526

Since temperature and potential temperature are only used for output purpose, they do not need to be calculated every dynamic time step. The subroutine below loops over all grids and calculates the value of temperature and potential temperature before the output time step. Particularly, the pointer to the Thermodynamics class pthermo calculates the temperature and potential temperature using its own member function Thermodynamics::GetTemp and PotentialTemp. pthermo, is a member of a higher level management class MeshBlock. So, you can use it directly inside a member function of class MeshBlock. In fact, all physics modules are managed by MeshBlock. As long as you have a pointer to a MeshBlock, you can access all physics in the simulation.

{
void UserWorkBeforeOutput(ParameterInput *pin)

Loop over the grids. js,je,is,ie are members of the MeshBlock class. They are integer values representing the start index and the end index of the grids in each dimension.

for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {

user_out_var stores the actual data. phydro is a pointer to the Hydro class, which has a member w that stores density, pressure, and velocities at each grid.

user_out_var(0,j,i) = pthermo->GetTemp(phydro->w.at(j,i));
user_out_var(1,j,i) = PotentialTemp(phydro->w.at(j,i), p0, pthermo);
}
}
Real PotentialTemp(T w, Real p0, Thermodynamics *pthermo)
Potential temperature.

### Forcing function

The sinking bubble is forced by the gravity and the viscous and thermal dissipation. The gravitational forcing is taken care of by other parts of the program and we only need to write a few lines of code to facilitate dissipation. The first step is to define a function that takes a pointer to the MeshBlock as an argument such that we can access all physics via this pointer. The name of this function is not important. It can be anything. But the order and types of the arguments must be (MeshBlock *, Real const, Real const, AthenaArray<Real> const&, AthenaArray<Real> const&, AthenaArray<Real> &). They are called the signature of the function.

void Diffusion(MeshBlock *pmb, Real const time, Real const dt,
{

pcoord is a pointer to the Coordinates class and it is a member of the MeshBlock class. We use the pointer to the MeshBlock class, pmb to access pcoord and use its member function to get the spacing of the grid.

Real dx = pmb->pcoord->dx1f(pmb->is);
Real dy = pmb->pcoord->dx2f(pmb->js);
AthenaArray< Real > dx1f
Definition: coordinates.hpp:41
AthenaArray< Real > dx2f
Definition: coordinates.hpp:41
int js
Definition: mesh.hpp:100
int is
Definition: mesh.hpp:100
Coordinates * pcoord
Definition: mesh.hpp:122

Similarly, we use pmb to find the pointer to the Thermodynamics class, pthermo.

Thermodynamics *pthermo = pmb->pthermo;
Thermodynamics * pthermo
Definition: mesh.hpp:135

Loop over the grids.

for (int j = pmb->js; j <= pmb->je; ++j)
for (int i = pmb->is; i <= pmb->ie; ++i) {

Similar to what we have done in MeshBlock::UserWorkBeforeOutput, we use the Thermodynamics class to calculate temperature and potential temperature.

Real temp = pthermo->GetTemp(w.at(j,i));
Real theta = PotentialTemp(w.at(j,i), p0, pthermo);
StrideIterator< T * > at(int i) const
Real GetTemp(T w) const
Temperature.
float temp
Definition: fit_ammonia.py:6

The thermal diffusion is applied to the potential temperature field, which is not exactly correct. But this is the setting of the test program.

Real theta_ip1_j = PotentialTemp(w.at(j+1,i), p0, pthermo);
Real theta_im1_j = PotentialTemp(w.at(j-1,i), p0, pthermo);
Real theta_i_jp1 = PotentialTemp(w.at(j,i+1), p0, pthermo);
Real theta_i_jm1 = PotentialTemp(w.at(j,i-1), p0, pthermo);

Add viscous dissipation to the velocities. Now you encounter another variable called u. u stands for conserved variables, which are density, momentums and total energy (kinetic + internal). In contrast, w stands for primitive variables, which are density, velocities, and pressure. Their relation is handled by the EquationOfState class. The conserved variables are meant for internal calculation. Solving for conserved variables guarantees the conservation properties. However, conserved variables are not easy to use for diagnostics. Therefore, another group of variables called the primitive variables are introduced for external calculations, like calculating transport fluxes, radiative fluxes and interacting with other physical components of the system. In this case, the diffusion is calculated by using the primitive variabes and the result is added to the conserved variables to ensure conservation properties.

u(IM1,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*(
w(IV1,j,i-1) + w(IV1,j,i+1) + w(IV1,j-1,i) + w(IV1,j+1,i) - 4.*w(IV1,j,i));
u(IM2,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*(
w(IV2,j,i-1) + w(IV2,j,i+1) + w(IV2,j-1,i) + w(IV2,j+1,i) - 4.*w(IV2,j,i));
@ IV2
Definition: athena.hpp:140
@ IV1
Definition: athena.hpp:140
@ IM1
Definition: athena.hpp:135
@ IDN
Definition: athena.hpp:135
@ IM2
Definition: athena.hpp:135

u(IEN,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*cp*temp/theta*(theta_ip1_j + theta_im1_j +
theta_i_jp1 + theta_i_jm1 - 4.*theta);
}
}
@ IEN
Definition: athena.hpp:135

### Program specific variables

This is the place where program specific variables are initialized. Not that the function is a member function of the Mesh class rather than the MeshBlock class we have been working with. The difference between class Mesh and class MeshBlock is that class Mesh is an all-encompassing class that manages multiple MeshBlocks while class MeshBlock manages all physics modules. During the instantiation of the classes. class Mesh is instantiated first and then it instantiates all MeshBlocks inside it. Therefore, this subroutine runs before any MeshBlock.

{
void InitUserMeshData(ParameterInput *pin)

The program specific forcing parameters are set here. They come from the input file, which is parsed by the ParameterInput class

Real gamma = pin->GetReal("hydro", "gamma");
K = pin->GetReal("problem", "K");
p0 = pin->GetReal("problem", "p0");
Rd = pin->GetReal("thermodynamics", "Rd");
cp = gamma/(gamma - 1.)*Rd;
Real GetReal(std::string block, std::string name)

This line code enrolls the forcing function we wrote in section Forcing function

EnrollUserExplicitSourceFunction(Diffusion);
}

### Initial condition

This is the final part of the program, in which we set the initial condition. It is called the problem generator in a general Athena++ application.

{
void ProblemGenerator(ParameterInput *pin)

Get the value of gravity. Positive is upward and negative is downward. We take the negative to get the absolute value.

Real grav = -phydro->hsrc.GetG1();

These lines read the parameter values in the input file, which are organized into sections. The problem specific parameters are usually placed under section problem.

Real Ts = pin->GetReal("problem", "Ts");
Real xc = pin->GetReal("problem", "xc");
Real xr = pin->GetReal("problem", "xr");
Real zc = pin->GetReal("problem", "zc");
Real zr = pin->GetReal("problem", "zr");
Real dT = pin->GetReal("problem", "dT");

Loop over the grids. The purpose is to set the temperature, pressure and density fields at each cell-centered grid.

for (int j = js; j <= je; ++j) {
for (int i = is; i <= ie; ++i) {

Get the Cartesian coordinates in the vertical and horizontal directions. x1 is usually the vertical direction and x2 is the horizontal direction. The meaning of x1 and x2 may change with the coordinate system. x1v and x2v retrieve the coordinate at cell centers.

Real x1 = pcoord->x1v(i);
Real x2 = pcoord->x2v(j);

Distance to the center of a cold air bubble at (xc,zc).

Real L = sqrt(math::sqr((x2 - xc)/xr) + math::sqr((x1 - zc)/zr));
double sqr(double x)
Definition: core.h:9

Real temp = Ts - grav*x1/cp;

Once we know the temperature, we can calculate the adiabatic and hydrostatic pressure as $$p=p_0(\frac{T}{T_s})^{c_p/R_d}$$, where $$p_0$$ is the surface pressure and $$T_s$$ is the surface temperature.

phydro->w(IPR,j,i) = p0*pow(temp/Ts, cp/Rd);
if (L <= 1.)
temp += dT*(cos(M_PI*L) + 1.)/2.;
@ IPR
Definition: athena.hpp:139
#define M_PI
Definition: cdisort.h:470

Set density using ideal gas law.

phydro->w(IDN,j,i) = phydro->w(IPR,j,i)/(Rd*temp);

Initialize velocities to zero.

phydro->w(IV1,j,i) = phydro->w(IV2,j,i) = 0.;
}
}

We have set all primitive variables. The last step is to calculate the conserved variables based on the primitive variables. This is done by calling the member function EquationOfState::PrimitiveToConserved. Note that pfield->bcc is a placeholder because there is not MHD involved in this example.

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

# Results

We are going to make a plot of the result using python. The program is make_plots.py. Execute it by:

python3 make_plots.py

A file named 2d.straka-theta.png will be generated. It shows the potential temperature anomaly of the environment with respect to 300 K at four time steps. The maximum potential temperature anomaly should be zero because the motion is adiabatic and mimum potential temperature anomaly increases with time due to diffusion. The cold air spreads out as expected and develops rolling vortices. The result can be compared to Figure 3 of [2].

# The plain program

/* -------------------------------------------------------------------------------------
* SNAP Example Program
*
* Contributer: Cheng Li, University of Michigan, 2021
* Contact: chengcli@umich.edu
* Reference: Straka et al., 1993
* -------------------------------------------------------------------------------------
*/
#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "../parameter_input.hpp"
#include "../coordinates/coordinates.hpp"
#include "../eos/eos.hpp"
#include "../hydro/hydro.hpp"
#include "../field/field.hpp"
#include "../mesh/mesh.hpp"
#include "../thermodynamics/thermodynamic_funcs.hpp"
namespace math {
#include "../math/core.h"
}
Real K, p0, cp, Rd;
{
SetUserOutputVariableName(0, "temp", "temperature", "K");
SetUserOutputVariableName(1, "theta", "potential temperature", "K");
}
{
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
user_out_var(0,j,i) = pthermo->GetTemp(phydro->w.at(j,i));
user_out_var(1,j,i) = PotentialTemp(phydro->w.at(j,i), p0, pthermo);
}
}
void Diffusion(MeshBlock *pmb, Real const time, Real const dt,
{
Real dx = pmb->pcoord->dx1f(pmb->is);
Real dy = pmb->pcoord->dx2f(pmb->js);
Thermodynamics *pthermo = pmb->pthermo;
for (int j = pmb->js; j <= pmb->je; ++j)
for (int i = pmb->is; i <= pmb->ie; ++i) {
Real temp = pthermo->GetTemp(w.at(j,i));
Real theta = PotentialTemp(w.at(j,i), p0, pthermo);
Real theta_ip1_j = PotentialTemp(w.at(j+1,i), p0, pthermo);
Real theta_im1_j = PotentialTemp(w.at(j-1,i), p0, pthermo);
Real theta_i_jp1 = PotentialTemp(w.at(j,i+1), p0, pthermo);
Real theta_i_jm1 = PotentialTemp(w.at(j,i-1), p0, pthermo);
u(IM1,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*(
w(IV1,j,i-1) + w(IV1,j,i+1) + w(IV1,j-1,i) + w(IV1,j+1,i) - 4.*w(IV1,j,i));
u(IM2,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*(
w(IV2,j,i-1) + w(IV2,j,i+1) + w(IV2,j-1,i) + w(IV2,j+1,i) - 4.*w(IV2,j,i));
u(IEN,j,i) += dt*K*w(IDN,j,i)/(dx*dy)*cp*temp/theta*(theta_ip1_j + theta_im1_j +
theta_i_jp1 + theta_i_jm1 - 4.*theta);
}
}
{
Real gamma = pin->GetReal("hydro", "gamma");
K = pin->GetReal("problem", "K");
p0 = pin->GetReal("problem", "p0");
Rd = pin->GetReal("thermodynamics", "Rd");
cp = gamma/(gamma - 1.)*Rd;
}
{
Real grav = -phydro->hsrc.GetG1();
Real Ts = pin->GetReal("problem", "Ts");
Real xc = pin->GetReal("problem", "xc");
Real xr = pin->GetReal("problem", "xr");
Real zc = pin->GetReal("problem", "zc");
Real zr = pin->GetReal("problem", "zr");
Real dT = pin->GetReal("problem", "dT");
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(math::sqr((x2 - xc)/xr) + math::sqr((x1 - zc)/zr));
Real temp = Ts - grav*x1/cp;
phydro->w(IPR,j,i) = p0*pow(temp/Ts, cp/Rd);
if (L <= 1.)
temp += dT*(cos(M_PI*L) + 1.)/2.;
phydro->w(IDN,j,i) = phydro->w(IPR,j,i)/(Rd*temp);
phydro->w(IV1,j,i) = phydro->w(IV2,j,i) = 0.;
}
}
}
AthenaArray< Real > x1v
Definition: coordinates.hpp:42
AthenaArray< Real > x2v
Definition: coordinates.hpp:42
void PrimitiveToConserved(const AthenaArray< Real > &prim, const AthenaArray< Real > &bc, AthenaArray< Real > &cons, Coordinates *pco, int il, int iu, int jl, int ju, int kl, int ku)
Real GetG1() const
HydroSourceTerms hsrc
Definition: hydro.hpp:71
AthenaArray< Real > w
Definition: hydro.hpp:46
AthenaArray< Real > u
Definition: hydro.hpp:46
Hydro * phydro
Definition: mesh.hpp:128
EquationOfState * peos
Definition: mesh.hpp:133
int ie
Definition: mesh.hpp:100
AthenaArray< Real > user_out_var
Definition: mesh.hpp:112
int je
Definition: mesh.hpp:100
int ke
Definition: mesh.hpp:100
int ks
Definition: mesh.hpp:100
Field * pfield
Definition: mesh.hpp:129
void EnrollUserExplicitSourceFunction(SrcTermFunc my_func)
Definition: mesh.cpp:1198