diff --git a/inputs/lw_implode.in b/inputs/lw_implode.in new file mode 100644 index 00000000..3a8718f1 --- /dev/null +++ b/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) + +problem = Liska Wendroff implosion + + +problem_id = 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 + + +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 + + +nx1 = 100 +nx2 = 100 +nx3 = 1 + + +file_type = hdf5 +dt = 0.1 +variables = prim +id = prim + + +file_type = hst +dt = 0.1 + + +integrator = vl2 +cfl = 0.4 +tlim = 2.5 +nlim = -1 + + +fluid = euler +eos = adiabatic +riemann = hllc +reconstruction = plm +gamma = 1.4 diff --git a/src/main.cpp b/src/main.cpp index 93f984bb..a4e1c459 100644 --- a/src/main.cpp +++ b/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; diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 382b5b7b..ec398e78 100644 --- a/src/pgen/CMakeLists.txt +++ b/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 diff --git a/src/pgen/lw_implode.cpp b/src/pgen/lw_implode.cpp new file mode 100644 index 00000000..8c36814f --- /dev/null +++ b/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 +#include + +// 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::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 diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index b3e55aad..aef501a8 100644 --- a/src/pgen/pgen.hpp +++ b/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;