Merge pull request #124 from lanl/jmm/update-submodules

update submodules

CMakeLists.txt

@@ -99,6 +99,9 @@ if(PHOEBUS_ENABLE_CUDA)
   set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "" FORCE)
   set(Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE ON CACHE BOOL "" FORCE)
   set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "" FORCE)
+  # Get around a bug with atomics on Cray
+  # This option will eventually be deprecated, so care required.
+  set(Kokkos_ENABLE_IMPL_DESUL_ATOMICS OFF CACHE BOOL "" FORCE)
 endif()
 
 # HDF5
@@ -159,6 +162,10 @@ set(SINGULARITY_SUBMODULE_MODE ON CACHE BOOL "" FORCE)
 set(SINGULARITY_USE_KOKKOSKERNELS OFF CACHE BOOL "" FORCE)
 set(SINGULARITY_BUILD_CLOSURE OFF CACHE BOOL "" FORCE)
 
+# Kill cmake's package registry because it can interfere
+set(CMAKE_FIND_PACKAGE_NO_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+set(CMAKE_FIND_PACKAGE_NO_SYSTEM_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+
 # singularity eos
 message("Configuring singularity-eos")
 add_subdirectory(external/singularity-eos singularity-eos)
@@ -173,14 +180,6 @@ include(cmake/Geometry.cmake)
 # Fluid
 include(cmake/Fluid.cmake)
 
-# Print RHS
-option(PHOEBUS_PRINT_RHS "Print right hand sides" OFF)
-if (PHOEBUS_PRINT_RHS)
-   set(PRINT_RHS 1)
-else()
-   set(PRINT_RHS 0)
-endif()
-
 # Phoebus src
 message("\nConfiguring src")
 add_subdirectory(src)

src/CMakeLists.txt

@@ -156,8 +156,6 @@ endif()
 # target_compile_options(phoebus PUBLIC -G -lineinfo)
 
 target_link_libraries(phoebus PRIVATE parthenon
-  singularity-eos::flags
-  singularity-eos
+  singularity-eos::singularity-eos
   singularity-opac::singularity-opac
-  $<$<BOOL:${PHOEBUS_ENABLE_HDF5}>:singularity-eos::hdf5>
 )

src/compile_constants.hpp.in

@@ -28,6 +28,4 @@
 
 #define CON2PRIM_STATISTICS @CON2PRIM_STATISTICS@
 
-#define PRINT_RHS @PRINT_RHS@
-
 #endif //COMPILE_CONSTANTS_HPP_

src/fluid/fluid.cpp

@@ -493,40 +493,6 @@ TaskStatus ConservedToPrimitiveClassic(T *rc, const IndexRange &ib, const IndexR
   return TaskStatus::complete;
 }
 
-#if SET_FLUX_SRC_DIAGS
-TaskStatus CopyFluxDivergence(MeshBlockData<Real> *rc) {
-  auto *pmb = rc->GetParentPointer().get();
-  auto &fluid = pmb->packages.Get("fluid");
-  if (!fluid->Param<bool>("active")) return TaskStatus::complete;
-
-  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
-  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
-  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
-
-  std::vector<std::string> vars(
-      {fluid_cons::density, fluid_cons::momentum, fluid_cons::energy});
-  PackIndexMap imap;
-  auto divf = rc->PackVariables(vars, imap);
-  const int crho = imap[fluid_cons::density].first;
-  const int cmom_lo = imap[fluid_cons::momentum].first;
-  const int cmom_hi = imap[fluid_cons::momentum].second;
-  const int ceng = imap[fluid_cons::energy].first;
-  std::vector<std::string> diag_vars({diagnostic_variables::divf});
-  auto diag = rc->PackVariables(diag_vars);
-
-  parthenon::par_for(
-      DEFAULT_LOOP_PATTERN, "CopyDivF", DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s,
-      ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) {
-        diag(0, k, j, i) = divf(crho, k, j, i);
-        diag(1, k, j, i) = divf(cmom_lo, k, j, i);
-        diag(2, k, j, i) = divf(cmom_lo + 1, k, j, i);
-        diag(3, k, j, i) = divf(cmom_lo + 2, k, j, i);
-        diag(4, k, j, i) = divf(ceng, k, j, i);
-      });
-  return TaskStatus::complete;
-}
-#endif
-
 TaskStatus CalculateFluidSourceTerms(MeshBlockData<Real> *rc,
                                      MeshBlockData<Real> *rc_src) {
   constexpr int ND = Geometry::NDFULL;

src/fluid/fluid.hpp

@@ -20,6 +20,8 @@
 #include <utils/error_checking.hpp>
 using namespace parthenon::package::prelude;
 
+#include "phoebus_utils/variables.hpp"
+
 namespace fluid {
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
@@ -40,18 +42,56 @@ TaskStatus ConservedToPrimitiveVanDerHolst(T *rc, const IndexRange &ib,
                                            const IndexRange &jb, const IndexRange &kb);
 TaskStatus CalculateFluidSourceTerms(MeshBlockData<Real> *rc,
                                      MeshBlockData<Real> *rc_src);
-TaskStatus CopyFluxDivergence(MeshBlockData<Real> *rc);
 TaskStatus CalculateFluxes(MeshBlockData<Real> *rc);
 TaskStatus FluxCT(MeshBlockData<Real> *rc);
 TaskStatus CalculateDivB(MeshBlockData<Real> *rc);
 Real EstimateTimestepBlock(MeshBlockData<Real> *rc);
 
+template <class T>
+TaskStatus CopyFluxDivergence(T *rc);
+
 template <typename T>
 using c2p_type = TaskStatus (*)(T *, const IndexRange &, const IndexRange &,
                                 const IndexRange &);
 using c2p_meshblock_type = c2p_type<MeshBlockData<Real>>;
 using c2p_mesh_type = c2p_type<MeshData<Real>>;
 
+#if SET_FLUX_SRC_DIAGS
+template <class T>
+TaskStatus CopyFluxDivergence(T *rc) {
+  auto pmb = rc->GetParentPointer();
+  auto &fluid = pmb->packages.Get("fluid");
+  const Params &params = fluid->AllParams();
+  if (!params.Get<bool>("active")) return TaskStatus::complete;
+
+  std::vector<std::string> vars(
+      {fluid_cons::density, fluid_cons::momentum, fluid_cons::energy});
+  PackIndexMap imap;
+  auto divf = rc->PackVariables(vars, imap);
+  const int crho = imap[fluid_cons::density].first;
+  const int cmom_lo = imap[fluid_cons::momentum].first;
+  const int cmom_hi = imap[fluid_cons::momentum].second;
+  const int ceng = imap[fluid_cons::energy].first;
+  std::vector<std::string> diag_vars({diagnostic_variables::divf});
+  auto diag = rc->PackVariables(diag_vars);
+
+  // TODO(JMM): If we expose a way to get cellbounds from the mesh or
+  // meshdata object, that would be better.
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "CopyDivF", DevExecSpace(), 0, divf.GetDim(5) - 1, 0,
+      divf.GetDim(3) - 1, 0, divf.GetDim(2) - 1, 0, divf.GetDim(1) - 1,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        diag(b, 0, k, j, i) = divf(b, crho, k, j, i);
+        diag(b, 1, k, j, i) = divf(b, cmom_lo, k, j, i);
+        diag(b, 2, k, j, i) = divf(b, cmom_lo + 1, k, j, i);
+        diag(b, 3, k, j, i) = divf(b, cmom_lo + 2, k, j, i);
+        diag(b, 4, k, j, i) = divf(b, ceng, k, j, i);
+      });
+  return TaskStatus::complete;
+}
+
+#endif
+
 } // namespace fluid
 
 #endif // FLUID_HPP_

src/geometry/boyer_lindquist.hpp

@@ -77,7 +77,6 @@ class BoyerLindquist {
     const Real th = X2;
     Real r2, r3, sth, cth, DD, mu;
     ComputeDeltaMu_(r, th, r2, r3, sth, cth, DD, mu);
-    const Real r3dm = r3 * DD * mu;
     const Real sth2 = sth * sth;
     LinearAlgebra::SetZero(g, NDFULL, NDFULL);
     g[0][0] = -(1.0 - ratio(2.0, r * mu));

src/monopole_gr/interp_3d_to_1d.hpp

@@ -225,9 +225,9 @@ GetMonopoleVarsHelper(const EnergyMomentum &tmunu, const Geometry_t &geom,
     const Real X2 = coords.x2v(j);
     const Real X3 = coords.x3v(k);
     transform(X1, X2, X3, C, c2s, s2c);
-    const Real r = C[0];
-    const Real th = C[1];
-    const Real ph = C[2];
+    // const Real r = C[0];
+    // const Real th = C[1];
+    // const Real ph = C[2];
 
     // J
     jr = 0;

src/pgen/homogeneous_sphere.cpp

@@ -52,7 +52,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const Real vx = pin->GetOrAddReal("homogeneous_sphere", "vx", 0.0);
 
   auto &coords = pmb->coords;
-  auto pmesh = pmb->pmy_mesh;
 
   IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire);
   IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire);

src/pgen/radiation_advection.cpp

@@ -91,10 +91,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto rad = pmb->packages.Get("radiation").get();
   auto species = rad->Param<std::vector<singularity::RadiationType>>("species");
   auto num_species = rad->Param<int>("num_species");
-  singularity::RadiationType species_d[3] = {};
-  for (int s = 0; s < num_species; s++) {
-    species_d[s] = species[s];
-  }
 
   const Real W = 1 / sqrt(1 - vx * vx);
   const Real t0p = 1.5 * kappa * width * width;

src/pgen/radiation_equilibration.cpp

@@ -52,10 +52,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const auto specB = idJ.GetBounds(1);
   const Real J = pin->GetOrAddReal("radiation_equilibration", "J", 0.0);
 
-  auto &coords = pmb->coords;
-  auto pmesh = pmb->pmy_mesh;
-  int ndim = pmesh->ndim;
-
   IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire);
   IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire);
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire);

src/pgen/rotor.cpp

@@ -64,9 +64,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto emax = pmb->packages.Get("eos")->Param<Real>("sie_max");
   auto gpkg = pmb->packages.Get("geometry");
   auto geom = Geometry::GetCoordinateSystem(rc.get());
-  Real a_snake, k_snake, alpha, betax, betay, betaz;
+  Real a_snake, k_snake, alpha, betay;
   alpha = 1;
-  a_snake = k_snake = betax = betay = betaz = 0;
+  a_snake = k_snake = betay = 0;
   Real tf = pin->GetReal("parthenon/time", "tlim");
   if (is_snake) {
     a_snake = gpkg->Param<Real>("a");