Merge branch 'develop' into lroberts36/add-combined-buffer-communication

.github/workflows/docs.yml

@@ -12,9 +12,12 @@ on:
 jobs:
     docs:
       name: build and deploy docs
-      runs-on: ubuntu-latest
+      # not using latest due to issues with pip user packages, see
+      # https://github.com/actions/runner-images/issues/10781 and
+      # https://github.com/actions/runner-images/issues/10636
+      runs-on: ubuntu-22.04
 
-      steps:      
+      steps:
         - name: Checkout code
           uses: actions/checkout@v2
           with:
@@ -23,9 +26,9 @@ jobs:
           run: export DEBIAN_FRONTEND=noninteractive
         - name: install dependencies
           run: |
-            pip install --break-system-packages sphinx
-            pip install --break-system-packages sphinx-rtd-theme
-            pip install --break-system-packages sphinx-multiversion
+            pip install sphinx
+            pip install sphinx-rtd-theme
+            pip install sphinx-multiversion
         - name: build docs
           run: |
             echo "Repo = ${GITHUB_REPOSITORY}"

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1103]](https://github.com/parthenon-hpc-lab/parthenon/pull/1103) Add sparsity to vector wave equation test
 - [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185/files) Bugfix to particle defragmentation
 - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D
 - [[PR 1183]](https://github.com/parthenon-hpc-lab/parthenon/pull/1183) Fix particle leapfrog example initialization data
@@ -11,12 +12,17 @@
 - [[PR 1161]](https://github.com/parthenon-hpc-lab/parthenon/pull/1161) Make flux field Metadata accessible, add Metadata::CellMemAligned flag, small perfomance upgrades
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 1206]](https://github.com/parthenon-hpc-lab/parthenon/pull/1206) Leapfrog fix
+- [[PR1203]](https://github.com/parthenon-hpc-lab/parthenon/pull/1203) Pin Ubuntu CI image
+- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag
 - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize
-- [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1
+- [[Issue 1165]](https://github.com/parthenon-hpc-lab/parthenon/issues/1165) Bump Kokkos submodule to 4.4.1
 - [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option
 - [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1170]](https://github.com/parthenon-hpc-lab/parthenon/pull/1170) Fixed incorrect initialization of array by a const not constexpr
+- [[PR 1189]](https://github.com/parthenon-hpc-lab/parthenon/pull/1189) Address CUDA MPI/ICP issue with Kokkos <=4.4.1
 - [[PR 1178]](https://github.com/parthenon-hpc-lab/parthenon/pull/1178) Fix issue with mesh pointer when using relative residual tolerance in BiCGSTAB solver.
 - [[PR 1173]](https://github.com/parthenon-hpc-lab/parthenon/pull/1173) Make debugging easier by making parthenon throw an error if ParameterInput is different on multiple MPI ranks.
 
@@ -27,7 +33,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
-
+- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag
 
 ## Release 24.08
 Date: 2024-08-30

CMakeLists.txt

@@ -291,6 +291,13 @@ if (Kokkos_ENABLE_CUDA)
   if(CHECK_REGISTRY_PRESSURE)
     add_compile_options(-Xptxas=-v)
   endif()
+
+  # Async malloc are enabled by default until Kokkos <= 4.4.1 but
+  # cause issues with IPCs and MPI,
+  # see https://github.com/parthenon-hpc-lab/athenapk/pull/114#issuecomment-2358981857
+  # Given that it's also not expected to be a significant performance gain
+  # it's hard disabled now.
+  set(Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC OFF CACHE BOOL "Disable Cuda async malloc (to address MPI/IPC issues)")
 endif()
 # Note that these options may not play nice with nvcc wrapper
 if (TEST_INTEL_OPTIMIZATION)

doc/sphinx/src/boundary_conditions.rst

@@ -1,4 +1,4 @@
-.. _sphinx-doc:
+.. _boundary-conds:
 
 Boundary Conditions
 ===================
@@ -10,7 +10,6 @@ Natively, Parthenon supports three kinds of boundary conditions:
 
 - ``periodic``
 - ``outflow``
-- ``reflecting``
 
 which are all imposed on variables with the ``Metadata::FillGhost``
 metadata flag. To set the boundaries in each direction, set the
@@ -22,8 +21,8 @@ metadata flag. To set the boundaries in each direction, set the
    ix1_bc = outflow
    ox1_bc = outflow
 
-   ix2_bc = reflecting
-   ox2_bc = reflecting
+   ix2_bc = outflow
+   ox2_bc = outflow
 
    ix3_bc = periodic
    ox3_bc = periodic
@@ -40,7 +39,9 @@ for your ``parthenon_manager``. e.g.,
 
 .. code:: c++
 
-   pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x1] = MyBoundaryInnerX1;
+   pman.app_input->RegisterBoundaryCondition(
+	  parthenon::BoundaryFace::inner_x1,
+	  "my_bc_name", MyBoundaryInnerX1);
 
 where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as
 
@@ -58,13 +59,13 @@ where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as
      outer_x3 = 5
    };
 
-You can then set this boundary condition via the ``user`` flag in the
-input file:
+You can then set this boundary condition by using the name you
+registered in the input file:
 
 ::
 
    <parthenon/mesh>
-   ix1_bc = user
+   ix1_bc = my_bc_name
 
 Boundary conditions so defined should look roughly like
 
@@ -100,9 +101,19 @@ Other than these requirements, the ``Boundary`` object can do whatever
 you like. Reference implementations of the standard boundary conditions
 are available `here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop/src/bvals/boundary_conditions.cpp>`__.
 
+.. note::
 
-Per package user-defined boundary conditions
---------------------------------------------
+  A per-variable reflecting boundary condition is available, but you
+  must register it manually. To do so, simply call
+  ``app_in->RegisterDefaultReflectingBoundaryConditions()`` and it
+  will be available as a mesh boundary with the name ``reflecting``.
+  The reason manual registration is required is to support custom
+  reflecting boundary conditions in the case where a single variable
+  is used as the state vector.
+
+
+Per package user-defined boundary conditions.
+---------------------------------
 
 In addition to user defined *global* boundary conditions, Parthenon also supports 
 registration of boundary conditions at the *per package* level. These per package 
@@ -133,4 +144,4 @@ for a more complete example):
     pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetMyBC<X1DIR, BCSide::Inner>());
     pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetMyBC<X2DIR, BCSide::Inner>());
     ...
-  }
+  }

example/advection/main.cpp

@@ -25,6 +25,7 @@ int main(int argc, char *argv[]) {
   pman.app_input->ProblemGenerator = advection_example::ProblemGenerator;
   pman.app_input->UserWorkAfterLoop = advection_example::UserWorkAfterLoop;
   pman.app_input->UserMeshWorkBeforeOutput = advection_example::UserMeshWorkBeforeOutput;
+  pman.app_input->RegisterDefaultReflectingBoundaryConditions();
 
   // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
   auto manager_status = pman.ParthenonInitEnv(argc, argv);

example/boundary_exchange/boundary_exchange_driver.cpp

@@ -76,23 +76,23 @@ int main(int argc, char *argv[]) {
                      ar3_t{1.0, 1.0, 1.0});
   // forest_def.AddFace(0, {n[0], n[1], n[3], n[2]}, ar3_t{0.0, 0.0, 0.0},
   // ar3_t{1.0, 1.0, 1.0});
-  forest_def.AddBC(edge_t({n[0], n[1]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[0], n[3]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[0], n[1]}));
+  forest_def.AddBC(edge_t({n[0], n[3]}));
 
   forest_def.AddFace(1, {n[1], n[4], n[2], n[5]}, ar3_t{2.0, 0.0, 0.0},
                      ar3_t{3.0, 1.0, 1.0});
-  forest_def.AddBC(edge_t({n[1], n[4]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[4], n[5]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[1], n[4]}));
+  forest_def.AddBC(edge_t({n[4], n[5]}));
 
   forest_def.AddFace(3, {n[3], n[2], n[6], n[7]}, ar3_t{0.0, 2.0, 0.0},
                      ar3_t{1.0, 3.0, 1.0});
-  forest_def.AddBC(edge_t({n[6], n[7]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[3], n[6]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[6], n[7]}));
+  forest_def.AddBC(edge_t({n[3], n[6]}));
 
   forest_def.AddFace(4, {n[2], n[5], n[7], n[8]}, ar3_t{2.0, 2.0, 0.0},
                      ar3_t{3.0, 3.0, 1.0});
-  forest_def.AddBC(edge_t({n[5], n[8]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[7], n[8]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[5], n[8]}));
+  forest_def.AddBC(edge_t({n[7], n[8]}));
 
   forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, 1, 0, 0, 0));
   pman.ParthenonInitPackagesAndMesh(forest_def);

example/fine_advection/advection_package.cpp

@@ -225,9 +225,9 @@ TaskStatus FillDerived(MeshData<Real> *md) {
           const int fi = (ndim > 0) ? (i - nghost) * 2 + nghost : i;
           pack(b, Conserved::phi_fine_restricted(), k, j, i) = 0.0;
           Real ntot = 0.0;
-          for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+          for (int koff = 0; koff <= (ndim > 2); ++koff)
             for (int joff = 0; joff <= (ndim > 1); ++joff)
-              for (int koff = 0; koff <= (ndim > 2); ++koff) {
+              for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                 ntot += 1.0;
                 pack(b, Conserved::phi_fine_restricted(), k, j, i) +=
                     pack(b, Conserved::phi_fine(), kf + koff, jf + joff, fi + ioff);

example/fine_advection/advection_package.hpp

@@ -80,11 +80,12 @@ TaskStatus CalculateFluxes(pack_desc_t &desc, parthenon::TopologicalElement FACE
   IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior, FACE);
   IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior, FACE);
   IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior, FACE);
-  constexpr int scratch_size = 0;
-  constexpr int scratch_level = 1;
+  constexpr int unused_scratch_size = 0;
+  constexpr int unused_scratch_level = 1;
   parthenon::par_for_outer(
-      PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, kb.s,
-      kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
+      PARTHENON_AUTO_LABEL, unused_scratch_size, unused_scratch_level, 0,
+      pack.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});
         for (int l = pack.GetLowerBound(b); l <= pack.GetUpperBound(b); ++l) {
           parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) {

example/fine_advection/parthenon_app_inputs.cpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-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
@@ -119,26 +119,26 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
             Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) +
                        coords.Xc<2>(j) * coords.Xc<2>(j) +
                        coords.Xc<3>(k) * coords.Xc<3>(k);
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) =
                       1. + amp * exp(-100.0 * rsq);
                 }
           } else if (profile_type == 2) {
             Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) +
                        coords.Xc<2>(j) * coords.Xc<2>(j) +
                        coords.Xc<3>(k) * coords.Xc<3>(k);
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) =
                       (rsq < 0.15 * 0.15 ? 1.0 : 0.0);
                 }
           } else {
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = 0.0;
                 }
           }

example/fine_advection/parthinput.advection

@@ -57,6 +57,13 @@ profile = hard_sphere
 refine_tol = 0.3    # control the package specific refinement tagging function
 derefine_tol = 0.03
 
+# Include advection equation for regular resolution CC fields
+do_regular_advection = true
+# Include advection equation for fine resolution CC fields
+do_fine_advection = true
+# Include vector wave equation using CT
+do_CT_advection = true
+
 <parthenon/output1>
 file_type = rst
 dt = 0.05

example/particles/main.cpp

@@ -49,12 +49,14 @@ int main(int argc, char *argv[]) {
 
   // Redefine parthenon defaults
   pman.app_input->ProcessPackages = particles_example::ProcessPackages;
-  pman.app_input->boundary_conditions[BoundaryFace::inner_x1] =
-      parthenon::BoundaryFunction::OutflowInnerX1;
-  pman.app_input->boundary_conditions[BoundaryFace::outer_x1] =
-      parthenon::BoundaryFunction::OutflowOuterX1;
-  pman.app_input->swarm_boundary_conditions[BoundaryFace::inner_x1] = SwarmUserInnerX1;
-  pman.app_input->swarm_boundary_conditions[BoundaryFace::outer_x1] = SwarmUserOuterX1;
+  pman.app_input->RegisterBoundaryCondition(BoundaryFace::inner_x1,
+                                            parthenon::BoundaryFunction::OutflowInnerX1);
+  pman.app_input->RegisterBoundaryCondition(BoundaryFace::outer_x1,
+                                            parthenon::BoundaryFunction::OutflowOuterX1);
+  pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1,
+                                                 SwarmUserInnerX1);
+  pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1,
+                                                 SwarmUserOuterX1);
   // Note that this example does not use a ProblemGenerator
   pman.ParthenonInitPackagesAndMesh();
 

example/particles/particles.cpp

@@ -288,17 +288,18 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
           } while (r > 0.5);
 
           // Randomly sample direction perpendicular to origin
-          Real theta = acos(2. * rng_gen.drand() - 1.);
-          Real phi = 2. * M_PI * rng_gen.drand();
-          v(0, n) = sin(theta) * cos(phi);
-          v(1, n) = sin(theta) * sin(phi);
-          v(2, n) = cos(theta);
+          const Real mu = 2.0 * rng_gen.drand() - 1.0;
+          const Real phi = 2. * M_PI * rng_gen.drand();
+          const Real stheta = std::sqrt(1.0 - mu * mu);
+          v(0, n) = stheta * cos(phi);
+          v(1, n) = stheta * sin(phi);
+          v(2, n) = mu;
           // Project v onto plane normal to sphere
           Real vdN = v(0, n) * x(n) + v(1, n) * y(n) + v(2, n) * z(n);
-          Real NdN = r * r;
-          v(0, n) = v(0, n) - vdN / NdN * x(n);
-          v(1, n) = v(1, n) - vdN / NdN * y(n);
-          v(2, n) = v(2, n) - vdN / NdN * z(n);
+          Real inverse_NdN = 1. / (r * r);
+          v(0, n) = v(0, n) - vdN * inverse_NdN * x(n);
+          v(1, n) = v(1, n) - vdN * inverse_NdN * y(n);
+          v(2, n) = v(2, n) - vdN * inverse_NdN * z(n);
 
           // Normalize
           Real v_tmp = sqrt(v(0, n) * v(0, n) + v(1, n) * v(1, n) + v(2, n) * v(2, n));
@@ -335,11 +336,12 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
           z(n) = minx_k + nx_k * dx_k * rng_gen.drand();
 
           // Randomly sample direction on the unit sphere, fixing speed
-          Real theta = acos(2. * rng_gen.drand() - 1.);
-          Real phi = 2. * M_PI * rng_gen.drand();
-          v(0, n) = vel * sin(theta) * cos(phi);
-          v(1, n) = vel * sin(theta) * sin(phi);
-          v(2, n) = vel * cos(theta);
+          const Real mu = 2.0 * rng_gen.drand() - 1.0;
+          const Real phi = 2. * M_PI * rng_gen.drand();
+          const Real stheta = std::sqrt(1.0 - mu * mu);
+          v(0, n) = vel * stheta * cos(phi);
+          v(1, n) = vel * stheta * sin(phi);
+          v(2, n) = vel * mu;
 
           // Create particles at the beginning of the timestep
           t(n) = t0;

example/sparse_advection/main.cpp

@@ -27,6 +27,7 @@ int main(int argc, char *argv[]) {
   pman.app_input->ProblemGenerator = sparse_advection_example::ProblemGenerator;
   pman.app_input->PostStepDiagnosticsInLoop =
       sparse_advection_example::PostStepDiagnosticsInLoop;
+  pman.app_input->RegisterDefaultReflectingBoundaryConditions();
 
   // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
   auto manager_status = pman.ParthenonInitEnv(argc, argv);

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

@@ -211,14 +211,17 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
         if one:
             return ret_code
 
-    # Compare the names of attributes in /Info, except "Time"
+    # Compare the names of attributes in /Info, except those we know
+    # may vary safely
     f0_Info = {
         key: value
         for key, value in f0.Info.items()
         if key != "Time"
         and key != "BlocksPerPE"
         and key != "WallTime"
         and key != "OutputFormatVersion"
+        and key != "BoundaryConditions"
+        and key != "SwarmBoundaryConditions"
     }
     f1_Info = {
         key: value
@@ -227,6 +230,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
         and key != "BlocksPerPE"
         and key != "WallTime"
         and key != "OutputFormatVersion"
+        and key != "BoundaryConditions"
+        and key != "SwarmBoundaryConditions"
     }
     if sorted(f0_Info.keys()) != sorted(f1_Info.keys()):
         print("Names of attributes in '/Info' of differ")

src/CMakeLists.txt

@@ -95,6 +95,9 @@ set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are
 configure_file(config.hpp.in generated/config.hpp @ONLY)
 
 add_library(parthenon
+  application_input.cpp
+  application_input.hpp
+
   bvals/comms/bvals_in_one.hpp
   bvals/comms/bvals_utils.hpp
   bvals/comms/build_boundary_buffers.cpp