Implement RTI problem

inputs/rt.in

@@ -0,0 +1,69 @@
+<comment>
+problem = Rayleigh-Taylor instability
+
+<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    = outflow     # inner-X2 boundary flag
+ox2_bc    = outflow     # 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

src/main.cpp

@@ -84,6 +84,8 @@ int main(int argc, char *argv[]) {
     pman.app_input->ProblemGenerator = rand_blast::ProblemGenerator;
     Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData;
     Hydro::ProblemSourceFirstOrder = rand_blast::RandomBlasts;
+  } else if (problem == "rt") {
+    pman.app_input->ProblemGenerator = rt::ProblemGenerator;
   } else if (problem == "cluster") {
     pman.app_input->ProblemGenerator = cluster::ProblemGenerator;
     pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput;

src/pgen/CMakeLists.txt

@@ -19,6 +19,7 @@ target_sources(athenaPK PRIVATE
     linear_wave_mhd.cpp
     orszag_tang.cpp
     rand_blast.cpp
+    rt.cpp
     sod.cpp
     turbulence.cpp
     )

src/pgen/pgen.hpp

@@ -121,4 +121,9 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin);
 void Cleanup();
 } // namespace turbulence
 
+namespace rt {
+  using namespace parthenon::driver::prelude;
+  void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
+}
+
 #endif // PGEN_PGEN_HPP_

src/pgen/rt.cpp

@@ -0,0 +1,111 @@
+//! \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;
+            }
+        }
+      }
+    }
+    u_dev.DeepCopy(u);
+  }
+
+  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