Merge branch 'develop' into BenWibking/enable-asan-ci

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1099]](https://github.com/parthenon-hpc-lab/parthenon/pull/1099) Functionality for outputting task graphs in GraphViz format.
 - [[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 1106]](https://github.com/parthenon-hpc-lab/parthenon/pull/1106) Add CMake options for turning on ASAN and HWASAN
@@ -27,6 +28,7 @@
 - [[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 1111]](https://github.com/parthenon-hpc-lab/parthenon/pull/1111) Fix undefined behavior due to bitshift of negative number in LogicalLocation
 - [[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
@@ -51,6 +53,7 @@
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR 1114]](https://github.com/parthenon-hpc-lab/parthenon/pull/1114) Enable sanitizers for extended CI host build
+- [[PR 1116]](https://github.com/parthenon-hpc-lab/parthenon/pull/1116) Fix NumPy 2.0 test script breakage
 - [[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

example/poisson_gmg/poisson_driver.cpp

@@ -87,15 +87,15 @@ TaskCollection PoissonDriver::MakeTaskCollection(BlockList_t &blocks) {
     // known when we solve A.u = rhs
     auto get_rhs = none;
     if (use_exact_rhs) {
-      auto copy_exact = tl.AddTask(get_rhs, solvers::utils::CopyData<exact, u>, md);
+      auto copy_exact = tl.AddTask(get_rhs, TF(solvers::utils::CopyData<exact, u>), md);
       auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(copy_exact, tl, md, true);
       PoissonEquation eqs;
       eqs.do_flux_cor = flux_correct;
       get_rhs = eqs.Ax<u, rhs>(tl, comm, md);
     }
 
     // Set initial solution guess to zero
-    auto zero_u = tl.AddTask(get_rhs, solvers::utils::SetToZero<u>, md);
+    auto zero_u = tl.AddTask(get_rhs, TF(solvers::utils::SetToZero<u>), md);
 
     auto solve = zero_u;
     if (solver == "BiCGSTAB") {
@@ -111,8 +111,8 @@ TaskCollection PoissonDriver::MakeTaskCollection(BlockList_t &blocks) {
     // If we are using a rhs to which we know the exact solution, compare our computed
     // solution to the exact solution
     if (use_exact_rhs) {
-      auto diff = tl.AddTask(solve, solvers::utils::AddFieldsAndStore<exact, u, u>, md,
-                             1.0, -1.0);
+      auto diff = tl.AddTask(solve, TF(solvers::utils::AddFieldsAndStore<exact, u, u>),
+                             md, 1.0, -1.0);
       auto get_err = solvers::utils::DotProduct<u, u>(diff, tl, &err, md);
       tl.AddTask(
           get_err,
@@ -127,7 +127,6 @@ TaskCollection PoissonDriver::MakeTaskCollection(BlockList_t &blocks) {
           this, i);
     }
   }
-
   return tc;
 }
 

example/sparse_advection/sparse_advection_driver.cpp

@@ -88,7 +88,8 @@ TaskCollection SparseAdvectionDriver::MakeTaskCollection(BlockList_t &blocks,
     // effectively, sc1 = sc0 + dudt*dt
     auto &sc1 = pmb->meshblock_data.Get(stage_name[stage]);
 
-    auto advect_flux = tl.AddTask(none, sparse_advection_package::CalculateFluxes, sc0);
+    auto advect_flux =
+        tl.AddTask(none, TF(sparse_advection_package::CalculateFluxes), sc0);
   }
 
   const int num_partitions = pmesh->DefaultNumPartitions();
@@ -103,29 +104,29 @@ TaskCollection SparseAdvectionDriver::MakeTaskCollection(BlockList_t &blocks,
     auto &mdudt = pmesh->mesh_data.GetOrAdd("dUdt", i);
 
     const auto any = parthenon::BoundaryType::any;
-    auto start_flxcor = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mc0);
-    auto start_bound = tl.AddTask(none, parthenon::StartReceiveBoundBufs<any>, mc1);
+    auto start_flxcor = tl.AddTask(none, TF(parthenon::StartReceiveFluxCorrections), mc0);
+    auto start_bound = tl.AddTask(none, TF(parthenon::StartReceiveBoundBufs<any>), mc1);
 
     auto set_flxcor =
         parthenon::AddFluxCorrectionTasks(start_flxcor, tl, mc0, pmesh->multilevel);
 
     // compute the divergence of fluxes of conserved variables
-    auto flux_div =
-        tl.AddTask(set_flxcor, FluxDivergence<MeshData<Real>>, mc0.get(), mdudt.get());
+    auto flux_div = tl.AddTask(set_flxcor, TF(FluxDivergence<MeshData<Real>>), mc0.get(),
+                               mdudt.get());
 
-    auto avg_data = tl.AddTask(flux_div, AverageIndependentData<MeshData<Real>>,
+    auto avg_data = tl.AddTask(flux_div, TF(AverageIndependentData<MeshData<Real>>),
                                mc0.get(), mbase.get(), beta);
     // apply du/dt to all independent fields in the container
-    auto update = tl.AddTask(avg_data, UpdateIndependentData<MeshData<Real>>, mc0.get(),
-                             mdudt.get(), beta * dt, mc1.get());
+    auto update = tl.AddTask(avg_data, TF(UpdateIndependentData<MeshData<Real>>),
+                             mc0.get(), mdudt.get(), beta * dt, mc1.get());
 
     // do boundary exchange
     auto boundary =
         parthenon::AddBoundaryExchangeTasks(update, tl, mc1, pmesh->multilevel);
 
     // if this is the last stage, check if we can deallocate any sparse variables
     if (stage == integrator->nstages) {
-      tl.AddTask(boundary, SparseDealloc, mc1.get());
+      tl.AddTask(boundary, TF(SparseDealloc), mc1.get());
     }
   }
 
@@ -138,20 +139,20 @@ TaskCollection SparseAdvectionDriver::MakeTaskCollection(BlockList_t &blocks,
     auto &sc1 = pmb->meshblock_data.Get(stage_name[stage]);
 
     // set physical boundaries
-    auto set_bc = tl.AddTask(none, parthenon::ApplyBoundaryConditions, sc1);
+    auto set_bc = tl.AddTask(none, TF(parthenon::ApplyBoundaryConditions), sc1);
 
     // estimate next time step
     if (stage == integrator->nstages) {
-      auto new_dt = tl.AddTask(set_bc, EstimateTimestep<MeshBlockData<Real>>, sc1.get());
+      auto new_dt =
+          tl.AddTask(set_bc, TF(EstimateTimestep<MeshBlockData<Real>>), sc1.get());
 
       // Update refinement
       if (pmesh->adaptive) {
         auto tag_refine = tl.AddTask(
-            set_bc, parthenon::Refinement::Tag<MeshBlockData<Real>>, sc1.get());
+            set_bc, TF(parthenon::Refinement::Tag<MeshBlockData<Real>>), sc1.get());
       }
     }
   }
-
   return tc;
 }
 

scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py

@@ -648,11 +648,11 @@ def Get(self, variable, flatten=True, interior=False, average_to_cell_centers=Tr
         if flatten:
             nblocks = vShape[0]
             if self.varTopology[variable] == "None":
-                remaining_size = np.product(vShape[1:])
+                remaining_size = np.prod(vShape[1:])
                 return self.varData[variable].reshape(nblocks, remaining_size)
             else:
                 preserved_shape = vShape[:-3]
-                remaining_size = np.product(vShape[-3:])
+                remaining_size = np.prod(vShape[-3:])
                 return self.varData[variable].reshape(*preserved_shape, remaining_size)
 
         if self.IncludesGhost and interior:

src/bvals/comms/boundary_communication.cpp

@@ -413,16 +413,16 @@ TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl,
 
   // auto out = (pro_local | pro);
 
-  auto send = tl.AddTask(dependency, SendBoundBufs<bounds>, md);
-  auto recv = tl.AddTask(dependency, ReceiveBoundBufs<bounds>, md);
-  auto set = tl.AddTask(recv, SetBounds<bounds>, md);
+  auto send = tl.AddTask(dependency, TF(SendBoundBufs<bounds>), md);
+  auto recv = tl.AddTask(dependency, TF(ReceiveBoundBufs<bounds>), md);
+  auto set = tl.AddTask(recv, TF(SetBounds<bounds>), md);
 
   auto pro = set;
   if (md->GetMeshPointer()->multilevel) {
-    auto cbound = tl.AddTask(set, ApplyBoundaryConditionsOnCoarseOrFineMD, md, true);
-    pro = tl.AddTask(cbound, ProlongateBounds<bounds>, md);
+    auto cbound = tl.AddTask(set, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, true);
+    pro = tl.AddTask(cbound, TF(ProlongateBounds<bounds>), md);
   }
-  auto fbound = tl.AddTask(pro, ApplyBoundaryConditionsOnCoarseOrFineMD, md, false);
+  auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false);
 
   return fbound;
 }
@@ -437,8 +437,9 @@ AddBoundaryExchangeTasks<BoundaryType::gmg_same>(TaskID, TaskList &,
 TaskID AddFluxCorrectionTasks(TaskID dependency, TaskList &tl,
                               std::shared_ptr<MeshData<Real>> &md, bool multilevel) {
   if (!multilevel) return dependency;
-  tl.AddTask(dependency, SendBoundBufs<BoundaryType::flxcor_send>, md);
-  auto receive = tl.AddTask(dependency, ReceiveBoundBufs<BoundaryType::flxcor_recv>, md);
-  return tl.AddTask(receive, SetBounds<BoundaryType::flxcor_recv>, md);
+  tl.AddTask(dependency, TF(SendBoundBufs<BoundaryType::flxcor_send>), md);
+  auto receive =
+      tl.AddTask(dependency, TF(ReceiveBoundBufs<BoundaryType::flxcor_recv>), md);
+  return tl.AddTask(receive, TF(SetBounds<BoundaryType::flxcor_recv>), md);
 }
 } // namespace parthenon

src/mesh/forest/logical_location.cpp

@@ -118,10 +118,10 @@ bool LogicalLocation::IsNeighbor(const LogicalLocation &in) const {
 
   bool neighbors = true;
   for (int dir = 0; dir < 3; ++dir) {
-    auto low = (l(dir) << level_shift_this) - 1;
+    auto low = l(dir) * block_size_this - 1;
     auto hi = low + block_size_this + 1;
 
-    auto low_in = (in.l(dir) << level_shift_in);
+    auto low_in = in.l(dir) * block_size_in;
     auto hi_in = low_in + block_size_in - 1;
     neighbors = neighbors && !(hi < low_in || low > hi_in);
   }
@@ -140,7 +140,7 @@ bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in,
 
   bool neighbors = true;
   for (int dir = 0; dir < 3; ++dir) {
-    auto low = (l(dir) << level_shift_this);
+    auto low = l(dir) * block_size_this;
     auto hi = low + block_size_this - 1;
     if (te_offset[dir] == -1) {
       low -= 1;
@@ -150,7 +150,7 @@ bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in,
       low = hi - 1;
     }
 
-    auto low_in = (in.l(dir) << level_shift_in);
+    auto low_in = in.l(dir) * block_size_in;
     auto hi_in = low_in + block_size_in - 1;
     neighbors = neighbors && !(hi < low_in || low > hi_in);
   }

src/solvers/bicgstab_solver.hpp

@@ -103,15 +103,16 @@ class BiCGSTABSolver {
 
     // Initialization: x <- 0, r <- rhs, rhat0 <- rhs,
     // rhat0r_old <- (rhat0, r), p <- r, u <- 0
-    auto zero_x = tl.AddTask(dependence, SetToZero<x>, md);
-    auto zero_u_init = tl.AddTask(dependence, SetToZero<u>, md);
-    auto copy_r = tl.AddTask(dependence, CopyData<rhs, r>, md);
-    auto copy_p = tl.AddTask(dependence, CopyData<rhs, p>, md);
-    auto copy_rhat0 = tl.AddTask(dependence, CopyData<rhs, rhat0>, md);
+    auto zero_x = tl.AddTask(dependence, TF(SetToZero<x>), md);
+    auto zero_u_init = tl.AddTask(dependence, TF(SetToZero<u>), md);
+    auto copy_r = tl.AddTask(dependence, TF(CopyData<rhs, r>), md);
+    auto copy_p = tl.AddTask(dependence, TF(CopyData<rhs, p>), md);
+    auto copy_rhat0 = tl.AddTask(dependence, TF(CopyData<rhs, rhat0>), md);
     auto get_rhat0r_init = DotProduct<rhat0, r>(dependence, tl, &rhat0r, md);
     auto initialize = tl.AddTask(
         TaskQualifier::once_per_region | TaskQualifier::local_sync,
         zero_x | zero_u_init | copy_r | copy_p | copy_rhat0 | get_rhat0r_init,
+        "zero factors",
         [](BiCGSTABSolver *solver) {
           solver->rhat0r_old = solver->rhat0r.val;
           solver->rhat0r.val = 0.0;
@@ -123,7 +124,7 @@ class BiCGSTABSolver {
           return TaskStatus::complete;
         },
         this);
-    tl.AddTask(TaskQualifier::once_per_region, dependence, [&]() {
+    tl.AddTask(TaskQualifier::once_per_region, dependence, "print to screen", [&]() {
       if (Globals::my_rank == 0 && params_.print_per_step)
         printf("# [0] v-cycle\n# [1] rms-residual\n# [2] rms-error\n");
       return TaskStatus::complete;
@@ -135,12 +136,12 @@ class BiCGSTABSolver {
     // 1. u <- M p
     auto precon1 = none;
     if (params_.precondition) {
-      auto set_rhs = itl.AddTask(precon1, CopyData<p, rhs>, md);
-      auto zero_u = itl.AddTask(precon1, SetToZero<u>, md);
+      auto set_rhs = itl.AddTask(precon1, TF(CopyData<p, rhs>), md);
+      auto zero_u = itl.AddTask(precon1, TF(SetToZero<u>), md);
       precon1 =
           preconditioner.AddLinearOperatorTasks(itl, set_rhs | zero_u, partition, pmesh);
     } else {
-      precon1 = itl.AddTask(none, CopyData<p, u>, md);
+      precon1 = itl.AddTask(none, TF(CopyData<p, u>), md);
     }
 
     // 2. v <- A u
@@ -153,7 +154,7 @@ class BiCGSTABSolver {
 
     // 4. h <- x + alpha u (alpha = rhat0r_old / rhat0v)
     auto correct_h = itl.AddTask(
-        get_rhat0v,
+        get_rhat0v, "h <- x + alpha u",
         [](BiCGSTABSolver *solver, std::shared_ptr<MeshData<Real>> &md) {
           Real alpha = solver->rhat0r_old / solver->rhat0v.val;
           return AddFieldsAndStore<x, u, h>(md, 1.0, alpha);
@@ -162,7 +163,7 @@ class BiCGSTABSolver {
 
     // 5. s <- r - alpha v (alpha = rhat0r_old / rhat0v)
     auto correct_s = itl.AddTask(
-        get_rhat0v,
+        get_rhat0v, "s <- r - alpha v",
         [](BiCGSTABSolver *solver, std::shared_ptr<MeshData<Real>> &md) {
           Real alpha = solver->rhat0r_old / solver->rhat0v.val;
           return AddFieldsAndStore<r, v, s>(md, 1.0, -alpha);
@@ -185,12 +186,12 @@ class BiCGSTABSolver {
     // 6. u <- M s
     auto precon2 = correct_s;
     if (params_.precondition) {
-      auto set_rhs = itl.AddTask(precon2, CopyData<s, rhs>, md);
-      auto zero_u = itl.AddTask(precon2, SetToZero<u>, md);
+      auto set_rhs = itl.AddTask(precon2, TF(CopyData<s, rhs>), md);
+      auto zero_u = itl.AddTask(precon2, TF(SetToZero<u>), md);
       precon2 =
           preconditioner.AddLinearOperatorTasks(itl, set_rhs | zero_u, partition, pmesh);
     } else {
-      precon2 = itl.AddTask(precon2, CopyData<s, u>, md);
+      precon2 = itl.AddTask(precon2, TF(CopyData<s, u>), md);
     }
 
     // 7. t <- A u
@@ -204,7 +205,7 @@ class BiCGSTABSolver {
 
     // 9. x <- h + omega u
     auto correct_x = itl.AddTask(
-        TaskQualifier::local_sync, get_tt | get_ts,
+        TaskQualifier::local_sync, get_tt | get_ts, "x <- h + omega u",
         [](BiCGSTABSolver *solver, std::shared_ptr<MeshData<Real>> &md) {
           Real omega = solver->ts.val / solver->tt.val;
           return AddFieldsAndStore<h, u, x>(md, 1.0, omega);
@@ -213,7 +214,7 @@ class BiCGSTABSolver {
 
     // 10. r <- s - omega t
     auto correct_r = itl.AddTask(
-        get_tt | get_ts,
+        get_tt | get_ts, "r <- s - omega t",
         [](BiCGSTABSolver *solver, std::shared_ptr<MeshData<Real>> &md) {
           Real omega = solver->ts.val / solver->tt.val;
           return AddFieldsAndStore<s, t, r>(md, 1.0, -omega);
@@ -240,6 +241,7 @@ class BiCGSTABSolver {
     // 13. p <- r + beta * (p - omega * v)
     auto update_p = itl.AddTask(
         TaskQualifier::local_sync, get_rhat0r | get_res2,
+        "p <- r + beta * (p - omega * v)",
         [](BiCGSTABSolver *solver, std::shared_ptr<MeshData<Real>> &md) {
           Real alpha = solver->rhat0r_old / solver->rhat0v.val;
           Real omega = solver->ts.val / solver->tt.val;
@@ -256,7 +258,7 @@ class BiCGSTABSolver {
     auto check = itl.AddTask(
         TaskQualifier::completion | TaskQualifier::once_per_region |
             TaskQualifier::global_sync,
-        update_p | correct_x,
+        update_p | correct_x, "rhat0r_old <- rhat0r",
         [](BiCGSTABSolver *solver, Mesh *pmesh, int max_iter, Real *res_tol) {
           solver->iter_counter++;
           Real rms_res = std::sqrt(solver->residual.val / pmesh->GetTotalCells());
@@ -277,7 +279,7 @@ class BiCGSTABSolver {
         },
         this, pmesh, params_.max_iters, ptol);
 
-    return tl.AddTask(solver_id, CopyData<x, u>, md);
+    return tl.AddTask(solver_id, TF(CopyData<x, u>), md);
   }
 
   Real GetSquaredResidualSum() const { return residual.val; }