Make a global variable for whether simulation is a restart (#1179)

* Set is_restart

* Add soft equivalence

* soft_equiv used for testing

* CHANGELOG, test_kokkos_abstraction

* copyright

* Oops wrong var

* Oops comparing error and reference

* Update src/utils/robust.hpp

Co-authored-by: Jonah Miller <jonah.maxwell.miller@gmail.com>

* style

* naming

* Make code more self-documenting

---------

Co-authored-by: Jonah Miller <jonah.maxwell.miller@gmail.com>

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1179]](https://github.com/parthenon-hpc-lab/parthenon/pull/1179) Make a global variable for whether simulation is a restart
 - [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option
 - [[PR 1161]](https://github.com/parthenon-hpc-lab/parthenon/pull/1161) Make flux field Metadata accessible, add Metadata::CellMemAligned flag, small perfomance upgrades
 

src/globals.cpp

@@ -3,7 +3,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (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
@@ -30,8 +30,9 @@ namespace Globals {
 int nghost;
 
 // all of these global variables are set at the start of main():
-int my_rank; // MPI rank of this process
-int nranks;  // total number of MPI ranks
+int my_rank;     // MPI rank of this process
+int nranks;      // total number of MPI ranks
+bool is_restart; // Whether this simulation is restarted from a checkpoint file
 
 // sparse configuration values that are needed in various places
 SparseConfig sparse_config;

src/globals.hpp

@@ -3,7 +3,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (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
@@ -36,6 +36,7 @@ struct SparseConfig {
 };
 
 extern int my_rank, nranks, nghost;
+extern bool is_restart;
 
 extern SparseConfig sparse_config;
 

src/parthenon_manager.cpp

@@ -70,6 +70,8 @@ ParthenonStatus ParthenonManager::ParthenonInitEnv(int argc, char *argv[]) {
   Globals::nranks = 1;
 #endif // MPI_PARALLEL
 
+  Globals::is_restart = IsRestart();
+
   Kokkos::initialize(argc, argv);
 
   // pgrete: This is a hack to disable allocation tracking until the Kokkos

src/utils/robust.hpp

@@ -60,6 +60,18 @@ KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) {
   const B sgn = b >= 0 ? 1 : -1;
   return a / (b + sgn * SMALL<B>());
 }
+
+// Return true equivalence if value and reference differ by less than precision
+// Optionally return true if reference value is close to zero
+KOKKOS_FORCEINLINE_FUNCTION
+bool SoftEquiv(const Real &val, const Real &ref,
+               const Real eps = 10. * std::numeric_limits<Real>::epsilon(),
+               const bool pass_on_small = true) {
+  const bool is_close = std::abs(val - ref) < eps * std::abs(ref);
+  const bool is_small = std::abs(ref) < std::numeric_limits<Real>::min();
+  return (is_close || (is_small && pass_on_small));
+}
+
 } // namespace robust
 } // namespace parthenon
 #endif // UTILS_ROBUST_HPP_

tst/unit/CMakeLists.txt

@@ -25,7 +25,7 @@ list(APPEND unit_tests_SOURCES
     test_unit_constants.cpp
     test_unit_domain.cpp
     test_unit_sort.cpp
-    kokkos_abstraction.cpp
+    test_kokkos_abstraction.cpp
     test_index_split.cpp
     test_logical_location.cpp
     test_forest.cpp

tst/unit/kokkos_abstraction.cpp → tst/unit/test_kokkos_abstraction.cpp

@@ -26,12 +26,14 @@
 
 #include "basic_types.hpp"
 #include "kokkos_abstraction.hpp"
+#include "utils/robust.hpp"
 
 using parthenon::DevExecSpace;
 using parthenon::ParArray1D;
 using parthenon::ParArray2D;
 using parthenon::ParArray3D;
 using parthenon::ParArray4D;
+using parthenon::robust::SoftEquiv;
 using Real = double;
 
 template <class T>
@@ -316,22 +318,22 @@ bool test_wrapper_nested_3d(OuterLoopPattern outer_loop_pattern,
   // Copy array back from device to host
   Kokkos::deep_copy(host_du, dev_du);
 
-  Real max_rel_err = -1;
   const Real rel_tol = std::numeric_limits<Real>::epsilon();
 
   // compare data on the host
   for (int k = 0; k < N; k++) {
     for (int j = 0; j < N; j++) {
       for (int i = 1; i < N - 1; i++) {
         const Real analytic = 2.0 * (i + 1) * pow((j + 2) * (k + 3), 2.0);
-        const Real err = host_du(k, j, i - 1) - analytic;
 
-        max_rel_err = fmax(fabs(err / analytic), max_rel_err);
+        if (!SoftEquiv(host_du(k, j, i - 1), analytic, rel_tol)) {
+          return false;
+        }
       }
     }
   }
 
-  return max_rel_err < rel_tol;
+  return true;
 }
 
 template <class OuterLoopPattern, class InnerLoopPattern>
@@ -385,7 +387,6 @@ bool test_wrapper_nested_4d(OuterLoopPattern outer_loop_pattern,
   // Copy array back from device to host
   Kokkos::deep_copy(host_du, dev_du);
 
-  Real max_rel_err = -1;
   const Real rel_tol = std::numeric_limits<Real>::epsilon();
 
   // compare data on the host
@@ -394,15 +395,16 @@ bool test_wrapper_nested_4d(OuterLoopPattern outer_loop_pattern,
       for (int j = 0; j < N; j++) {
         for (int i = 1; i < N - 1; i++) {
           const Real analytic = 2.0 * (i + 1) * pow((j + 2) * (k + 3) * (n + 4), 2.0);
-          const Real err = host_du(n, k, j, i - 1) - analytic;
 
-          max_rel_err = fmax(fabs(err / analytic), max_rel_err);
+          if (!SoftEquiv(host_du(n, k, j, i - 1), analytic, rel_tol)) {
+            return false;
+          }
         }
       }
     }
   }
 
-  return max_rel_err < rel_tol;
+  return true;
 }
 
 TEST_CASE("nested par_for loops", "[wrapper]") {