parthenon-hpc-lab · BenWibking · Feb 21, 2024 · Aug 29, 2023 · Aug 29, 2023 · Sep 26, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -0,0 +1,23 @@
+# Changelog
+
+## Current develop (i.e., `main` branch)
+
+### Added (new features/APIs/variables/...)
+
+### Changed (changing behavior/API/variables/...)
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
+
+### Fixed (not changing behavior/API/variables/...)
+
+### Infrastructure
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Added `CHANGELOG.md`
+
+### Removed (removing behavior/API/varaibles/...)
+
+### Incompatibilities (i.e. breaking changes)
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
+  - Updated access to block dimension: `pmb->block_size.nx1` -> `pmb->block_size.nx(X1DIR)` (and similarly x2 and x3)
+  - Update access to mesh size: `pmesh->mesh_size.x1max` -> `pmesh->mesh_size.xmax(X1DIR)` (and similarly x2, x3, and min)
+  - Updated Parthenon `GradMinMod` signature for custom prolongation ops
+  - `GetBlockPointer` returns a raw pointer not a shared one (and updated interfaces to use raw pointers rather than shared ones)
+
diff --git a/external/parthenon b/external/parthenon
diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp
@@ -775,8 +775,8 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
   int il, iu, jl, ju, kl, ku;
   jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e;
   // TODO(pgrete): are these looop limits are likely too large for 2nd order
-  if (pmb->block_size.nx2 > 1) {
-    if (pmb->block_size.nx3 == 1) // 2D
+  if (pmb->block_size.nx(X2DIR) > 1) {
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e;
     else // 3D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1;
@@ -848,7 +848,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
         parthenon::ScratchPad2D<Real>::shmem_size(num_scratch_vars, nx1) * 3;
     // set the loop limits
     il = ib.s - 1, iu = ib.e + 1, kl = kb.s, ku = kb.e;
-    if (pmb->block_size.nx3 == 1) // 2D
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       kl = kb.s, ku = kb.e;
     else // 3D
       kl = kb.s - 1, ku = kb.e + 1;

diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp
@@ -81,34 +81,40 @@ struct ProlongateCellMinModMultiD {
 
     Real dx1fm = 0;
     Real dx1fp = 0;
+    [[maybe_unused]] Real gx1m = 0, gx1p = 0;
     Real gx1c = 0;
     if constexpr (INCLUDE_X1) {
       Real dx1m, dx1p;
       GetGridSpacings<1, el>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, &dx1fm,
                              &dx1fp);
-      gx1c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1),
-                        coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p);
+      gx1c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1),
+                     coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p, gx1m, gx1p);
     }
 
     Real dx2fm = 0;
     Real dx2fp = 0;
+    [[maybe_unused]] Real gx2m = 0, gx2p = 0;
     Real gx2c = 0;
     if constexpr (INCLUDE_X2) {
       Real dx2m, dx2p;
       GetGridSpacings<2, el>(coords, coarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, &dx2fm,
                              &dx2fp);
-      gx2c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i),
-                        coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p);
+      gx2c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i),
+                     coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p, gx2m, gx2p);
     }
     Real dx3fm = 0;
     Real dx3fp = 0;
+    [[maybe_unused]] Real gx3m = 0, gx3p = 0;
     Real gx3c = 0;
     if constexpr (INCLUDE_X3) {
       Real dx3m, dx3p;
       GetGridSpacings<3, el>(coords, coarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, &dx3fm,
                              &dx3fp);
-      gx3c = GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i),
-                        coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p);
+      gx3c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i),
+                     coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p, gx3m, dx3p);
     }
 
     // Max. expected total difference. (dx#fm/p are positive by construction)

diff --git a/src/pgen/cloud.cpp b/src/pgen/cloud.cpp
@@ -213,7 +213,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
 }
 
 void InflowWindX2(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse) {
-  std::shared_ptr<MeshBlock> pmb = mbd->GetBlockPointer();
+  auto pmb = mbd->GetBlockPointer();
   auto cons = mbd->PackVariables(std::vector<std::string>{"cons"}, coarse);
   // TODO(pgrete) Add par_for_bndry to Parthenon without requiring nb
   const auto nb = IndexRange{0, 0};

diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp
@@ -25,6 +25,7 @@
 #include <string>    // c_str()
 
 // Parthenon headers
+#include "Kokkos_MathematicalFunctions.hpp"
 #include "kokkos_abstraction.hpp"
 #include "mesh/domain.hpp"
 #include "mesh/mesh.hpp"
@@ -347,6 +348,11 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
   hydro_pkg->AddField("mach_sonic", m);
   // temperature
   hydro_pkg->AddField("temperature", m);
+  // radial velocity
+  hydro_pkg->AddField("v_r", m);
+
+  // spherical theta
+  hydro_pkg->AddField("theta_sph", m);
 
   if (hydro_pkg->Param<Cooling>("enable_cooling") == Cooling::tabular) {
     // cooling time
@@ -359,6 +365,9 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
 
     // plasma beta
     hydro_pkg->AddField("plasma_beta", m);
+
+    // plasma beta
+    hydro_pkg->AddField("B_mag", m);
   }
 
   /************************************************************
@@ -563,7 +572,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
        ************************************************************/
       const auto &magnetic_tower = hydro_pkg->Param<MagneticTower>("magnetic_tower");
 
-      magnetic_tower.AddInitialFieldToPotential(pmb.get(), a_kb, a_jb, a_ib, A);
+      magnetic_tower.AddInitialFieldToPotential(pmb, a_kb, a_jb, a_ib, A);
 
       /************************************************************
        * Add dipole magnetic field to the magnetic potential
@@ -666,7 +675,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
     // Init phases on all blocks
     for (int b = 0; b < md->NumBlocks(); b++) {
       auto pmb = md->GetBlockData(b)->GetBlockPointer();
-      few_modes_ft.SetPhases(pmb.get(), pin);
+      few_modes_ft.SetPhases(pmb, pin);
     }
     // As for t_corr in few_modes_ft, the choice for dt is
     // in principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain
@@ -705,9 +714,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
                                       MPI_SUM, MPI_COMM_WORLD));
 #endif // MPI_PARALLEL
 
-    const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min;
-    const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min;
-    const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min;
+    const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
+    const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
+    const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
     auto v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v)));
 
     pmb->par_for(
@@ -734,7 +743,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
     // Init phases on all blocks
     for (int b = 0; b < md->NumBlocks(); b++) {
       auto pmb = md->GetBlockData(b)->GetBlockPointer();
-      few_modes_ft.SetPhases(pmb.get(), pin);
+      few_modes_ft.SetPhases(pmb, pin);
     }
     // As for t_corr in few_modes_ft, the choice for dt is
     // in principle arbitrary because the inital b_hat is 0 and the b_hat_new will contain
@@ -783,9 +792,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
                                       MPI_SUM, MPI_COMM_WORLD));
 #endif // MPI_PARALLEL
 
-    const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min;
-    const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min;
-    const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min;
+    const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
+    const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
+    const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
     auto b_norm = std::sqrt(b2_sum / (Lx * Ly * Lz) / (SQR(sigma_b)));
 
     pmb->par_for(
@@ -818,6 +827,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   auto &entropy = data->Get("entropy").data;
   auto &mach_sonic = data->Get("mach_sonic").data;
   auto &temperature = data->Get("temperature").data;
+  auto &v_r = data->Get("v_r").data;
+  auto &theta_sph = data->Get("theta_sph").data;
 
   // for computing temperature from primitives
   auto units = pkg->Param<Units>("units");
@@ -844,8 +855,12 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
         const Real x = coords.Xc<1>(i);
         const Real y = coords.Xc<2>(j);
         const Real z = coords.Xc<3>(k);
-        const Real r2 = SQR(x) + SQR(y) + SQR(z);
-        log10_radius(k, j, i) = 0.5 * std::log10(r2);
+        const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z));
+        log10_radius(k, j, i) = std::log10(r);
+
+        v_r(k, j, i) = ((v1 * x) + (v2 * y) + (v3 * z)) / r;
+
+        theta_sph(k, j, i) = std::acos(z / r);
 
         // compute entropy
         const Real K = P / std::pow(rho / mbar, gam);
@@ -884,6 +899,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   if (pkg->Param<Fluid>("fluid") == Fluid::glmmhd) {
     auto &plasma_beta = data->Get("plasma_beta").data;
     auto &mach_alfven = data->Get("mach_alfven").data;
+    auto &b_mag = data->Get("B_mag").data;
 
     pmb->par_for(
         "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
@@ -896,6 +912,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
           const Real Bz = prim(IB3, k, j, i);
           const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz));
 
+          b_mag(k, j, i) = Kokkos::sqrt(B2);
+
           // compute Alfven mach number
           const Real v_A = std::sqrt(B2 / rho);
           const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS

diff --git a/src/pgen/cpaw.cpp b/src/pgen/cpaw.cpp
@@ -34,6 +34,8 @@
 
 namespace cpaw {
 using namespace parthenon::driver::prelude;
+using namespace parthenon::package::prelude;
+
 // Parameters which define initial solution -- made global so that they can be shared
 // with functions A1,2,3 which compute vector potentials
 Real den, pres, gm1, b_par, b_perp, v_perp, v_par;
@@ -213,8 +215,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
   }
 
   // write errors
-  std::fprintf(pfile, "%d  %d", mesh->mesh_size.nx1, mesh->mesh_size.nx2);
-  std::fprintf(pfile, "  %d  %d  %e", mesh->mesh_size.nx3, tm.ncycle, rms_err);
+  std::fprintf(pfile, "%d  %d", mesh->mesh_size.nx(X1DIR), mesh->mesh_size.nx(X2DIR));
+  std::fprintf(pfile, "  %d  %d  %e", mesh->mesh_size.nx(X3DIR), tm.ncycle, rms_err);
   std::fprintf(pfile, "  %e  %e  %e  %e", err[IDN], err[IM1], err[IM2], err[IM3]);
   std::fprintf(pfile, "  %e", err[IEN]);
   std::fprintf(pfile, "  %e  %e  %e", err[IB1], err[IB2], err[IB3]);

diff --git a/src/pgen/field_loop.cpp b/src/pgen/field_loop.cpp
@@ -135,13 +135,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   int iprob = pin->GetInteger("problem/field_loop", "iprob");
   Real ang_2, cos_a2(0.0), sin_a2(0.0), lambda(0.0);
 
-  Real x1size = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min;
-  Real x2size = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min;
+  Real x1size =
+      pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR);
+  Real x2size =
+      pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR);
 
   const bool two_d = pmb->pmy_mesh->ndim < 3;
   // for 2D sim set x3size to zero so that v_z is 0 below
   Real x3size =
-      two_d ? 0 : pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min;
+      two_d ? 0
+            : pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR);
 
   // For (iprob=4) -- rotated cylinder in 3D -- set up rotation angle and wavelength
   if (iprob == 4) {

diff --git a/src/pgen/linear_wave.cpp b/src/pgen/linear_wave.cpp
@@ -32,6 +32,7 @@
 
 namespace linear_wave {
 using namespace parthenon::driver::prelude;
+using namespace parthenon::package::prelude;
 
 // TODO(pgrete) temp fix to address removal in Parthenon. Update when merging with MHD
 constexpr int NWAVE = 5;
@@ -280,9 +281,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
   if (parthenon::Globals::my_rank == 0) {
     // normalize errors by number of cells
     const auto mesh_size = mesh->mesh_size;
-    const auto vol = (mesh_size.x1max - mesh_size.x1min) *
-                     (mesh_size.x2max - mesh_size.x2min) *
-                     (mesh_size.x3max - mesh_size.x3min);
+    const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
+                     (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
+                     (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
     for (int i = 0; i < (NHYDRO + NFIELD); ++i)
       l1_err[i] = l1_err[i] / vol;
     // compute rms error
@@ -320,8 +321,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
     }
 
     // write errors
-    std::fprintf(pfile, "%d  %d", mesh_size.nx1, mesh_size.nx2);
-    std::fprintf(pfile, "  %d  %d", mesh_size.nx3, tm.ncycle);
+    std::fprintf(pfile, "%d  %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
+    std::fprintf(pfile, "  %d  %d", mesh_size.nx(X2DIR), tm.ncycle);
     std::fprintf(pfile, "  %e  %e", rms_err, l1_err[IDN]);
     std::fprintf(pfile, "  %e  %e  %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]);
     std::fprintf(pfile, "  %e", l1_err[IEN]);

diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp
@@ -32,6 +32,7 @@
 
 namespace linear_wave_mhd {
 using namespace parthenon::driver::prelude;
+using namespace parthenon::package::prelude;
 
 constexpr int NMHDWAVE = 7;
 // Parameters which define initial solution -- made global so that they can be shared
@@ -300,9 +301,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
   if (parthenon::Globals::my_rank == 0) {
     // normalize errors by number of cells
     const auto mesh_size = mesh->mesh_size;
-    const auto vol = (mesh_size.x1max - mesh_size.x1min) *
-                     (mesh_size.x2max - mesh_size.x2min) *
-                     (mesh_size.x3max - mesh_size.x3min);
+    const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
+                     (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
+                     (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
     for (int i = 0; i < (NGLMMHD); ++i)
       l1_err[i] = l1_err[i] / vol;
     // compute rms error
@@ -342,8 +343,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
     }
 
     // write errors
-    std::fprintf(pfile, "%d  %d", mesh_size.nx1, mesh_size.nx2);
-    std::fprintf(pfile, "  %d  %d", mesh_size.nx3, tm.ncycle);
+    std::fprintf(pfile, "%d  %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
+    std::fprintf(pfile, "  %d  %d", mesh_size.nx(X3DIR), tm.ncycle);
     std::fprintf(pfile, "  %e  %e", rms_err, l1_err[IDN]);
     std::fprintf(pfile, "  %e  %e  %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]);
     std::fprintf(pfile, "  %e", l1_err[IEN]);

diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp
@@ -217,10 +217,10 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
   const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0;
   const auto p0 = pin->GetReal("problem/turbulence", "p0");
   const auto rho0 = pin->GetReal("problem/turbulence", "rho0");
-  const auto x3min = pmesh->mesh_size.x3min;
-  const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min;
-  const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min;
-  const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min;
+  const auto x3min = pmesh->mesh_size.xmin(X3DIR);
+  const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
+  const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
+  const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
   const auto kz = 2.0 * M_PI / Lz;
 
   // already pack data here to get easy access to coords in kernels
@@ -402,9 +402,12 @@ void Perturb(MeshData<Real> *md, const Real dt) {
                                     MPI_SUM, MPI_COMM_WORLD));
 #endif // MPI_PARALLEL
 
-  const auto Lx = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min;
-  const auto Ly = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min;
-  const auto Lz = pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min;
+  const auto Lx =
+      pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR);
+  const auto Ly =
+      pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR);
+  const auto Lz =
+      pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR);
   const auto accel_rms = hydro_pkg->Param<Real>("turbulence/accel_rms");
   auto norm = accel_rms / std::sqrt(sums[0] / (Lx * Ly * Lz));