Merge branch 'develop' into lroberts36/add-sparse-vector-wave-test

.github/workflows/ci-macos.yml

@@ -24,23 +24,27 @@ jobs:
     - uses: actions/setup-python@v5
       with:
         python-version: '3.11'
-        cache: 'pip' # caching pip dependencies
+        cache: 'pip'
+        cache-dependency-path: '**/requirements.txt'
+    - run: pip install -r requirements.txt
+
+    - name: Install dependencies
+      run: |
+        brew install openmpi hdf5-mpi adios2 || true
 
     - name: Configure
-        run: |
-          cmake -B build \
-            -DCMAKE_BUILD_TYPE=Release \
+      run: cmake -B build -DCMAKE_BUILD_TYPE=Release
 
     - name: Build
-        run: cmake --build build
+      run: cmake --build build
 
     - name: Test
-        run: |
-          cd build
-          ctest -LE 'regression'
-
-     - uses: actions/upload-artifact@v3
-        with:
-          name: configure-log-unit-macos
-          path: build/CMakeFiles/CMakeOutput.log
-          retention-days: 3
+      run: |
+        cd build
+        ctest -LE mpi
+
+    - uses: actions/upload-artifact@v3
+      with:
+        name: configure-log-unit-macos
+        path: build/CMakeFiles/CMakeOutput.log
+        retention-days: 3

CHANGELOG.md

@@ -5,7 +5,11 @@
 ### Added (new features/APIs/variables/...)
 - [[PR 1091]](https://github.com/parthenon-hpc-lab/parthenon/pull/1091) Add vector wave equation example.
 - [[PR 991]](https://github.com/parthenon-hpc-lab/parthenon/pull/991) Add fine fields.
-- [[PR 1084]](https://github.com/parthenon-hpc-lab/parthenon/pull/1084) Properly free swarm boundary MPI requests 
+- [[PR 1106]](https://github.com/parthenon-hpc-lab/parthenon/pull/1106) Add CMake options for turning on ASAN and HWASAN
+- [[PR 1100]](https://github.com/parthenon-hpc-lab/parthenon/pull/1100) Custom refinement ops propagated to fluxes
+- [[PR 1090]](https://github.com/parthenon-hpc-lab/parthenon/pull/1090) SMR with swarms
+- [[PR 1079]](https://github.com/parthenon-hpc-lab/parthenon/pull/1079) Address XDMF/Visit Issues
+- [[PR 1084]](https://github.com/parthenon-hpc-lab/parthenon/pull/1084) Properly free swarm boundary MPI requests
 - [[PR 1020]](https://github.com/parthenon-hpc-lab/parthenon/pull/1020) Add bi- and trilinear interpolation routines
 - [[PR 1081]](https://github.com/parthenon-hpc-lab/parthenon/pull/1081) Add GetSize and GetIndex to sparse pack
 - [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
@@ -23,6 +27,15 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1092]](https://github.com/parthenon-hpc-lab/parthenon/pull/1092) Updates to DataCollection and MeshData to remove requirement of predefining MeshBlockData
+- [[PR 1113]](https://github.com/parthenon-hpc-lab/parthenon/pull/1113) Prevent division by zero
+- [[PR 1112]](https://github.com/parthenon-hpc-lab/parthenon/pull/1112) Remove shared_ptr cycle in forest::Tree
+- [[PR 1104]](https://github.com/parthenon-hpc-lab/parthenon/pull/1104) Fix reading restarts due to hidden ghost var
+- [[PR 1098]](https://github.com/parthenon-hpc-lab/parthenon/pull/1098) Move to symmetrized logical coordinates and fix SMR bug
+- [[PR 1095]](https://github.com/parthenon-hpc-lab/parthenon/pull/1095) Add missing include guards in hdf5 restart
+- [[PR 1093]](https://github.com/parthenon-hpc-lab/parthenon/pull/1093) Fix forest size for symmetry dimensions
+- [[PR 1089]](https://github.com/parthenon-hpc-lab/parthenon/pull/1089) Fix loading restart files without derefinement counter
+- [[PR 1079]](https://github.com/parthenon-hpc-lab/parthenon/pull/1079) Address XDMF/Visit Issues
 - [[PR 1088]](https://github.com/parthenon-hpc-lab/parthenon/pull/1088) Correctly fill fluxes for non-cell variables in SparsePacks
 - [[PR 1083]](https://github.com/parthenon-hpc-lab/parthenon/pull/1083) Correctly fill VariableFluxPack for edge fluxes in 2D
 - [[PR 1087]](https://github.com/parthenon-hpc-lab/parthenon/pull/1087) Make sure InnerLoopPatternTVR is resolved on device properly when it is the default loop pattern
@@ -37,6 +50,7 @@
 - [[PR 1031]](https://github.com/parthenon-hpc-lab/parthenon/pull/1031) Fix bug in non-cell centered AMR
 
 ### Infrastructure (changes irrelevant to downstream codes)
+- [[PR 1055]](https://github.com/parthenon-hpc-lab/parthenon/pull/1055) Refactor mesh constructors
 - [[PR 1066]](https://github.com/parthenon-hpc-lab/parthenon/pull/1066) Re-introduce default loop patterns and exec spaces
 - [[PR 1064]](https://github.com/parthenon-hpc-lab/parthenon/pull/1064) Forbid erroneous edge case when adding MeshData on a partition
 - [[PR 1035]](https://github.com/parthenon-hpc-lab/parthenon/pull/1035) Fix multigrid infrastructure to work with forest
@@ -46,9 +60,10 @@
 
 
 ### Removed (removing behavior/API/varaibles/...)
-
+- [[PR 1108]](https://github.com/parthenon-hpc-lab/parthenon/pull/1108) Remove NaN payload tags infrastructure
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 1108]](https://github.com/parthenon-hpc-lab/parthenon/pull/1108) Remove NaN payload tags infrastructure
 - [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
 - [[PR 1037]](https://github.com/parthenon-hpc-lab/parthenon/pull/1037) Add SwarmPacks
 - [[PR 1042]](https://github.com/parthenon-hpc-lab/parthenon/pull/1042) Use Offset class and clean up of NeighborBlock

CMakeLists.txt

@@ -57,6 +57,8 @@ option(CHECK_REGISTRY_PRESSURE "Check the registry pressure for Kokkos CUDA kern
 option(TEST_INTEL_OPTIMIZATION "Test intel optimization and vectorization" OFF)
 option(TEST_ERROR_CHECKING "Enables the error checking unit test. This test will FAIL" OFF)
 option(CODE_COVERAGE "Enable code coverage reporting" OFF)
+option(ENABLE_ASAN "Turn on ASAN" OFF)
+option(ENABLE_HWASAN "Turn on HWASAN (currently ARM-only)" OFF)
 
 include(cmake/Format.cmake)
 include(cmake/Lint.cmake)
@@ -290,7 +292,19 @@ if (Kokkos_ENABLE_CUDA AND "${PARTHENON_ENABLE_GPU_MPI_CHECKS}" )
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cmake/CTestCustom.cmake.in ${CMAKE_BINARY_DIR}/CTestCustom.cmake @ONLY)
 endif()
 
-
+# option to turn on AddressSanitizer for debugging
+if(ENABLE_ASAN)
+  message(STATUS "Compiling with AddressSanitizer and UndefinedBehaviorSanitizer *enabled*")
+  add_compile_options(-fsanitize=address -fsanitize=undefined -fno-sanitize-recover=all -fsanitize=float-divide-by-zero -fsanitize=float-cast-overflow -fno-sanitize=null -fno-sanitize=alignment)
+  add_link_options(-fsanitize=address -fsanitize=undefined)
+endif(ENABLE_ASAN)
+
+# option to turn on HWAddressSanitizer for debugging
+if(ENABLE_HWASAN)
+  message(STATUS "Compiling with HWAddressSanitizer *enabled*")
+  add_compile_options(-fsanitize=hwaddress)
+  add_link_options(-fsanitize=hwaddress)
+endif(ENABLE_HWASAN)
 
 
 # Build Tests and download Catch2

doc/sphinx/src/interface/refinement_operations.rst

@@ -110,3 +110,10 @@ You must register both prolongation and restriction together. You may,
 however, use the default Parthenon structs if desired. Then any variable
 registered with this metadata object will use your custom prolongation
 and restriction operations.
+
+When a variable with custom operations is enrolled and marked
+``Metadata::WithFluxes``, the resulting flux variables that are created will
+also have the same custom operations enrolled. In general the custom operations
+will need to be different for the variable and for its fluxes; these can be
+distinguished inside the custom operations by referring to the
+``TopologicalElement`` template parameter.

example/fine_advection/CMakeLists.txt

@@ -1,5 +1,5 @@
 #=========================================================================================
-# (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+# (C) (or copyright) 2024. Triad National Security, LLC. All rights reserved.
 #
 # This program was produced under U.S. Government contract 89233218CNA000001 for Los
 # Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC

example/fine_advection/README.md

@@ -0,0 +1,5 @@
+This example implements upwind advection of a cell-centered scalar variable defined 
+on the regular grid and for another cell-centered variable on the fine grid (which
+is twice the resolution and is selected using Metadata::Fine). The newer type-based
+`SparsePack`s are used throughout and machinery for doing a generalized Stoke's 
+theorem based update is included.

example/fine_advection/advection_driver.cpp

@@ -75,17 +75,6 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in
   // Build MeshBlockData containers that will be included in MeshData containers. It is
   // gross that this has to be done by hand.
   const auto &stage_name = integrator->stage_name;
-  if (stage == 1) {
-    for (int i = 0; i < blocks.size(); i++) {
-      auto &pmb = blocks[i];
-      // first make other useful containers
-      auto &base = pmb->meshblock_data.Get();
-      pmb->meshblock_data.Add("dUdt", base);
-      for (int s = 1; s < integrator->nstages; s++)
-        pmb->meshblock_data.Add(stage_name[s], base);
-    }
-  }
-
   const Real beta = integrator->beta[stage - 1];
   const Real dt = integrator->dt;
 
@@ -102,6 +91,10 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in
     auto start_send = tl.AddTask(none, parthenon::StartReceiveBoundaryBuffers, mc1);
     auto start_flxcor = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mc0);
 
+    // Make a sparse variable pack descriptors that can be used to build packs
+    // including some subset of the fields in this example. This will be passed
+    // to the Stokes update routines, so that they can internally create variable
+    // packs that operate on only the desired set of variables.
     using namespace advection_package::Conserved;
     static auto desc = parthenon::MakePackDescriptor<phi>(
         pmesh->resolved_packages.get(), {parthenon::Metadata::WithFluxes},

example/fine_advection/advection_package.cpp

@@ -39,7 +39,6 @@ using namespace parthenon::package::prelude;
 // *************************************************//
 
 namespace advection_package {
-using parthenon::UserHistoryOperation;
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto pkg = std::make_shared<StateDescriptor>("advection_package");

example/fine_advection/parthinput.advection

@@ -1,8 +1,4 @@
 # ========================================================================================
-#  Athena++ astrophysical MHD code
-#  Copyright(C) 2014 James M. Stone <[email protected]> and other code contributors
-#  Licensed under the 3-clause BSD License, see LICENSE file for details
-# ========================================================================================
 #  (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 #
 #  This program was produced under U.S. Government contract 89233218CNA000001 for Los

example/fine_advection/stokes.hpp

@@ -43,6 +43,14 @@ TaskStatus WeightedSumDataElement(parthenon::CellLevel cl,
       PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack1.GetNBlocks() - 1, kb.s,
       kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
         parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e});
+        PARTHENON_REQUIRE(pack1.GetLowerBound(b) == pack2.GetLowerBound(b),
+                          "Packs are different size.");
+        PARTHENON_REQUIRE(pack1.GetLowerBound(b) == pack_out.GetLowerBound(b),
+                          "Packs are different size.");
+        PARTHENON_REQUIRE(pack1.GetUpperBound(b) == pack2.GetUpperBound(b),
+                          "Packs are different size.");
+        PARTHENON_REQUIRE(pack1.GetUpperBound(b) == pack_out.GetUpperBound(b),
+                          "Packs are different size.");
         for (int l = pack1.GetLowerBound(b); l <= pack1.GetUpperBound(b); ++l) {
           parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) {
             const auto [j, i] = idxer(idx);
@@ -115,6 +123,10 @@ void StokesComponent(Real fac, parthenon::CellLevel cl,
       kb.s, kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
         auto &coords = pack_in.GetCoordinates(b);
         parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e});
+        PARTHENON_REQUIRE(pack_in.GetLowerBound(b) == pack_out.GetLowerBound(b),
+                          "Packs are different size.");
+        PARTHENON_REQUIRE(pack_in.GetUpperBound(b) == pack_out.GetUpperBound(b),
+                          "Packs are different size.");
         for (int l = pack_out.GetLowerBound(b); l <= pack_out.GetUpperBound(b); ++l) {
           parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) {
             const auto [j, i] = idxer(idx);

example/particle_leapfrog/parthinput.particle_leapfrog

@@ -3,7 +3,7 @@
 #  Copyright(C) 2021 The Parthenon collaboration
 #  Licensed under the 3-clause BSD License, see LICENSE file for details
 # ========================================================================================
-#  (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved.
+#  (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
 #
 #  This program was produced under U.S. Government contract 89233218CNA000001 for Los
 #  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -19,7 +19,7 @@
 problem_id = particles
 
 <parthenon/mesh>
-refinement = none
+refinement = static
 
 nx1 = 16
 x1min = -0.5
@@ -44,6 +44,15 @@ nx1 = 8
 nx2 = 8
 nx3 = 8
 
+<parthenon/static_refinement1>
+level = 1
+x1min = -0.5
+x1max = 0
+x2min = -0.5
+x2max = 0
+x3min = -0.5
+x3max = 0
+
 <parthenon/time>
 tlim = 1.0
 nlim = 100000

example/poisson/poisson_package.cpp

@@ -22,6 +22,7 @@
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 #include <solvers/solver_utils.hpp>
+#include <utils/robust.hpp>
 
 #include "defs.hpp"
 #include "kokkos_abstraction.hpp"
@@ -296,7 +297,8 @@ TaskStatus CheckConvergence(T *u, T *du) {
       parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0,
       v.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &eps) {
-        Real reps = std::abs(dv(b, idphi, k, j, i) / v(b, iphi, k, j, i));
+        Real reps = std::abs(
+            parthenon::robust::ratio(dv(b, idphi, k, j, i), v(b, iphi, k, j, i)));
         Real aeps = std::abs(dv(b, idphi, k, j, i));
         eps = std::max(eps, std::min(reps, aeps));
       },

requirements.txt

@@ -0,0 +1,3 @@
+h5py
+numpy
+matplotlib

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

@@ -285,7 +285,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
             for key in f0.fid[var].keys():
                 if var == "Blocks" and key == "loc.level-gid-lid-cnghost-gflag":
                     continue  # depends on number of MPI ranks and distribution of blocks among ranks
-
+                if var == "Blocks" and key == "derefinement_count":
+                    continue  # Gold files have not been updated to include derefinement_count
                 # Compare raw data of these variables
                 val0 = f0.fid[var][key]
                 val1 = f1.fid[var][key]

src/CMakeLists.txt

@@ -258,7 +258,6 @@ add_library(parthenon
   utils/morton_number.hpp
   utils/mpi_types.hpp
   utils/multi_pointer.hpp
-  utils/nan_payload_tag.hpp
   utils/object_pool.hpp
   utils/partition_stl_containers.hpp
   utils/reductions.hpp

src/bvals/comms/bnd_info.cpp

@@ -267,6 +267,8 @@ BndInfo::BndInfo(MeshBlock *pmb, const NeighborBlock &nb,
     var = v->data.Get();
   }
 
+  coords = pmb->coords;
+
   auto elements = v->GetTopologicalElements();
   if (v->IsSet(Metadata::Flux)) elements = GetFluxCorrectionElements(v, nb.offsets);
   ntopological_elements = elements.size();

src/defs.hpp

@@ -94,8 +94,9 @@ struct RegionSize {
 
   // Returns global coordinate position within a block based on block local
   // coordinate u running from zero to one
-  Real LogicalToActualPosition(Real u, CoordinateDirection dir) const {
-    return u * (xmax_[dir - 1] - xmin_[dir - 1]) + xmin_[dir - 1];
+  Real SymmetrizedLogicalToActualPosition(Real u, CoordinateDirection dir) const {
+    return static_cast<Real>(0.5) * (xmin_[dir - 1] + xmax_[dir - 1]) +
+           (u * xmax_[dir - 1] - u * xmin_[dir - 1]);
   }
   // A "symmetry" direction is a a direction that posesses a translational symmetry
   // (or rotational symmetry, etc. for non-cartesian coordinate systems) in the given

src/interface/data_collection.cpp

@@ -32,30 +32,26 @@ std::shared_ptr<T> &DataCollection<T>::Add(const std::string &label) {
   return containers_[label];
 }
 
+template <>
 std::shared_ptr<MeshData<Real>> &
-GetOrAdd_impl(Mesh *pmy_mesh_,
-              std::map<std::string, std::shared_ptr<MeshData<Real>>> &containers_,
-              BlockList_t &block_list, const std::string &mbd_label,
-              const int &partition_id, const std::optional<int> gmg_level) {
-  std::string label = mbd_label + "_part-" + std::to_string(partition_id);
-  if (gmg_level) label = label + "_gmg-" + std::to_string(*gmg_level);
+DataCollection<MeshData<Real>>::GetOrAdd_impl(const std::string &mbd_label,
+                                              const int &partition_id,
+                                              const std::optional<int> gmg_level) {
+  std::string label = GetKey(mbd_label, partition_id, gmg_level);
   auto it = containers_.find(label);
   if (it == containers_.end()) {
     // TODO(someone) add caching of partitions to Mesh at some point
     const int pack_size = pmy_mesh_->DefaultPackSize();
+    auto &block_list =
+        gmg_level ? pmy_mesh_->gmg_block_lists[*gmg_level] : pmy_mesh_->block_list;
     auto partitions = partition::ToSizeN(block_list, pack_size);
     // Account for possibly empty block_list
     if (partitions.size() == 0) partitions = std::vector<BlockList_t>(1);
     for (auto i = 0; i < partitions.size(); i++) {
-      std::string md_label = mbd_label + "_part-" + std::to_string(i);
-      if (gmg_level) md_label = md_label + "_gmg-" + std::to_string(*gmg_level);
+      std::string md_label = GetKey(mbd_label, partition_id, gmg_level);
       containers_[md_label] = std::make_shared<MeshData<Real>>(mbd_label);
-      containers_[md_label]->Set(partitions[i], pmy_mesh_);
-      if (gmg_level) {
-        containers_[md_label]->grid = GridIdentifier::two_level_composite(*gmg_level);
-      } else {
-        containers_[md_label]->grid = GridIdentifier::leaf();
-      }
+      containers_[md_label]->Initialize(partitions[i], pmy_mesh_, gmg_level);
+      containers_[md_label]->partition = i;
     }
   }
   return containers_[label];
@@ -65,16 +61,14 @@ template <>
 std::shared_ptr<MeshData<Real>> &
 DataCollection<MeshData<Real>>::GetOrAdd(const std::string &mbd_label,
                                          const int &partition_id) {
-  return GetOrAdd_impl(pmy_mesh_, containers_, pmy_mesh_->block_list, mbd_label,
-                       partition_id, {});
+  return GetOrAdd_impl(mbd_label, partition_id, {});
 }
 
 template <>
 std::shared_ptr<MeshData<Real>> &
 DataCollection<MeshData<Real>>::GetOrAdd(int gmg_level, const std::string &mbd_label,
                                          const int &partition_id) {
-  return GetOrAdd_impl(pmy_mesh_, containers_, pmy_mesh_->gmg_block_lists[gmg_level],
-                       mbd_label, partition_id, gmg_level);
+  return GetOrAdd_impl(mbd_label, partition_id, gmg_level);
 }
 
 template class DataCollection<MeshData<Real>>;