Bump Kokkos to 5.1.1 and Parthenon to upcoming 26.xx version (#167)

* Bump Kokkos to 5.1.1 and Parthenon to upcoming 26.xx version

* Remove incorrect enforced symmetry with race condition and add explicit checks

* Remove complex conjugate modes from sample input

* Make REQUIRE a WARN for backwards compatibility

* Remove Hermitian modes from gen fmturb modes

* Add Changelog

* Update reference data given new sample vectors

* use ref data created on CI machine

* Update ref data following latest PArtheon fixes

* Enforce stronger limits on error in particle advection

* Add info on breaking changes

* Fix compiluation given updated Parthenon API

* Use ULP logic for comparison in test

* Make formatter happy

* Bump parth to use standard compliant labels

* Fix parth compile and print warn only on rank 0

* Update relevant info in changelog

Author: pgrete

CHANGELOG.md

@@ -3,6 +3,23 @@
 ## Current develop (i.e., `main` branch)
 
 ### General notes
+With the latest update of the Parthenon submodule several new features are now available, e.g.,
+- [OpenPMD output](https://parthenon-hpc-lab.github.io/parthenon/pgrete/pmd-output/src/outputs.html#openpmd) including support for slices, data compression and coarse graining\
+Note, the original naming convention (in the development branch) used for labeling components was not fully standard compliant.
+New outputs are automatically written in the standard compliant version.
+However, to keep the old scheme (e.g., when restarting from existing simulation data with the intent to keep the timeseries consistent add `openpmd_format_version=1` to the corresponding openpmd output blocks.
+- A watchdog (to kill a simulation that hangs for whatever reason). Just run with `-w HH:MM:SS`.
+- Support for (tracer) particles with AMR.
+- Support for using more parallelism for the ghost zone exchange kernels, which speeds up simulation with few (<10) blocks per rank on device, see [here](https://github.com/parthenon-hpc-lab/parthenon/pull/1271).
+A good starting point is `parthenon/mesh/minimum_number_of_teams_for_boundary_kernel=8`.
+
+With the update of Kokkos to version 5.1.1 (or >5.0 in general) a performance regression was identified.
+This is likely related to the Kokkos-internal use of int64 indices in the new `mdspan` based `View`s resulting more register usage (which. in turn, results in lower occupancy on devices).
+The Kokkos team is aware of this and working on a fix.
+If the current performance is (significantly) below expectation, one can try to use "legacy" views via `Kokkos_ENABLE_IMPL_VIEW_LEGACY=ON`.
+
+**IMPORTANT** The latest Parthenon submodule includes a fix for a race condition in the flux correction communication routine when run with mesh refinement on GPUs, see [here](https://github.com/parthenon-hpc-lab/parthenon/pull/1405).
+Please update immediately or rebuild AthenaPK with `PARTHENON_DISABLE_SPARSE=OFF` to mitigate the race condition.
 
 ### Added (new features/APIs/variables/...)
 
@@ -12,11 +29,14 @@
 - [[PR 172]](https://github.com/parthenon-hpc-lab/athenapk/pull/172) Fixed data race condition in few modes IFT (no practical implication)
 
 ### Infrastructure
+- [[PR 167]](https://github.com/parthenon-hpc-lab/athenapk/pull/167) Bump Kokkos to 5.1.1 and `Parthenon` to upcoming 26.xx version (incl OpenPMD support)
 - [[PR 168]](https://github.com/parthenon-hpc-lab/athenapk/pull/168) Document agentic coding guidelines (and add PR template)
 
 ### Removed (removing behavior/API/varaibles/...)
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 167]](https://github.com/parthenon-hpc-lab/athenapk/pull/167) C++20 is now the minimum standard
+- [[PR 167]](https://github.com/parthenon-hpc-lab/athenapk/pull/167) New Parthenon submodule changed input file parsing (removed `*pib = pin->pfirst_block;`), see [here](https://github.com/parthenon-hpc-lab/parthenon/pull/1385)
 
 ## Release 26.05
 

CMakeLists.txt

@@ -19,7 +19,7 @@ project(athenaPK VERSION 26.05 LANGUAGES CXX)
 
 set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 
 # make paths in target_source() absolute
 # used for sources in src/pgen folder
@@ -70,6 +70,7 @@ set(PARTHENON_ENABLE_TESTING OFF CACHE BOOL "Disable Parthenon testing.")
 set(PARTHENON_DISABLE_OPENMP ON CACHE BOOL "Disable OpenMP")
 set(PARTHENON_DISABLE_EXAMPLES ON CACHE BOOL "Don't build Parthenon examples.")
 set(PARTHENON_DISABLE_SPARSE ON CACHE BOOL "Disable sparse (not used in AthenaPK yet)")
+set(PARTHENON_DISABLE_OPENPMD OFF CACHE BOOL "Enable OpenPMD outputs in AthenaPK by default")
 
 if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/external/parthenon/CMakeLists.txt)
   add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/parthenon parthenon)

README.md

@@ -58,7 +58,7 @@ please either
 ##### Required
 
 * CMake 3.13 or greater
-* C++17 compatible compiler
+* C++20 compatible compiler
 * Parthenon (using the submodule version provided by AthenaPK)
 * Kokkos (using the submodule version provided by AthenaPK)
 
@@ -67,6 +67,7 @@ please either
 * MPI
 * OpenMP (for host parallelism. Note that MPI is the recommended option for on-node parallelism.)
 * HDF5 (for outputs)
+* OpenPMD and ADIOS2 (for outputs, recommended future standard but no `yt` support yet. Can be disabled via `PARTHENON_DISABLE_OPENPMD=ON` during configure)
 * Python3 (for regressions tests with numpy, scipy, matplotlib, unyt, and h5py modules)
 * Ascent (for in situ visualization and analysis)
 
@@ -125,22 +126,36 @@ Some example input files are provided in the [inputs](inputs/) folder.
 
 #### Data Analysis
 
-There exit several options to read/process data written by AthenaPK -- specifically in
-the `file_type = hdf5` format, see
+There exit several options to read/process data written by AthenaPK, see
 [Parthenon doc](https://parthenon-hpc-lab.github.io/parthenon/develop/src/outputs.html):
 
 1. With [ParaView](https://www.paraview.org/) and
 [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/).
-In ParaView, select the "XDMF Reader" when prompted.
+For HDF5 in ParaView, select the "XDMF Reader" when prompted.
+For openPMD outputs, the native BP5 reader can be used to open the files.
+However, this may perform poorly (due to lack of automated mesh decomposition
+and load balancing).
+Thus, it is recommended to use the specific
+[openpmd-visit-reader](https://github.com/BenWibking/openpmd-visit-reader)
+plugin for VisIt, which has robustly handled larger, multiple TB datasets.
 
 2. With [yt](https://yt-project.org/)
 As of versions >=4.4 `*.phdf` files can be read as usual with `yt.load()`.
+OpenPMD/ADIOS2 output cannot be properly read by `yt` yet.
 
-3. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis).
+3. OpenPMD outputs in general.
+Outputs written following the openPMD standard can be processed by any tool implemting
+the standard.
+For Python based analysis the most straightforward way is via the
+[openPMD-api](https://github.com/openPMD/openPMD-api)
+or even more easily via the [openPMD-viewer](https://github.com/openPMD/openPMD-viewer)
+that builds on the API (both can be installed via `pip`).
+
+4. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis).
 This requires Ascent to be installed/available at compile time of AthenaPK.
 To enable set `PARTHENON_ENABLE_ASCENT=ON`.
 
-4. (Not recommended) Using the integrated Python script called "`phdf`" provided by Parthenon,
+5. (Not recommended) Using the integrated Python script called "`phdf`" provided by Parthenon,
 i.e., the either install `parthenon_tools`
 (located in `external/parthenon/scripts/python/packages/parthenon/tools`) or add
 that directory to your Python path.

external/Kokkos

@@ -119,13 +119,10 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
     // rescale sim time limit
     pin->SetReal("parthenon/time", "tlim", tlim_rescaled);
     // rescale dt of each output block
-    parthenon::InputBlock *pib = pin->pfirst_block;
-    while (pib != nullptr) {
-      if (pib->block_name.compare(0, 16, "parthenon/output") == 0) {
-        auto dt = pin->GetReal(pib->block_name, "dt");
-        pin->SetReal(pib->block_name, "dt", dt * t_cc);
-      }
-      pib = pib->pnext; // move to next input block name
+    auto blocks = pin->GetBlockNamesWithPrefix("parthenon/output");
+    for (const auto &block_name : blocks) {
+      auto dt = pin->GetReal(block_name, "dt");
+      pin->SetReal(block_name, "dt", dt * t_cc);
     }
 
     msg << "## Interpreted time limits (partenon/time/tlim and dt for outputs) as in "

external/parthenon

@@ -163,13 +163,10 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
     // rescale sim time limit
     pin->SetReal("parthenon/time", "tlim", tlim_orig * tlim);
     // rescale dt of each output block
-    parthenon::InputBlock *pib = pin->pfirst_block;
-    while (pib != nullptr) {
-      if (pib->block_name.compare(0, 16, "parthenon/output") == 0) {
-        auto dt = pin->GetReal(pib->block_name, "dt");
-        pin->SetReal(pib->block_name, "dt", dt * tlim);
-      }
-      pib = pib->pnext; // move to next input block name
+    auto blocks = pin->GetBlockNamesWithPrefix("parthenon/output");
+    for (const auto &block_name : blocks) {
+      auto dt = pin->GetReal(block_name, "dt");
+      pin->SetReal(block_name, "dt", dt * tlim);
     }
 
     // Now disable rescaling of times so that this is done only once and not for restarts

src/pgen/cloud.cpp

@@ -1,6 +1,6 @@
 //========================================================================================
 // AthenaPK - a performance portable block structured AMR astrophysical MHD code.
-// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved.
+// Copyright (c) 2023-2026, Athena-Parthenon Collaboration. All rights reserved.
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
 //========================================================================================
@@ -74,17 +74,19 @@ FewModesFT::FewModesFT(parthenon::ParameterInput *pin, parthenon::StateDescripto
       if (k_vec_host(0, i) == 0 && k_vec_host(0, i) == k_vec_host(0, j) &&
           k_vec_host(1, i) == -k_vec_host(1, j) &&
           k_vec_host(2, i) == -k_vec_host(2, j)) {
-        PARTHENON_WARN(
-            "The given set of k_vec include complex conjugate partners for mode " +
-            std::to_string(i) + " with components " + std::to_string(k_vec_host(0, i)) +
-            " " + std::to_string(k_vec_host(1, i)) + " " +
-            std::to_string(k_vec_host(2, i)) +
-            "."
-            "In theory, this results in these modes being counted double. In practice, "
-            "this will have little to no effect as the normalization is done separate "
-            "for the real space field rather than the spectral field. However, if this "
-            "is a fresh simulation (and not a restarted one with given/fixed k_vec) it "
-            "is recommended to update the set of k_vec in the input file.");
+        if (parthenon::Globals::my_rank == 0) {
+          PARTHENON_WARN(
+              "The given set of k_vec include complex conjugate partners for mode " +
+              std::to_string(i) + " with components " + std::to_string(k_vec_host(0, i)) +
+              " " + std::to_string(k_vec_host(1, i)) + " " +
+              std::to_string(k_vec_host(2, i)) +
+              "."
+              "In theory, this results in these modes being counted double. In practice, "
+              "this will have little to no effect as the normalization is done separate "
+              "for the real space field rather than the spectral field. However, if this "
+              "is a fresh simulation (and not a restarted one with given/fixed k_vec) it "
+              "is recommended to update the set of k_vec in the input file.");
+        }
       }
       PARTHENON_REQUIRE_THROWS(!(k_vec_host(0, i) == k_vec_host(0, j) &&
                                  k_vec_host(1, i) == k_vec_host(1, j) &&

src/pgen/shattering.cpp

@@ -95,8 +95,7 @@ def Analyse(self, parameters):
             relabs = 0.5 * np.abs(((a - b) + 1) % 2 - 1)
 
             # In principle, the difference should be 0
-            # TODO(PG) investigate where this comes from. Integrator? Interpolation? ...?
-            if relabs.max() > 0.003:
+            if relabs.max() > 1e-14:
                 print(
                     f"ERROR: difference between intial and final position to large for {pos}: {relabs.max()}"
                 )