Feature/subcycle dev (#1)

* CI: Update to AMReX 24.09

* Arith: Implement efficient SIMD trigonometic functions

* Arith: Correct SIMD tests

* Arith: Correct SIMD tests

* Arith: Remove `constexpr` that some compilers don't like

* Arith: Remove `constexpr` that some compilers don't like

* TestDerivs: Correct C++ calls to trigonometic functions

* CI: Disable SIMD the new way

* CI: Update Docker images

* CI: Correct Docker image names

* CarpetX: Update to nvtx3

* CI: Reduce optimization level for OneAPI

* CI: Switch to production images

* CarpetXRegrid: Move regrid_error to a new thorn.

* CI: Add CarpetXRegrid to thorn list

* CarpetX: Change the implementation to Driver

This involves adding a parameter "periodic" whose default value is "no."
It also involves modifying parameter files, file names, and source files
where the name "CarpetX" needs to be replaced by "Driver."

* Algo: remove parameter file and changing data from test output

* BenchX: sort thorns alphabetically

* BenchX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* BoxInBox: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* ErrorEstimator: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* FluxWaveToyX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* MovingBoxToy: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* PoissonX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* SiMDWaveToyX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* SpacetimeWaveToyX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* StaggeredWaveToyX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestArrayGroup: remove parameter file and changing data from test output

* TestBoundaries: adapt to moving regrid error into its own thorn

* TestBoxInBox: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestDerivs: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestInterpolate: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestLoopX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestNorms: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestODESolvers: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestODESolvers2: remove parameter file and changing data from test output

* TestOutput: adapt to moving regrid error into its own thorn

* TestProlongate: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* TestSymmetries: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* WaveToyX: adapt to CarpetX implementation change

required due to CarpetX implemting "Driver"

* ODESolvers: Schedule `ODESolvers_PostStep` after recovering

* Loop: Update AMReX include files

* ADMBaseX: Clean up schedule

* HydroBaseX: Clean up schedule

* CarpetX: Correct ENO operators

* CarpetX: Correct ENO operators

* CarpetX: Correct poison check on HIP systems

* CarpetX: Implement `minmod` prolongation operator

* CarpetX: Support openPMD_api 0.16.0

* CarpetX: Use Blosc compression by default for openPMD output

* CarpetX: Reformat code

* CarpetX: ADIOS2: Avoid too much screen output

* CarpetX: ADIOS2: Use Blosc compression by default

* Arith: Clean up code

* BoxInBox: Don't allow zero radius

* Loop: Inline more aggressively

* MovingBoxToy: Code cleanup

* TestSymmetries: Beautify code

* Loop: New macro CCTK_KERNEL

* CI: Update Docker images

* FluxWaveToyX: Correct RHS handling after ODESolver update

* SIMDWaveToyX: Correct RHS handling after ODESolver update

* WaveToyX: Correct RHS handling after ODESolver update

* StaggeredWaveToyX: Correct RHS handling after ODESolver update

* CI: Update Docker images to AMReX 24.11

* FluxWaveToyX: Correct schedule and test output

* TestODESolvers: Correct schedule

* TestProlongate: Remove unsupporte 1st-order ENO operators

* CarpetX: only use BLOSC2 if compiled into ADIOS2 (EinsteinToolkit#314)

* CarpetX: Reduce termination condition

Reduce the final termination condition, not just the run time. This guarantees that all processes will terminate at the same iteration.

* CarpetX: Correct error check in interpolator

* CarpetX: Simplify error output in interpolator

* CarpetX: Slight performance improvement in interpolator

* CarpetX: Improve safety/reduce warnings in interpolator

* TestInterpolate: Test cell-centred interpolation

* Loop: explicitly create function args to avoid compiler issue

* ODESolvers: rename CarpetX to Driver

* Subcycling: rename CarpetX to Driver

* TestDerivs: remove Arith:: for cos

* TestSubcyclingMC: rename CarpetX to Driver

* TestSubcyclingMC2: remove INHERITS: CarpetX

* TestDerivs: update derivs-8th-order.par

* TestSubcycling: change implementation to Driver

* TestSubcycling: change implementation to Driver

* TestSubcyclingMC2: change implementation to Driver

---------

Co-authored-by: Erik Schnetter <[email protected]>
Co-authored-by: Steven R. Brandt <[email protected]>
Co-authored-by: Roland Haas <[email protected]>

ADMBaseX/schedule.ccl

@@ -1,41 +1,61 @@
 # Schedule definitions for thorn ADMBaseX
 
-SCHEDULE GROUP ADMBaseX_InitialData AT initial
-{
-} "Schedule group for calculating ADM initial data"
+if (CCTK_IsThornActive("ODESolvers")) {
+
+  SCHEDULE GROUP ADMBaseX_InitialData IN ODESolvers_Initial
+  {
+  } "Schedule group for calculating ADM initial data"
 
-SCHEDULE GROUP ADMBaseX_InitialGauge AT initial AFTER ADMBaseX_InitialData
-{
-} "Schedule group for the ADM initial gauge condition"
+  SCHEDULE GROUP ADMBaseX_InitialGauge IN ODESolvers_Initial AFTER ADMBaseX_InitialData
+  {
+  } "Schedule group for the ADM initial gauge condition"
 
-SCHEDULE GROUP ADMBaseX_PostInitial AT initial AFTER (ADMBaseX_InitialData ADMBaseX_InitialGauge)
-{
-} "Schedule group for modifying the ADM initial data, such as e.g. adding noise"
+  SCHEDULE GROUP ADMBaseX_PostInitial IN ODESolvers_Initial AFTER (ADMBaseX_InitialData ADMBaseX_InitialGauge)
+  {
+  } "Schedule group for modifying the ADM initial data, such as e.g. adding noise"
 
-if (CCTK_IsThornActive("ODESolvers")) {
   SCHEDULE GROUP ADMBaseX_SetADMVars IN ODESolvers_PostStep
   {
   } "Set ADM variables in this group"
 
   SCHEDULE GROUP ADMBaseX_SetADMRHS IN ODESolvers_PostStep
   {
   } "Set ADM RHS variables in this group"  
+
 } else {
+
+  SCHEDULE GROUP ADMBaseX_SetADMVars AT post_recover_variables
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP ADMBaseX_SetADMRHS AT post_recover_variables
+  {
+  } "Set ADM RHS variables in this group"
+
   SCHEDULE GROUP ADMBaseX_SetADMVars AT postregrid
   {
   } "Set ADM variables in this group"
 
-  SCHEDULE GROUP ADMBaseX_SetADMVars AT poststep
+  SCHEDULE GROUP ADMBaseX_SetADMRHS AT postregrid
+  {
+  } "Set ADM RHS variables in this group"
+
+  SCHEDULE GROUP ADMBaseX_SetADMVars AT postrestrict
   {
   } "Set ADM variables in this group"
 
-  SCHEDULE GROUP ADMBaseX_SetADMRHS AT postregrid
+  SCHEDULE GROUP ADMBaseX_SetADMRHS AT postrestrict
   {
   } "Set ADM RHS variables in this group"  
+
+  SCHEDULE GROUP ADMBaseX_SetADMVars AT poststep
+  {
+  } "Set ADM variables in this group"
 
   SCHEDULE GROUP ADMBaseX_SetADMRHS AT poststep
   {
   } "Set ADM RHS variables in this group"  
+
 }
 
 if (CCTK_EQUALS(initial_data, "Cartesian Minkowski")) {

Algo/test/test_roots.par

@@ -6,3 +6,9 @@ ActiveThorns = "
 
 Cactus::presync_mode = "mixed-error"
 Cactus::cctk_itlast = 0
+
+IO::out_fileinfo = "axis labels"
+IO::parfile_write = no
+
+CarpetX::out_metadata = no
+

Arith/src/defs.hxx

@@ -138,14 +138,16 @@ template <> struct nan<double> {
 
 // An explicitly unrolled for loop
 
-template <int imin, int imax, int istep = 1, typename F,
-          enable_if_t<(istep > 0 ? imin >= imax : imin <= imax)> * = nullptr>
+template <
+    int imin, int imax, int istep = 1, typename F,
+    std::enable_if_t<(istep > 0 ? imin >= imax : imin <= imax)> * = nullptr>
 constexpr ARITH_INLINE ARITH_DEVICE ARITH_HOST void unroll_for(const F &f) {
   // done: do nothing
 }
 
-template <int imin, int imax, int istep = 1, typename F,
-          enable_if_t<!(istep > 0 ? imin >= imax : imin <= imax)> * = nullptr>
+template <
+    int imin, int imax, int istep = 1, typename F,
+    std::enable_if_t<!(istep > 0 ? imin >= imax : imin <= imax)> * = nullptr>
 constexpr ARITH_INLINE ARITH_DEVICE ARITH_HOST void unroll_for(const F &f) {
   f(imin);
   unroll_for<imin + istep, imax, istep>(f);
@@ -206,6 +208,12 @@ constexpr ARITH_INLINE ARITH_DEVICE ARITH_HOST T flipsign(const T &x,
   return copysign(T(1), y) * x;
 }
 
+// Return 1/x
+template <typename T>
+constexpr ARITH_INLINE ARITH_DEVICE ARITH_HOST T inv(const T &x) {
+  return T(1) / x;
+}
+
 // A max function that returns nan when any argument is nan
 template <typename T>
 constexpr ARITH_INLINE ARITH_DEVICE ARITH_HOST T max1(const T &x, const T &y) {

Arith/src/simd.cxx

@@ -3,6 +3,7 @@
 #include <cctk.h>
 
 #include <cassert>
+#include <cmath>
 
 namespace Arith {
 
@@ -25,17 +26,34 @@ std::size_t get_memop_count() {
   return count;
 }
 
+template <typename T> constexpr bool isequal(const simd<T> x, const T a) {
+  return all(x == a);
+}
+template <typename T> constexpr bool isequal(const simdl<T> x, const bool a) {
+  return all(x == a);
+}
+template <typename T> constexpr bool isapprox(const simd<T> x, const T a) {
+  return all(fabs(x - a) < 1.0e-14);
+}
+
+void check(const bool isgood) {
+  if (isgood)
+    return;
+  CCTK_ERROR("Test failure");
+}
+
 void TestSIMD() {
   // nvcc V11.1.74 doesn't accept this as "constexpr" values
 #ifndef __CUDACC__
-  typedef simd<CCTK_REAL> realv;
+  using real = CCTK_REAL;
+  using realv = simd<real>;
 
   realv x;
   realv y = 0;
   realv z = zero<realv>();
 
-  assert(all(y == 0));
-  assert(all(y == z));
+  check(all(y == 0));
+  check(all(y == z));
 
   realv a = 2;
   realv b = 3;
@@ -45,10 +63,131 @@ void TestSIMD() {
   realv r1 = mulsub(a, b, c);
   realv r2 = negmuladd(a, b, c);
   realv r3 = negmulsub(a, b, c);
-  assert(all(r0 == muladd(2, 3, 4)));
-  assert(all(r1 == mulsub(2, 3, 4)));
-  assert(all(r2 == negmuladd(2, 3, 4)));
-  assert(all(r3 == negmulsub(2, 3, 4)));
+  check(all(r0 == muladd(2, 3, 4)));
+  check(all(r1 == mulsub(2, 3, 4)));
+  check(all(r2 == negmuladd(2, 3, 4)));
+  check(all(r3 == negmulsub(2, 3, 4)));
+
+  real s = 2;
+  real t = 2;
+  real u = 4;
+
+  check(isequal(+a, +s));
+  check(isequal(-a, -s));
+
+  check(isequal(a + b, s + t));
+  check(isequal(a - b, s - t));
+  check(isequal(a * b, s * t));
+  check(isequal(a / b, s / t));
+  check(isequal(s + b, s + t));
+  check(isequal(s - b, s - t));
+  check(isequal(s * b, s * t));
+  check(isequal(s / b, s / t));
+  check(isequal(a + t, s + t));
+  check(isequal(a - t, s - t));
+  check(isequal(a * t, s * t));
+  check(isequal(a / t, s / t));
+
+  check(isequal(a == b, s == t));
+  check(isequal(a != b, s != t));
+  check(isequal(a < b, s < t));
+  check(isequal(a > b, s > t));
+  check(isequal(a <= b, s <= t));
+  check(isequal(a >= b, s >= t));
+  check(isequal(s == b, s == t));
+  check(isequal(s != b, s != t));
+  check(isequal(s < b, s < t));
+  check(isequal(s > b, s > t));
+  check(isequal(s <= b, s <= t));
+  check(isequal(s >= b, s >= t));
+  check(isequal(a == t, s == t));
+  check(isequal(a != t, s != t));
+  check(isequal(a < t, s < t));
+  check(isequal(a > t, s > t));
+  check(isequal(a <= t, s <= t));
+  check(isequal(a >= t, s >= t));
+
+  check(isapprox(abs(a), abs(s)));
+  check(isapprox(acos(a), acos(s)));
+  check(isapprox(acosh(a), acosh(s)));
+  check(allisfinite(a) == allisfinite(s));
+  check(anyisnan(a) == anyisnan(s));
+  check(isapprox(asin(a), asin(s)));
+  check(isapprox(asinh(a), asinh(s)));
+  check(isapprox(atan(a), atan(s)));
+  check(isapprox(atanh(a), atanh(s)));
+  check(isapprox(cbrt(a), cbrt(s)));
+  // check(isapprox(cis(a), cis(s)));
+  // check(isapprox(cispi(a), cispi(s)));
+  check(isapprox(copysign(a, b), copysign(s, t)));
+  check(isapprox(cos(a), cos(s)));
+  check(isapprox(cosh(a), cosh(s)));
+  check(isapprox(cospi(a), cos(CCTK_REAL(M_PI) * s)));
+  check(isapprox(exp(a), exp(s)));
+  check(isapprox(exp10(a), pow(CCTK_REAL(10), s)));
+  check(isapprox(exp2(a), exp2(s)));
+  check(isapprox(fabs(a), fabs(s)));
+  check(isapprox(flipsign(a, b), flipsign(s, t)));
+  check(isapprox(fmax(a, b), fmax(s, t)));
+  check(isapprox(fmax(a, t), fmax(s, t)));
+  check(isapprox(fmax(s, b), fmax(s, t)));
+  check(isapprox(fmin(a, b), fmin(s, t)));
+  check(isapprox(fmin(a, t), fmin(s, t)));
+  check(isapprox(fmin(s, b), fmin(s, t)));
+  check(isapprox(hypot(a, b), hypot(s, t)));
+  check(isapprox(inv(a), inv(s)));
+  check(isequal(isfinite(a), isfinite(s)));
+  check(isequal(isinf(a), isinf(s)));
+  check(isequal(isnan(a), isnan(s)));
+  check(isapprox(log(a), log(s)));
+  check(isapprox(log10(a), log10(s)));
+  check(isapprox(log2(a), log2(s)));
+  check(isapprox(max(a, b), max(s, t)));
+  check(isapprox(max(a, t), max(s, t)));
+  check(isapprox(max(s, b), max(s, t)));
+  check(isapprox(min(a, b), min(s, t)));
+  check(isapprox(min(a, t), min(s, t)));
+  check(isapprox(min(s, b), min(s, t)));
+  check(isapprox(muladd(a, b, c), muladd(s, t, u)));
+  check(isapprox(muladd(a, b, u), muladd(s, t, u)));
+  check(isapprox(muladd(a, t, c), muladd(s, t, u)));
+  check(isapprox(muladd(a, t, u), muladd(s, t, u)));
+  check(isapprox(muladd(s, b, c), muladd(s, t, u)));
+  check(isapprox(muladd(s, b, u), muladd(s, t, u)));
+  check(isapprox(muladd(s, t, c), muladd(s, t, u)));
+  check(isapprox(mulsub(a, b, c), mulsub(s, t, u)));
+  check(isapprox(mulsub(a, b, u), mulsub(s, t, u)));
+  check(isapprox(mulsub(a, t, c), mulsub(s, t, u)));
+  check(isapprox(mulsub(a, t, u), mulsub(s, t, u)));
+  check(isapprox(mulsub(s, b, c), mulsub(s, t, u)));
+  check(isapprox(mulsub(s, b, u), mulsub(s, t, u)));
+  check(isapprox(mulsub(s, t, c), mulsub(s, t, u)));
+  check(isapprox(negmuladd(a, b, c), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(a, b, u), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(a, t, c), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(a, t, u), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(s, b, c), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(s, b, u), negmuladd(s, t, u)));
+  check(isapprox(negmuladd(s, t, c), negmuladd(s, t, u)));
+  check(isapprox(negmulsub(a, b, c), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(a, b, u), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(a, t, c), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(a, t, u), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(s, b, c), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(s, b, u), negmulsub(s, t, u)));
+  check(isapprox(negmulsub(s, t, c), negmulsub(s, t, u)));
+  check(isapprox(pow(a, b), pow(s, t)));
+  check(isapprox(pow(a, t), pow(s, t)));
+  check(isapprox(pow(s, b), pow(s, t)));
+  check(isequal(signbit(a), signbit(s)));
+  check(isapprox(sin(a), sin(s)));
+  check(isapprox(sinh(a), sinh(s)));
+  check(isapprox(sinpi(a), sin(CCTK_REAL(M_PI) * s)));
+  check(isapprox(sqrt(a), sqrt(s)));
+  check(isapprox(tan(a), tan(s)));
+  check(isapprox(tanh(a), tanh(s)));
+  // check(isapprox(tanpi(a), tanpi(s)));
+
 #endif
 }