Add reflecting boundary conditions (#140)

* bump Parthenon

* Bump Parth to include restart ghost fix

* update parthenon to openpmd branch

* Bump Parthenon and Kokkos

* Bump Kokkos 4.4.0 and Parthenon 24.08

* Add LW implosion problem

* Fix history filename

* Update interface (cleanup and workaround)

* cpp-py-formatter

* Fix typo

* Remove unused vars

* Add hydro bvals and LW implosion pgen

* Add more details on outflow bound. conditions

* Add reflection symmetry regression test

---------

Co-authored-by: Ben Wibking <ben@wibking.com>
Co-authored-by: par-hermes <par-hermes@automation.bot>

Author: 3 people

CHANGELOG.md

@@ -10,6 +10,7 @@ In between a subtle bug was introduced that resulted in inconsistent divergence
 refinement.
 
 ### Added (new features/APIs/variables/...)
+- [[PR 140]](https://github.com/parthenon-hpc-lab/athenapk/pull/140) Add hydro reflecting boundary conditions
 - [[PR102]](https://github.com/parthenon-hpc-lab/athenapk/pull/102) Add support for tracer particles
 - [[PR 89]](https://github.com/parthenon-hpc-lab/athenapk/pull/89) Add viscosity and resistivity
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping

docs/input.md

@@ -407,4 +407,24 @@ Following "restrictions" apply to the current tracer implementation:
 - Default tracer values (such as the primitive fields) are only updated right before writing an output file.
 - Ids are only unique if tracers are seeded at the beginning at the simulations and no new tracer particles are added dynamically while the simulation is running.
 
-Please get in touch, if you interested in running simulations that require lifting one (or more) of those restrictions.
+Please get in touch, if you interested in running simulations that require lifting one (or more) of those restrictions.
+
+## Boundary conditions
+
+In addition to enrolling custom boundary conditions, three general options are currently supported by default:
+
+- `periodic`: Work without restrictions (implemented in upstream Parthenon)
+- `outflow`: Last layer of active cells is copied into ghost cells. Work (technicaly) without restrictions (implemented in upstream Parthenon). Physically care must be taken when using outflow conditions in subsonic flows (as characteristics can travel into the domain, which is not properly handled by copying the active cell data).
+It is recommended to implement custom, problem specific boundary conditions for subsonic outflows, see large body of engineering literature.
+For supersonic outflows there is no issue as the flow itself is faster than the fastest characteristic.
+- `reflecting`: Work only for cell-centered *hydro* fluid variables
+(implemented in AthenaPK currently without support for particles or `Metadata::Fine` fields,
+ where the latter are currently not being used in AthenaPK).
+ Note, that we specifically do *not* provide "reflecting" boundary conditions for the magnetic
+ fields as these kind of configurations are typically highly problem dependent and, thus,
+ require manual treatment.
+
+They are being used by specifying the name as parameter value for the
+`ix1_bc`, `ox1_bc`, `ix2_bc`, `ox2_bc`, `ix3_bc`, and `ox3_bc` paraemters
+in the `<parthenon/mesh>` block of the input file for the inner and outer
+(say left and right in x1 direction) in each coordinate direction, respectively.

inputs/lw_implode.in

@@ -0,0 +1,71 @@
+# AthenaPK - a performance portable block structured AMR MHD code
+# Copyright (c) 2024-2025, 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 = 3
+
+nx1 = 256
+x1min = 0.0
+x1max = 0.3
+ix1_bc = reflecting
+ox1_bc = reflecting
+
+nx2 = 256
+x2min = 0.0
+x2max = 0.3
+ix2_bc = reflecting
+ox2_bc = reflecting
+
+nx3 = 1
+x3min = -0.5
+x3max = 0.5
+ix3_bc = periodic
+ox3_bc = periodic
+
+<parthenon/meshblock>
+nx1 = 256
+nx2 = 256
+nx3 = 1
+
+<parthenon/output0>
+file_type = hdf5
+dt = 0.1
+#dn = 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/CMakeLists.txt

@@ -6,6 +6,7 @@ add_executable(
         main.cpp
         eos/adiabatic_glmmhd.cpp
         units.hpp
+        bvals/boundary_conditions_apk.hpp
         eos/adiabatic_hydro.cpp
         hydro/diffusion/conduction.cpp
         hydro/diffusion/diffusion.cpp

src/bvals/boundary_conditions_apk.hpp

@@ -0,0 +1,90 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2025, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+//! \file boundary_conditions_apk.chpp
+//  \brief AthenaPK specific boundary conditions
+//
+
+#ifndef BVALS_BOUNDARY_CONDITIONS_APK_HPP_
+#define BVALS_BOUNDARY_CONDITIONS_APK_HPP_
+
+#include <memory>
+#include <string>
+#include <vector>
+
+// Parthenon headers
+#include <parthenon/package.hpp>
+
+#include "basic_types.hpp"
+#include "bvals/boundary_conditions_generic.hpp"
+#include "mesh/domain.hpp"
+#include "mesh/mesh.hpp"
+#include "mesh/meshblock.hpp"
+#include "utils/error_checking.hpp"
+
+#include "../main.hpp"
+
+namespace Hydro {
+namespace BoundaryFunction {
+
+using namespace parthenon::package::prelude;
+using parthenon::CoordinateDirection;
+// using parthenon::MeshBlockData;
+// using parthenon::Real;
+using parthenon::BoundaryFunction::BCSide;
+
+template <CoordinateDirection DIR, BCSide SIDE>
+void ReflectBC(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse) {
+  // make sure DIR is X[123]DIR so we don't have to check again
+  static_assert(DIR == X1DIR || DIR == X2DIR || DIR == X3DIR, "DIR must be X[123]DIR");
+
+  MeshBlock *pmb = mbd->GetBlockPointer();
+
+  auto hydro_pkg = pmb->packages.Get("Hydro");
+  auto fluid = hydro_pkg->Param<Fluid>("fluid");
+  PARTHENON_REQUIRE_THROWS(
+      fluid == Fluid::euler,
+      "Reflecting boundary conditions for MHD need special treatment.");
+
+  // convenient shorthands
+  constexpr bool X1 = (DIR == X1DIR);
+  constexpr bool X2 = (DIR == X2DIR);
+  constexpr bool X3 = (DIR == X3DIR);
+  constexpr bool INNER = (SIDE == BCSide::Inner);
+
+  const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
+
+  const auto &range = X1 ? bounds.GetBoundsI(IndexDomain::interior)
+                         : (X2 ? bounds.GetBoundsJ(IndexDomain::interior)
+                               : bounds.GetBoundsK(IndexDomain::interior));
+  const int ref = INNER ? range.s : range.e;
+
+  constexpr IndexDomain domain =
+      INNER ? (X1 ? IndexDomain::inner_x1
+                  : (X2 ? IndexDomain::inner_x2 : IndexDomain::inner_x3))
+            : (X1 ? IndexDomain::outer_x1
+                  : (X2 ? IndexDomain::outer_x2 : IndexDomain::outer_x3));
+
+  // used for reflections
+  const int offset = (2 * ref) + (INNER ? -1 : 1);
+
+  auto cons = mbd->PackVariables(std::vector<std::string>{"cons"}, coarse);
+  const bool fine = false; // no usage of fine fields in AthenaPK for now
+
+  const auto nv = IndexRange{0, cons.GetDim(4) - 1};
+  pmb->par_for_bndry(
+      "ReflectBC", nv, domain, parthenon::TopologicalElement::CC, coarse, fine,
+      KOKKOS_LAMBDA(const int &v, const int &k, const int &j, const int &i) {
+        const bool reflect = v == DIR;
+        cons(v, k, j, i) =
+            (reflect ? -1.0 : 1.0) *
+            cons(v, X3 ? offset - k : k, X2 ? offset - j : j, X1 ? offset - i : i);
+      });
+}
+
+} // namespace BoundaryFunction
+} // namespace Hydro
+
+#endif // BVALS_BOUNDARY_CONDITIONS_APK_HPP_

src/hydro/hydro_driver.cpp

@@ -111,7 +111,6 @@ TaskStatus RKL2StepFirst(MeshData<Real> *md_Y0, MeshData<Real> *md_Yjm1,
   auto Yjm2 = md_Yjm2->PackVariablesAndFluxes(flags_ind);
   auto MY0 = md_MY0->PackVariablesAndFluxes(flags_ind);
 
-  const int ndim = pmb->pmy_mesh->ndim;
   // Using separate loops for each dim as the launch overhead should be hidden
   // by enough work over the entire pack and it allows to not use any conditionals.
   parthenon::par_for(

src/main.cpp

@@ -5,10 +5,13 @@
 #include <sstream>
 
 // Parthenon headers
+#include "bvals/boundary_conditions_generic.hpp"
+#include "defs.hpp"
 #include "globals.hpp"
 #include "parthenon_manager.hpp"
 
 // AthenaPK headers
+#include "bvals/boundary_conditions_apk.hpp"
 #include "hydro/hydro.hpp"
 #include "hydro/hydro_driver.hpp"
 #include "main.hpp"
@@ -89,6 +92,8 @@ int main(int argc, char *argv[]) {
     Hydro::ProblemInitPackageData = field_loop::ProblemInitPackageData;
   } else if (problem == "kh") {
     pman.app_input->MeshProblemGenerator = 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;
@@ -117,6 +122,23 @@ int main(int argc, char *argv[]) {
     PARTHENON_THROW(msg);
   }
 
+  const std::string REFLECTING = "reflecting";
+  using BF = parthenon::BoundaryFace;
+  using Hydro::BoundaryFunction::ReflectBC;
+  using parthenon::BoundaryFunction::BCSide;
+  pman.app_input->RegisterBoundaryCondition(BF::inner_x1, REFLECTING,
+                                            ReflectBC<X1DIR, BCSide::Inner>);
+  pman.app_input->RegisterBoundaryCondition(BF::outer_x1, REFLECTING,
+                                            ReflectBC<X1DIR, BCSide::Outer>);
+  pman.app_input->RegisterBoundaryCondition(BF::inner_x2, REFLECTING,
+                                            ReflectBC<X2DIR, BCSide::Inner>);
+  pman.app_input->RegisterBoundaryCondition(BF::outer_x2, REFLECTING,
+                                            ReflectBC<X2DIR, BCSide::Outer>);
+  pman.app_input->RegisterBoundaryCondition(BF::inner_x3, REFLECTING,
+                                            ReflectBC<X3DIR, BCSide::Inner>);
+  pman.app_input->RegisterBoundaryCondition(BF::outer_x3, REFLECTING,
+                                            ReflectBC<X3DIR, BCSide::Outer>);
+
   pman.ParthenonInitPackagesAndMesh();
 
   // Startup the corresponding driver for the integrator

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/cluster/agn_feedback.cpp

@@ -290,7 +290,6 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
 
   // Velocity of added gas
   const Real jet_velocity = kinetic_jet_velocity_;
-  const Real jet_specific_internal_e = kinetic_jet_e_;
 
   // Amount of momentum density ( density * velocity) to dump in each cell
   const Real jet_momentum = jet_density * jet_velocity;

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