Boundary bugfix(es) and exposing FluxDiv_ interface (#558)

* Expose FluxDiv_ call to that it can be used downstream

* Fix index bug in boundaries

* change kbounds for boundaries to be the same as the others. This potentially repeats work, but is more robust and allows the user to set only the interior in the problem generator

* Rename now public FluxDiv_ function to FluxDivHelper

* Adjust unit test to fixed k bounds

* make one generic boundary condition function

* update copyright

* Update src/bvals/boundary_conditions.cpp

Co-authored-by: Philipp Grete <gretephi@msu.edu>

* make inner/outer, outflow/reflect enums instead of bools

* Add block profile to advection example

* WIP Adding non-const v to advection example. Need dt, and v fluxes are off

* Add 3D fluxes for cell advection

* Use PackIndexMap to determien variable v

* Add bval convergence test

* Add check for variables not initialized

* Fix missing header

* Address review comments

Co-authored-by: Jonah Miller <jonahm@lanl.gov>
Co-authored-by: Jonah Miller <jonah.maxwell.miller@gmail.com>
Co-authored-by: Jonas Lippuner <jlippuner@lanl.gov>

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 558]](https://github.com/lanl/parthenon/pull/558) Boundary bugfix(es) incl. regression tests and exposing FluxDiv_ interface
 - [[PR 563]](https://github.com/lanl/parthenon/pull/563) Physical boundary options for particles
 - [[PR 582]](https://github.com/lanl/parthenon/pull/582) Adding global reductions and basic functionality needed for solvers.
 - [[PR 556]](https://github.com/lanl/parthenon/pull/556) Introduce iterative tasks and regionally dependent tasks

example/advection/advection_package.cpp

@@ -25,6 +25,7 @@
 #include "defs.hpp"
 #include "kokkos_abstraction.hpp"
 #include "reconstruct/dc_inline.hpp"
+#include "utils/error_checking.hpp"
 
 using namespace parthenon::package::prelude;
 
@@ -43,6 +44,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   Real cfl = pin->GetOrAddReal("Advection", "cfl", 0.45);
   pkg->AddParam<>("cfl", cfl);
+
+  // Use constant, uniform velocity or vector valued velocity.
+  // Latter is used for testing boundary conditions.
+  auto v_const = pin->GetOrAddBoolean("Advection", "v_const", true);
   Real vx = pin->GetOrAddReal("Advection", "vx", 1.0);
   Real vy = pin->GetOrAddReal("Advection", "vy", 1.0);
   Real vz = pin->GetOrAddReal("Advection", "vz", 1.0);
@@ -53,7 +58,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   auto profile_str = pin->GetOrAddString("Advection", "profile", "wave");
   if (!((profile_str == "wave") || (profile_str == "smooth_gaussian") ||
-        (profile_str == "hard_sphere"))) {
+        (profile_str == "hard_sphere") || (profile_str == "block"))) {
     PARTHENON_FAIL(("Unknown profile in advection example: " + profile_str).c_str());
   }
   pkg->AddParam<>("profile", profile_str);
@@ -171,6 +176,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
                  std::vector<int>({vec_size}), advected_labels);
     pkg->AddField(field_name, m);
   }
+  if (!v_const) {
+    m = Metadata({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
+                  Metadata::FillGhost, Metadata::Vector},
+                 std::vector<int>({3}), std::vector<std::string>{"vx", "vy", "vz"});
+    pkg->AddField(std::string("v"), m);
+  }
   if (fill_derived) {
     field_name = "one_minus_advected";
     m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
@@ -304,6 +315,25 @@ void SquareIt(MeshBlockData<Real> *rc) {
       KOKKOS_LAMBDA(const int n, const int k, const int j, const int i) {
         v(out + n, k, j, i) = v(in + n, k, j, i) * v(in + n, k, j, i);
       });
+
+  // The following block/logic is also just added for regression testing.
+  // More specifically, the "smooth_gaussian" profile is initially != 0 everywhere, but
+  // initialializes IndexDomain::interior.
+  // FillDerived (here, SquareIt) is called after the ghost cells are exchanged and over
+  // IndexDomain::entire.
+  // Thus, no 0 (from initializing Kokkos views) should be left if all faces/corners/edges
+  // are correct, which is what we check in the loop below.
+  auto pkg = pmb->packages.Get("advection_package");
+  const auto &profile = pkg->Param<std::string>("profile");
+  if (profile == "smooth_gaussian") {
+    const auto &advected = rc->Get("advected").data;
+    pmb->par_for(
+        "advection_package::SquareIt bval check", 0, num_vars - 1, kb.s, kb.e, jb.s, jb.e,
+        ib.s, ib.e, KOKKOS_LAMBDA(const int n, const int k, const int j, const int i) {
+          PARTHENON_REQUIRE(advected(n, k, j, i) != 0.0,
+                            "Advected not properly initialized.");
+        });
+  }
 }
 
 // demonstrate usage of a "post" fill derived routine
@@ -431,7 +461,14 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
   const auto &vy = pkg->Param<Real>("vy");
   const auto &vz = pkg->Param<Real>("vz");
 
-  auto v = rc->PackVariablesAndFluxes(std::vector<MetadataFlag>{Metadata::WithFluxes});
+  PackIndexMap index_map;
+  auto v = rc->PackVariablesAndFluxes(std::vector<MetadataFlag>{Metadata::WithFluxes},
+                                      index_map);
+
+  // For non constant velocity, we need the index of the velocity vector as it's part of
+  // the variable pack.
+  const auto idx_v = index_map["v"].first;
+  const auto v_const = idx_v < 0; // using "at own perill" magic number
 
   const int scratch_level = 1; // 0 is actual scratch (tiny); 1 is HBM
   const int nx1 = pmb->cellbounds.ncellsi(IndexDomain::entire);
@@ -449,15 +486,23 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
         member.team_barrier();
 
         for (int n = 0; n < nvar; n++) {
-          if (vx > 0.0) {
-            par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
-              v.flux(X1DIR, n, k, j, i) = ql(n, i) * vx;
-            });
-          } else {
-            par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
-              v.flux(X1DIR, n, k, j, i) = qr(n, i) * vx;
-            });
-          }
+          par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+            // standard avection with fixed, global vx
+            if (v_const) {
+              if (vx > 0.0) {
+                v.flux(X1DIR, n, k, j, i) = ql(n, i) * vx;
+              } else {
+                v.flux(X1DIR, n, k, j, i) = qr(n, i) * vx;
+              }
+              // Custom flux function to move isolated, cells around. Just used for
+              // bvals testing.
+            } else {
+              v.flux(X1DIR, n, k, j, i) =
+                  ql(idx_v, i) > 0.0 ? ql(n, i) * ql(idx_v, i) : 0.0;
+              v.flux(X1DIR, n, k, j, i) +=
+                  qr(idx_v, i) < 0.0 ? qr(n, i) * qr(idx_v, i) : 0.0;
+            }
+          });
         }
       });
 
@@ -481,15 +526,25 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
           // Sync all threads in the team so that scratch memory is consistent
           member.team_barrier();
           for (int n = 0; n < nvar; n++) {
-            if (vy > 0.0) {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X2DIR, n, k, j, i) = ql(n, i) * vy;
-              });
-            } else {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X2DIR, n, k, j, i) = qr(n, i) * vy;
-              });
-            }
+            par_for_inner(member, ib.s, ib.e, [&](const int i) {
+              // standard avection with fixed, global vy
+              if (v_const) {
+                if (vy > 0.0) {
+                  v.flux(X2DIR, n, k, j, i) = ql(n, i) * vy;
+                } else {
+                  v.flux(X2DIR, n, k, j, i) = qr(n, i) * vy;
+                }
+                // Custom flux function to move isolated, cells around. Just used for
+                // bvals testing.
+              } else {
+                v.flux(X2DIR, n, k, j, i) = ql(idx_v + X2DIR - 1, i) > 0.0
+                                                ? ql(n, i) * ql(idx_v + X2DIR - 1, i)
+                                                : 0.0;
+                v.flux(X2DIR, n, k, j, i) += qr(idx_v + X2DIR - 1, i) < 0.0
+                                                 ? qr(n, i) * qr(idx_v + X2DIR - 1, i)
+                                                 : 0.0;
+              }
+            });
           }
         });
   }
@@ -514,15 +569,25 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
           // Sync all threads in the team so that scratch memory is consistent
           member.team_barrier();
           for (int n = 0; n < nvar; n++) {
-            if (vz > 0.0) {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X3DIR, n, k, j, i) = ql(n, i) * vz;
-              });
-            } else {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X3DIR, n, k, j, i) = qr(n, i) * vz;
-              });
-            }
+            par_for_inner(member, ib.s, ib.e, [&](const int i) {
+              // standard avection with fixed, global vz
+              if (v_const) {
+                if (vz > 0.0) {
+                  v.flux(X3DIR, n, k, j, i) = ql(n, i) * vz;
+                } else {
+                  v.flux(X3DIR, n, k, j, i) = qr(n, i) * vz;
+                }
+                // Custom flux function to move isolated, cells around. Just used for
+                // bvals testing.
+              } else {
+                v.flux(X3DIR, n, k, j, i) = ql(idx_v + X3DIR - 1, i) > 0.0
+                                                ? ql(n, i) * ql(idx_v + X3DIR - 1, i)
+                                                : 0.0;
+                v.flux(X3DIR, n, k, j, i) += qr(idx_v + X3DIR - 1, i) < 0.0
+                                                 ? qr(n, i) * qr(idx_v + X3DIR - 1, i)
+                                                 : 0.0;
+              }
+            });
           }
         });
   }

example/advection/parthenon_app_inputs.cpp

@@ -20,6 +20,7 @@
 #include "advection_package.hpp"
 #include "config.hpp"
 #include "defs.hpp"
+#include "interface/variable_pack.hpp"
 #include "utils/error_checking.hpp"
 
 using namespace parthenon::package::prelude;
@@ -39,6 +40,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto pkg = pmb->packages.Get("advection_package");
   const auto &amp = pkg->Param<Real>("amp");
   const auto &vel = pkg->Param<Real>("vel");
+  const auto &vx = pkg->Param<Real>("vx");
+  const auto &vy = pkg->Param<Real>("vy");
+  const auto &vz = pkg->Param<Real>("vz");
   const auto &k_par = pkg->Param<Real>("k_par");
   const auto &cos_a2 = pkg->Param<Real>("cos_a2");
   const auto &cos_a3 = pkg->Param<Real>("cos_a3");
@@ -47,18 +51,26 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const auto &profile = pkg->Param<std::string>("profile");
 
   auto cellbounds = pmb->cellbounds;
-  IndexRange ib = cellbounds.GetBoundsI(IndexDomain::entire);
-  IndexRange jb = cellbounds.GetBoundsJ(IndexDomain::entire);
-  IndexRange kb = cellbounds.GetBoundsK(IndexDomain::entire);
+  IndexRange ib = cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = cellbounds.GetBoundsK(IndexDomain::interior);
 
   auto coords = pmb->coords;
-  auto q = data->PackVariables(std::vector<MetadataFlag>{Metadata::Independent});
+  PackIndexMap index_map;
+  auto q =
+      data->PackVariables(std::vector<MetadataFlag>{Metadata::Independent}, index_map);
   const auto num_vars = q.GetDim(4);
 
+  // For non constant velocity, we need the index of the velocity vector as it's part of
+  // the variable pack.
+  const auto idx_v = index_map["v"].first;
+  const auto idx_adv = index_map.get("advected").first;
+
   int profile_type;
   if (profile == "wave") profile_type = 0;
   if (profile == "smooth_gaussian") profile_type = 1;
   if (profile == "hard_sphere") profile_type = 2;
+  if (profile == "block") profile_type = 3;
 
   pmb->par_for(
       "Advection::ProblemGenerator", 0, num_vars - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
@@ -80,6 +92,29 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           q(n, k, j, i) = 0.0;
         }
       });
+
+  const auto block_id = pmb->gid;
+  // initialize some arbitrary cells in the first block that move in all 6 directions
+  if (profile_type == 3 && block_id == 0) {
+    pmb->par_for(
+        "Advection::ProblemGenerator bvals test", 0, 1,
+        KOKKOS_LAMBDA(const int /*unused*/) {
+          q(idx_adv, 4, 4, 4) = 10.0;
+          q(idx_v, 4, 4, 4) = vx;
+          q(idx_adv, 4, 6, 4) = 10.0;
+          q(idx_v, 4, 6, 4) = -vx;
+
+          q(idx_adv, 6, 4, 4) = 10.0;
+          q(idx_v + 1, 6, 4, 4) = vy;
+          q(idx_adv, 6, 6, 6) = 10.0;
+          q(idx_v + 1, 6, 6, 6) = -vy;
+
+          q(idx_adv, 4, 8, 8) = 10.0;
+          q(idx_v + 2, 4, 8, 8) = vz;
+          q(idx_adv, 4, 10, 8) = 10.0;
+          q(idx_v + 2, 4, 10, 8) = -vz;
+        });
+  }
 }
 
 //========================================================================================

scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py

@@ -520,7 +520,9 @@ def compare(
                 bad_idx = tuple(bad_idx)
 
                 # Find the bad location
-                bad_loc = loc[:, bad_idx[0], bad_idx[1], bad_idx[2], bad_idx[3]]
+                bad_loc = np.array(loc)[
+                    :, bad_idx[0], bad_idx[1], bad_idx[2], bad_idx[3]
+                ]
 
                 # TODO(forrestglines): Check that the bkji and zyx reported are the correct order
                 print(f"Diff in {var:20s}")