Merge branch 'pgrete/bump-parth-for-rr' into pgrete/turb-next

CHANGELOG.md

@@ -0,0 +1,29 @@
+# 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 109]](https://github.com/parthenon-hpc-lab/athenapk/pull/109) Bump Parthenon to latest develop (2024-05-29)
+- [[PR 105]](https://github.com/parthenon-hpc-lab/athenapk/pull/105) Bump Parthenon to latest develop (2024-03-13)
+- [[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 109]](https://github.com/parthenon-hpc-lab/athenapk/pull/109) Bump Parthenon to latest develop (2024-05-29)
+  - Changed signature of `UserWorkBeforeOutput` to include `SimTime` as last paramter
+  - Fixes bitwise idential restarts for AMR simulations (the derefinement counter is now included)
+  - Order of operations in flux-correction has changed (expect round-off error differences to previous results for AMR sims)
+- [[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)
+

docs/cluster.md

@@ -76,14 +76,14 @@ $$
     \frac{3 H_0^2}{8 \pi G}.
 $$
 
-Parameters for the HERNQUIST BCG are controlled via:
+Parameters for the `HERNQUIST` BCG are controlled via:
 ```
 <problem/cluster/gravity>
 
 m_bcg_s = 0.001 # in code_mass
 r_bcg_s = 0.004 # in code_length
 ```
-where a HERNQUIST profile adds a gravitational acceleration defined by
+where a `HERNQUIST` profile adds a gravitational acceleration defined by
 
 $$
  g_{BCG}(r) = G \frac{ M_{BCG} }{R_{BCG}^2} \frac{1}{\left( 1 + \frac{r}{R_{BCG}}\right)^2}
@@ -216,7 +216,7 @@ triggering_mode = COLD_GAS # or NONE, BOOSTED_BONDI, BONDI_SCHAYE
 ```
 where `triggering_mode=NONE` will disable AGN triggering. 
 
-With BOOSTED_BONDI accretion, the mass rate of accretion follows
+With `BOOSTED_BONDI` accretion, the mass rate of accretion follows
 
 $$
 \dot{M} = \alpha \frac { 2 \pi G^2 M^2_{SMBH} \hat {\rho} } {
@@ -235,7 +235,7 @@ m_smbh = 1.0e-06 # in code_mass
 accretion_radius = 0.001 # in code_length
 bondi_alpha= 100.0 # unitless
 ```
-With BONDI_SCHAYE accretion, the `$\alpha$` used for BOOSTED_BONDI accretion is modified to depend on the number density following:
+With `BONDI_SCHAYE` accretion, the $\alpha$ used for `BOOSTED_BONDI` accretion is modified to depend on the number density following:
 
 $$
 \alpha =
@@ -252,7 +252,7 @@ bondi_n0= 2.9379989445851786e+72 # in 1/code_length**3
 bondi_beta= 2.0 # unitless
 ```
 
-With both BOOSTED_BONDI and BONDI_SCHAYE accretion, mass is removed from each
+With both `BOOSTED_BONDI` and `BONDI_SCHAYE` accretion, mass is removed from each
 cell within the accretion zone at a mass weighted rate. E.g. the mass in each
 cell within the accretion region changes by
 ```
@@ -264,7 +264,7 @@ unchanged. Thus velocities and temperatures will increase where mass is
 removed.
 
 
-With COLD_GAS accretion, the accretion rate becomes the total mass within the accretion zone equal to or
+With `COLD_GAS` accretion, the accretion rate becomes the total mass within the accretion zone equal to or
 below a defined cold temperature threshold divided by a defined accretion
 timescale. The temperature threshold and accretion timescale are defined by
 ```
@@ -360,13 +360,17 @@ velocity $v_{jet}$ can be set via
 kinetic_jet_temperature = 1e7 # K
 ```
 However, $T_{jet}$ and $v_{jet}$ must be non-negative and fulfill
+
 $$
-v_{jet} = \sqrt{ 2 \left ( \epsilon c^2 - (1 - \epsilon) \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right} \right ) }
+v_{jet} = \sqrt{ 2 \left ( \epsilon c^2 - (1 - \epsilon) \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right) } \right ) }
 $$
+
 to ensure that the sum of rest mass energy, thermal energy, and kinetic energy of the new gas sums to $\dot{M} c^2$. Note that these equations places limits  on $T_{jet}$ and $v_{jet}$, specifically
+
 $$
-v_{jet} \leq c \sqrt{ 2 \epsilon } \qquad \text{and} \qquad \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right} \leq c^2 \frac{ \epsilon}{1 - \epsilon}
+v_{jet} \leq c \sqrt{ 2 \epsilon } \qquad \text{and} \qquad \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right) } \leq c^2 \frac{ \epsilon}{1 - \epsilon}
 $$
+
 If the above equations are not satified then an exception will be thrown at
 initialization. If neither $T_{jet}$ nor $v_{jet}$ are specified, then
 $v_{jet}$ will be computed assuming $T_{jet}=0$ and a warning will be given
@@ -429,19 +433,19 @@ where the injected magnetic field follows
 
 $$
 \begin{align}
-\mathcal{B}_r      &=\mathcal{B}_0 2 \frac{h r}{\ell^2} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
-\mathcal{B}_\theta &=\mathcal{B}_0 \alpha \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right ) } \\\\
-\mathcal{B}_h      &=\mathcal{B}_0 2 \left( 1 - \frac{r^2}{\ell^2} \right ) \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
+\mathcal{B}\_r      &= \mathcal{B}\_0 2 \frac{h r}{\ell^2} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\
+\mathcal{B}\_\theta &= \mathcal{B}\_0 \alpha \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right ) } \\
+\mathcal{B}\_h      &= \mathcal{B}\_0 2 \left( 1 - \frac{r^2}{\ell^2} \right ) \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
 \end{align}
 $$
 
 which has  the corresponding vector potential field
 
 $$
 \begin{align}
-\mathcal{A}_r &= 0 \\\\
-\mathcal{A}_{\theta} &= \mathcal{B}_0 \ell \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
-\mathcal{A}_h &= \mathcal{B}_0 \ell \frac{\alpha}{2}\exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
+\mathcal{A}\_r &= 0 \\
+\mathcal{A}\_{\theta} &= \mathcal{B}\_0 \ell \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\
+\mathcal{A}\_h &= \mathcal{B}\_0 \ell \frac{\alpha}{2}\exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
 \end{align}
 $$
 
@@ -463,7 +467,7 @@ $$
 where $\dot{\rho}_B$ is set to
 
 $$
-\dot{\rho}_B = \frac{3 \pi}{2} \frac{\dot{M} \left ( 1 - \epsilon \right ) f_{magnetic}}{\ell^3}
+\dot{\rho}_B = \frac{3 \pi}{2} \frac{\dot{M} \left ( 1 - \epsilon \right ) f\_\mathrm{magnetic}}{\ell^3}
 $$
 
 so that the total mass injected matches the accreted mass propotioned to magnetic feedback.
@@ -485,13 +489,17 @@ according to their ratio.
 #### Simple field loop (donut) feedback
 
 Magnetic energy is injected according to the following simple potential
+
 $$
-A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right) \qquad \mathrm{for} \qquad h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
 $$
+
 resultig in a magnetic field configuration of
+
 $$
-B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right) \qquad \mathrm{for} \qquad h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
 $$
+
 with all other components being zero.
 
 
@@ -526,7 +534,7 @@ fixed_mass_rate = 1.0
 
 ## SNIA Feedback
 
-Following [Prasad 2020](doi.org/10.1093/mnras/112.2.195), AthenaPK can inject
+Following [Prasad 2020](https://doi.org/10.3847/1538-4357/abc33c), AthenaPK can inject
 mass and energy from type Ia supernovae following the mass profile of the BCG.
 This SNIA feedback can be configured with
 ```
@@ -562,4 +570,4 @@ temperature_threshold = 2e4 # in K
 ```
 
 Note that all parameters need to be specified explicitly for the feedback to work
-(i.e., no hidden default values).
+(i.e., no hidden default values).

src/hydro/prolongation/custom_ops.hpp

@@ -87,8 +87,9 @@ struct ProlongateCellMinModMultiD {
       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, gx1m, gx1p);
+      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;
@@ -99,8 +100,9 @@ struct ProlongateCellMinModMultiD {
       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, gx2m, gx2p);
+      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;
@@ -110,8 +112,9 @@ struct ProlongateCellMinModMultiD {
       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, gx3m, 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)

src/hydro/rsolvers/hydro_hlle.hpp

@@ -50,7 +50,7 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
     Real gm1 = gamma - 1.0;
     Real igm1 = 1.0 / gm1;
     parthenon::par_for_inner(member, il, iu, [&](const int i) {
-      Real wli[(NHYDRO)], wri[(NHYDRO)];
+      Real wli[(NHYDRO)], wri[(NHYDRO)], wroe[(NHYDRO)];
       Real fl[(NHYDRO)], fr[(NHYDRO)], flxi[(NHYDRO)];
       //--- Step 1.  Load L/R states into local variables
       wli[IDN] = wl(IDN, i);
@@ -65,34 +65,37 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
       wri[IV3] = wr(ivz, i);
       wri[IPR] = wr(IPR, i);
 
-      //--- Step 2.  Compute middle state estimates with PVRS (Toro 10.5.2)
-      Real al, ar, el, er;
-      Real cl = eos.SoundSpeed(wli);
-      Real cr = eos.SoundSpeed(wri);
-      el = wli[IPR] * igm1 +
-           0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
-      er = wri[IPR] * igm1 +
-           0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
-      Real rhoa = .5 * (wli[IDN] + wri[IDN]); // average density
-      Real ca = .5 * (cl + cr);               // average sound speed
-      Real pmid = .5 * (wli[IPR] + wri[IPR] + (wli[IV1] - wri[IV1]) * rhoa * ca);
-
-      //--- Step 3.  Compute sound speed in L,R
-      Real ql, qr;
-      ql = (pmid <= wli[IPR])
-               ? 1.0
-               : (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wli[IPR] - 1.0));
-      qr = (pmid <= wri[IPR])
-               ? 1.0
-               : (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wri[IPR] - 1.0));
-
-      //--- Step 4. Compute the max/min wave speeds based on L/R states
-
-      al = wli[IV1] - cl * ql;
-      ar = wri[IV1] + cr * qr;
-
-      Real bp = ar > 0.0 ? ar : 0.0;
-      Real bm = al < 0.0 ? al : 0.0;
+      //--- Step 2.  Compute Roe-averaged state
+      Real sqrtdl = std::sqrt(wli[IDN]);
+      Real sqrtdr = std::sqrt(wri[IDN]);
+      Real isdlpdr = 1.0 / (sqrtdl + sqrtdr);
+
+      wroe[IDN] = sqrtdl * sqrtdr;
+      wroe[IV1] = (sqrtdl * wli[IV1] + sqrtdr * wri[IV1]) * isdlpdr;
+      wroe[IV2] = (sqrtdl * wli[IV2] + sqrtdr * wri[IV2]) * isdlpdr;
+      wroe[IV3] = (sqrtdl * wli[IV3] + sqrtdr * wri[IV3]) * isdlpdr;
+
+      // Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
+      // rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
+      const Real el = wli[IPR] * igm1 +
+                      0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
+      const Real er = wri[IPR] * igm1 +
+                      0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
+      const Real hroe = ((el + wli[IPR]) / sqrtdl + (er + wri[IPR]) / sqrtdr) * isdlpdr;
+
+      //--- Step 3.  Compute sound speed in L,R, and Roe-averaged states
+
+      const Real cl = eos.SoundSpeed(wli);
+      const Real cr = eos.SoundSpeed(wri);
+      Real q = hroe - 0.5 * (SQR(wroe[IV1]) + SQR(wroe[IV2]) + SQR(wroe[IV3]));
+      const Real a = (q < 0.0) ? 0.0 : std::sqrt(gm1 * q);
+
+      //--- Step 4. Compute the max/min wave speeds based on L/R and Roe-averaged values
+      const Real al = std::min((wroe[IV1] - a), (wli[IV1] - cl));
+      const Real ar = std::max((wroe[IV1] + a), (wri[IV1] + cr));
+
+      Real bp = ar > 0.0 ? ar : TINY_NUMBER;
+      Real bm = al < 0.0 ? al : TINY_NUMBER;
 
       //-- Step 5. Compute L/R fluxes along lines bm/bp: F_L - (S_L)U_L; F_R - (S_R)U_R
       Real vxl = wli[IV1] - bm;

src/pgen/cluster.cpp

@@ -348,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
@@ -808,7 +813,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
 }
 
 void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin,
-                          const parthenon::SimTime &) {
+                          const parthenon::SimTime & /*tm*/) {
   // get hydro
   auto pkg = pmb->packages.Get("Hydro");
   const Real gam = pin->GetReal("hydro", "gamma");
@@ -823,6 +828,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");
@@ -849,8 +856,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);

src/pgen/cluster/agn_feedback.cpp

@@ -361,8 +361,6 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
             ///////////////////////////////////////////////////////////////////
 
             eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
-            const Real old_specific_internal_e =
-                prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.));
 
             cons(IDN, k, j, i) += jet_density;
             cons(IM1, k, j, i) += jet_momentum * sign_jet * jet_axis_x;
@@ -379,12 +377,6 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
             }
 
             eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
-            const Real new_specific_internal_e =
-                prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.));
-            PARTHENON_REQUIRE(
-                new_specific_internal_e > jet_specific_internal_e ||
-                    new_specific_internal_e > old_specific_internal_e,
-                "Kinetic injection leads to temperature below jet and existing gas");
           }
 
           // Apply velocity ceiling

src/pgen/cluster/cluster_clips.cpp

@@ -70,6 +70,7 @@ void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
     const Real vAceil2 = SQR(vAceil);
     const Real gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
 
+    const bool magnetic_fields = (hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd);
     Real added_dfloor_mass = 0.0, removed_vceil_energy = 0.0, added_vAceil_mass = 0.0,
          removed_eceil_energy = 0.0;
 
@@ -123,7 +124,7 @@ void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
               }
             }
 
-            if (vAceil2 < std::numeric_limits<Real>::infinity()) {
+            if (magnetic_fields && vAceil2 < std::numeric_limits<Real>::infinity()) {
               // Apply Alfven velocity ceiling by raising density
               const Real rho = prim(IDN, k, j, i);
               const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) +

src/pgen/cluster/magnetic_tower.hpp

@@ -67,7 +67,7 @@ class MagneticTowerObj {
       // Compute the potential in jet cylindrical coordinates
       a_r = 0.0;
       a_theta = 0.0;
-      if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) {
+      if (fabs(h) >= offset_ && fabs(h) <= offset_ + thickness_) {
         a_h = field_ * l_scale_ * exp_r2_h2;
       } else {
         a_h = 0.0;
@@ -110,7 +110,7 @@ class MagneticTowerObj {
       const parthenon::Real exp_r2_h2 = exp(-pow(r / l_scale_, 2));
       // Compute the field in jet cylindrical coordinates
       b_r = 0.0;
-      if (fabs(h) >= 0.001 && fabs(h) <= offset_ + thickness_) {
+      if (fabs(h) >= offset_ && fabs(h) <= offset_ + thickness_) {
         b_theta = 2.0 * field_ * r / l_scale_ * exp_r2_h2;
       } else {
         b_theta = 0.0;