Add LW implosion problem

inputs/lw_implode.in

@@ -0,0 +1,70 @@
+# AthenaPK - a performance portable block structured AMR MHD code
+# Copyright (c) 2024, Athena Parthenon Collaboration. All rights reserved.
+# Licensed under the BSD 3-Clause License (the "LICENSE");
+
+# Problem generator for square implosion problem
+# REFERENCE: R. Liska & B. Wendroff, SIAM J. Sci. Comput., 25, 995 (2003)
+<comment>
+problem = Liska Wendroff implosion
+
+<job>
+problem_id = lw_implode
+
+<problem/lw_implode>
+# Interior Conditions
+d_in  = 0.125            # density
+p_in  = 0.14             # pressure
+
+# Exterior Conditions
+d_out  = 1.0             # density
+p_out  = 1.0             # pressure
+
+<parthenon/mesh>
+refinement = none
+nghost = 2
+
+nx1 = 200
+x1min = 0.0
+x1max = 0.3
+ix1_bc = reflecting
+ox1_bc = reflecting
+
+nx2 = 200
+x2min = 0.0
+x2max = 0.3
+ix2_bc = reflecting
+ox2_bc = reflecting
+
+nx3 = 1
+x3min = -0.5
+x3max = 0.5
+ix2_bc = periodic
+ox2_bc = periodic
+
+<parthenon/meshblock>
+nx1 = 100
+nx2 = 100
+nx3 = 1
+
+<parthenon/output0>
+file_type = hdf5
+dt = 0.1
+variables = prim
+id = prim
+
+<parthenon/output1>
+file_type = hst
+dt = 0.1
+
+<parthenon/time>
+integrator = vl2
+cfl = 0.4
+tlim = 2.5
+nlim = -1
+
+<hydro>
+fluid = euler
+eos = adiabatic
+riemann = hllc
+reconstruction = plm
+gamma = 1.4

src/main.cpp

@@ -81,6 +81,8 @@ int main(int argc, char *argv[]) {
     Hydro::ProblemInitPackageData = field_loop::ProblemInitPackageData;
   } else if (problem == "kh") {
     pman.app_input->ProblemGenerator = kh::ProblemGenerator;
+  } else if (problem == "lw_implode") {
+    pman.app_input->ProblemGenerator = lw_implode::ProblemGenerator;
   } else if (problem == "rand_blast") {
     pman.app_input->ProblemGenerator = rand_blast::ProblemGenerator;
     Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData;

src/pgen/CMakeLists.txt

@@ -20,6 +20,7 @@ target_sources(athenaPK PRIVATE
     kh.cpp
     linear_wave.cpp
     linear_wave_mhd.cpp
+    lw_implode.cpp
     orszag_tang.cpp
     rand_blast.cpp
     sod.cpp

src/pgen/lw_implode.cpp

@@ -0,0 +1,77 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR MHD code
+// Copyright (c) 2024, Athena Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-Clause License (the "LICENSE")
+//========================================================================================
+//! \file lw_implode.cpp
+//! \brief Problem generator for square implosion problem
+//!
+//! REFERENCE: R. Liska & B. Wendroff, SIAM J. Sci. Comput., 25, 995 (2003)
+//========================================================================================
+
+// Parthenon headers
+#include "mesh/mesh.hpp"
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+
+// Athena headers
+#include "../main.hpp"
+#include "utils/error_checking.hpp"
+
+namespace lw_implode {
+using namespace parthenon::driver::prelude;
+
+void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) {
+  auto hydro_pkg = pmb->packages.Get("Hydro");
+  auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  PARTHENON_REQUIRE_THROWS(
+      hydro_pkg->Param<Fluid>("fluid") == Fluid::euler,
+      "Only hydro runs are supported for LW implosion problem generator.");
+
+  auto d_in = pin->GetReal("problem/lw_implode", "d_in");
+  auto p_in = pin->GetReal("problem/lw_implode", "p_in");
+
+  auto d_out = pin->GetReal("problem/lw_implode", "d_out");
+  auto p_out = pin->GetReal("problem/lw_implode", "p_out");
+
+  auto gm1 = pin->GetReal("hydro", "gamma") - 1.0;
+
+  // initialize conserved variables
+  auto &mbd = pmb->meshblock_data.Get();
+  auto &cons = mbd->Get("cons").data;
+  auto &coords = pmb->coords;
+
+  // Following Athena++
+  // https://github.com/PrincetonUniversity/athena/blob/1591aab84ba7055e5b356a8f069695ea451af8a0/src/pgen/lw_implode.cpp#L43
+  // to make sure the ICs are symmetric, set y0 to be in between cell centers
+  const auto x2min = pin->GetReal("parthenon/mesh", "x2min");
+  const auto x2max = pin->GetReal("parthenon/mesh", "x2max");
+  Real y0 = 0.5 * (x2max + x2min);
+  for (int j = jb.s; j <= jb.e; j++) {
+    if (coords.Xc<2>(j) > y0) {
+      // TODO(felker): check this condition for multi-meshblock setups
+      // further adjust y0 to be between cell center and lower x2 face
+      y0 = coords.Xf<2>(j) + 0.5 * coords.Dxf<2>(j);
+      break;
+    }
+  }
+
+  pmb->par_for(
+      "Init lw_implode", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int k, const int j, const int i) {
+        cons(IM1, k, j, i) = 0.0;
+        cons(IM2, k, j, i) = 0.0;
+        cons(IM3, k, j, i) = 0.0;
+        if (coords.Xc<2>(j) > (y0 - coords.Xc<1>(i))) {
+          cons(IDN, k, j, i) = d_out;
+          cons(IEN, k, j, i) = p_out / gm1;
+        } else {
+          cons(IDN, k, j, i) = d_in;
+          cons(IEN, k, j, i) = p_in / gm1;
+        }
+      });
+}
+} // namespace lw_implode

src/pgen/pgen.hpp

@@ -85,6 +85,13 @@ using namespace parthenon::driver::prelude;
 
 void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
 } // namespace kh
+
+namespace lw_implode {
+using namespace parthenon::driver::prelude;
+
+void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
+} // namespace lw_implode
+
 namespace rand_blast {
 using namespace parthenon::driver::prelude;