Doubt in problem output (Not generating the output)

[tbhaxor]

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

[pgrete]

Hi @tbhaxor
it looks like you're initializing the data on the hos

    // initializing on host
    auto u = u_dev.GetHostMirrorAndCopy();

but do not copy the data back to the device.
Either you update the initialization to be done on the device (i.e., using a par_for) or you add a deep copy after the initialization loop to copy the data from the host to the device.
I recommend the former as the "init on host" is a leftover from refactoring Athena++ code and we should clean that up at some point.

[tbhaxor]

Hi @pgrete

Could you guide me into direction with some snippets?

I recommend the former..

So I am using auto u = u_dev.GetDeviceMirror();, but still same issue

[pgrete]

Here's a reference to "init on host" with the deep copy back to the device at the end

athenapk/src/pgen/blast.cpp

Line 202 in 8f04e18

u_dev.DeepCopy(u);

and here's one with the initialization on the device

athenapk/src/pgen/orszag_tang.cpp

Line 40 in 8f04e18

pmb->par_for(

[tbhaxor]

I tried both of them, neither of them are working in my case. Could it be some problem with problem setup?

[pgrete]

Could you please open a PR with the exact (full) code that you tried, i.e., including the modification to main.cpp etc? That' might help in narrowing down the issue.

[tbhaxor]

Hi @pgrete, thanks for waiting. Here is the PR #70