athenapk
athenapk copied to clipboard
Published
20 hours ago •
parthenon-hpc-lab
parthenon-hpc-lab
Doubt in problem output (Not generating the output)
I am facing problem with output of the Rayleigh-Taylor problem.
C++ Code
//! \file rt.cpp
//! \brief Rayleigh Taylor problem generator
//! Problem domain should be -.05 < x < .05; -.05 < y < .05, -.1 < z < .1, gamma=5/3 to
//! match Dimonte et al. Interface is at z=0; perturbation added to Vz. Gravity acts in
//! z-dirn. Special reflecting boundary conditions added in x3. A=1/2. Options:
//! - iprob = 1 -- Perturb V3 using single mode
//! - iprob = 2 -- Perturb V3 using multiple mode
//! - iprob = 3 -- B rotated by "angle" at interface, multimode perturbation
//!
// C headers
// C++ headers
#include <algorithm> // min, max
#include <cmath> // log
#include <cstring> // strcmp()
#include <sstream>
// Parthenon headers
#include "mesh/mesh.hpp"
#include "parthenon/driver.hpp"
#include "parthenon/package.hpp"
#include "utils/error_checking.hpp"
// AthenaPK headers
#include "../main.hpp"
namespace rt {
using namespace parthenon::driver::prelude;
void SetupSingleMode(MeshBlock *pmb, parthenon::ParameterInput *pin) {
if (pmb->pmy_mesh->ndim == 1) {
PARTHENON_THROW("This problem should be either in 2d or 3d.");
return;
}
Real kx = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min);
Real ky = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min);
Real kz = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min);
Real amp = pin->GetReal("problem/rt","amp");
Real drat = pin->GetOrAddReal("problem/rt","drat",3.0);
Real grav_acc = pin->GetReal("hydro","const_accel_val");
auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
auto gam = pin->GetReal("hydro", "gamma");
auto gm1 = (gam - 1.0);
Real p0 = 1.0/gam;
// initialize conserved variables
auto &rc = pmb->meshblock_data.Get();
auto &u_dev = rc->Get("cons").data;
auto &coords = pmb->coords;
// initializing on host
auto u = u_dev.GetHostMirrorAndCopy();
for (size_t k = kb.s; k < kb.e; k++) {
for (size_t j = jb.s; j < jb.e; j++) {
for (size_t i = ib.s; j < ib.e; i++) {
auto x1v = coords.Xc<1>(i);
auto x2v = coords.Xc<2>(j);
auto x3v = coords.Xc<3>(k);
switch (pmb->pmy_mesh->ndim) {
case 2:{
Real den=1.0;
if (x2v > 0.0) den *= drat;
u(IM2,k,j,i) = (1.0 + cos(kx*x1v))*(1.0 + cos(ky*x2v))/4.0;
u(IDN,k,j,i) = den;
u(IM1,k,j,i) = 0.0;
u(IM2,k,j,i) *= (den*amp);
u(IM3,k,j,i) = 0.0;
u(IEN,k,j,i) = (p0 + grav_acc*den*x2v)/gm1 + 0.5*SQR(u(IM2,k,j,i))/den;
}
break;
case 3: {
Real den=1.0;
if (x3v > 0.0) den *= drat;
u(IM3,k,j,i) = (1.0+cos(kx*x1v))*(1.0+cos(ky*x2v))*(1.0+cos(kz*x3v))/8.0;
u(IDN,k,j,i) = den;
u(IM1,k,j,i) = 0.0;
u(IM2,k,j,i) = 0.0;
u(IM3,k,j,i) *= (den*amp);
u(IEN,k,j,i) = (p0 + grav_acc*den*x3v)/gm1 + 0.5*SQR(u(IM3,k,j,i))/den;
}
break;
}
}
}
}
}
void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
auto iprob = pin->GetOrAddInteger("problem/rt", "iprob", 1);
switch (iprob) {
case 1:
SetupSingleMode(pmb, pin);
break;
default:
std::stringstream msg;
msg << "Problem type " << iprob << " is not supported.";
PARTHENON_THROW(msg.str());
}
}
} // namespace rt
Input file
<comment>
problem = Kelvin-Helmholtz instability
reference = Lecoanet et al., MNRAS 455, 4274-4288, 2016
<job>
problem_id = rt
<problem/rt>
iprob = 1
amp = 0.01
drat = 2.0
<parthenon/mesh>
refinement = adaptive
numlevel = 3
nghost = 4
nx1 = 100 # Number of zones in X1-direction
x1min = -0.16666667 # minimum value of X1
x1max = 0.16666667 # maximum value of X1
ix1_bc = periodic # inner-X1 boundary flag
ox1_bc = periodic # outer-X1 boundary flag
nx2 = 200 # Number of zones in X2-direction
x2min = -0.5 # minimum value of X2
x2max = 0.5 # maximum value of X2
ix2_bc = reflecting # inner-X2 boundary flag
ox2_bc = reflecting # outer-X2 boundary flag
nx3 = 1 # Number of zones in X3-direction
x3min = -0.5 # minimum value of X3
x3max = 0.5 # maximum value of X3
ix3_bc = periodic # inner-X3 boundary flag
ox3_bc = periodic # outer-X3 boundary flag
<parthenon/meshblock>
nx1 = 100 # Number of cells in each MeshBlock, X1-dir
nx2 = 200 # Number of cells in each MeshBlock, X2-dir
nx3 = 1 # Number of cells in each MeshBlock, X3-dir
<parthenon/time>
integrator = vl2
cfl = 0.4
tlim = 10.0
nlim = 100000
perf_cycle_offset = 2 # number of inital cycles not to be included in perf calc
ncycle_out_mesh = -1000
<hydro>
fluid = euler
eos = adiabatic
riemann = hlle
reconstruction = ppm
gamma = 1.666666666666667 # gamma = C_p/C_v
const_accel = true # add constant accleration source term
const_accel_val = -0.1 # value of constant acceleration
const_accel_dir = 2 # direction of constant acceleration
first_order_flux_correct = true
dfloor = 0.00000001 # unused, in [code units]
pfloor = 0.00000001 # unused, in [code units]
<parthenon/output0>
file_type = hdf5
dt = 1.0
variables = prim
<parthenon/output1>
file_type = hst
dt = 0.1
Contents of the history file are changing but neither paraview, nor movie2d script is showing the evolution of fluid.
# History data
# [1]=time [2]=dt [3]=mass [4]=1-mom [5]=2-mom [6]=3-mom [7]=KE [8]=tot-E
0.00000e+00 1.03280e-03 3.33333e-09 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 5.00000e-09
1.00181e-01 1.03280e-03 3.33333e-09 -5.31666e-25 -4.29764e-24 0.00000e+00 2.81286e-39 5.00000e-09
2.00362e-01 1.03280e-03 3.33333e-09 -1.06333e-24 -8.59528e-24 0.00000e+00 1.12514e-38 5.00000e-09
3.00544e-01 1.03280e-03 3.33333e-09 -1.59500e-24 -1.28929e-23 0.00000e+00 2.53157e-38 5.00000e-09
4.00725e-01 1.03280e-03 3.33333e-09 -2.12666e-24 -1.71906e-23 0.00000e+00 4.50058e-38 5.00000e-09
5.00906e-01 1.03280e-03 3.33333e-09 -2.65833e-24 -2.14882e-23 0.00000e+00 7.03215e-38 5.00000e-09
6.00054e-01 1.03280e-03 3.33333e-09 -3.18451e-24 -2.57415e-23 0.00000e+00 1.00915e-37 5.00000e-09
7.00235e-01 1.03280e-03 3.33333e-09 -3.71618e-24 -3.00392e-23 0.00000e+00 1.37424e-37 5.00000e-09
8.00417e-01 1.03280e-03 3.33333e-09 -4.24785e-24 -3.43368e-23 0.00000e+00 1.79559e-37 5.00000e-09
9.00598e-01 1.03280e-03 3.33333e-09 -4.77951e-24 -3.86345e-23 0.00000e+00 2.27320e-37 5.00000e-09
1.00078e+00 1.03280e-03 3.33333e-09 -5.31118e-24 -4.29321e-23 0.00000e+00 2.80706e-37 5.00000e-09
1.10096e+00 1.03280e-03 3.33333e-09 -5.84284e-24 -4.72298e-23 0.00000e+00 3.39718e-37 5.00000e-09
1.20011e+00 1.03280e-03 3.33333e-09 -6.36903e-24 -5.14831e-23 0.00000e+00 4.03661e-37 5.00000e-09
1.30029e+00 1.03280e-03 3.33333e-09 -6.90070e-24 -5.57807e-23 0.00000e+00 4.73867e-37 5.00000e-09
1.40047e+00 1.03280e-03 3.33333e-09 -7.43236e-24 -6.00784e-23 0.00000e+00 5.49698e-37 5.00000e-09
1.50065e+00 1.03280e-03 3.33333e-09 -7.96403e-24 -6.43760e-23 0.00000e+00 6.31155e-37 5.00000e-09
1.60083e+00 1.03280e-03 3.33333e-09 -8.49569e-24 -6.86737e-23 0.00000e+00 7.18237e-37 5.00000e-09
1.70101e+00 1.03280e-03 3.33333e-09 -9.02736e-24 -7.29713e-23 0.00000e+00 8.10946e-37 5.00000e-09
1.80016e+00 1.03280e-03 3.33333e-09 -9.55354e-24 -7.72246e-23 0.00000e+00 9.08237e-37 5.00000e-09
1.90034e+00 1.03280e-03 3.33333e-09 -1.00852e-23 -8.15223e-23 0.00000e+00 1.01214e-36 5.00000e-09
... stripped ...
It created 10 image files, and all of them has same plot
