Pull new Parthenon (solvers, fixes, outputs). Fork old BiCGStab

kharma/b_cleanup/b_cleanup.cpp

@@ -62,9 +62,6 @@ void B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md) {}
 using namespace parthenon;
 using namespace parthenon::solvers;
 
-// TODO get the transport manager working later
-// Needs a call every X steps option, probably return a TaskList or TaskRegion
-
 std::shared_ptr<KHARMAPackage> B_Cleanup::Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages)
 {
     auto pkg = std::make_shared<KHARMAPackage>("B_Cleanup");

kharma/b_cleanup/bicgstab_solver.hpp

@@ -22,10 +22,11 @@
 #include "interface/meshblock_data.hpp"
 #include "interface/state_descriptor.hpp"
 #include "kokkos_abstraction.hpp"
-#include "solvers/solver_utils.hpp"
 #include "tasks/task_id.hpp"
 #include "tasks/task_list.hpp"
 
+#include "solver_utils.hpp"
+
 namespace parthenon {
 
 namespace solvers {

kharma/b_cleanup/solver_utils.hpp

@@ -0,0 +1,176 @@
+//========================================================================================
+// (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef SOLVERS_SOLVER_UTILS_HPP_
+#define SOLVERS_SOLVER_UTILS_HPP_
+
+#include <string>
+#include <vector>
+
+#include "basic_types.hpp"
+#include "kokkos_abstraction.hpp"
+
+namespace parthenon {
+
+namespace solvers {
+
+struct SparseMatrixAccessor {
+  ParArray1D<int> ioff, joff, koff;
+  ParArray1D<int> ioff_inv, joff_inv, koff_inv;
+  ParArray1D<int> inv_entries;
+
+  const int nstencil;
+  int ndiag;
+  SparseMatrixAccessor() : nstencil(0), ndiag(0) {}
+  SparseMatrixAccessor(const SparseMatrixAccessor &sp)
+      : ioff(sp.ioff), joff(sp.joff), koff(sp.koff), nstencil(sp.nstencil),
+        ndiag(sp.ndiag), inv_entries(sp.inv_entries) {}
+  SparseMatrixAccessor(const std::string &label, const int n,
+                       std::vector<std::vector<int>> off)
+      : ioff(label + "_ioff", n), joff(label + "_joff", n), koff(label + "_koff", n),
+        ioff_inv(label + "_ioff_inv", n), joff_inv(label + "_joff_inv", n),
+        koff_inv(label + "_koff_inv", n), inv_entries(label + "_inv_ent", n),
+        nstencil(n) {
+    PARTHENON_REQUIRE_THROWS(off.size() == 3,
+                             "Offset array must have dimensions off[3][*]");
+    PARTHENON_REQUIRE_THROWS(off[0].size() >= n, "Offset array off[0][*] too small");
+    PARTHENON_REQUIRE_THROWS(off[1].size() >= n, "Offset array off[1][*] too small");
+    PARTHENON_REQUIRE_THROWS(off[2].size() >= n, "Offset array off[2][*] too small");
+    auto ioff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), ioff);
+    auto joff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), joff);
+    auto koff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), koff);
+
+    auto ioff_inv_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), ioff_inv);
+    auto joff_inv_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), joff_inv);
+    auto koff_inv_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), koff_inv);
+
+    auto inv_ent_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), inv_entries);
+
+    for (int i = 0; i < n; i++) {
+      ioff_h(i) = off[0][i];
+      joff_h(i) = off[1][i];
+      koff_h(i) = off[2][i];
+      // this is inverse.
+      ioff_inv_h(i) = -off[0][i];
+      joff_inv_h(i) = -off[1][i];
+      koff_inv_h(i) = -off[2][i];
+
+      if (off[0][i] == 0 && off[1][i] == 0 && off[2][i] == 0) {
+        ndiag = i;
+      }
+    }
+    for (int i = 0; i < n; i++) {
+      for (int j = 0; j < n; j++) {
+        if (ioff_h(i) == ioff_inv_h(j) && joff_h(i) == joff_inv_h(j) &&
+            koff_h(i) == koff_inv_h(j)) {
+          inv_entries(i) = j;
+          std::cout << "inv_entries:" << i << " " << j << std::endl;
+        }
+      } // j
+    }   // i
+
+    Kokkos::deep_copy(ioff, ioff_h);
+    Kokkos::deep_copy(joff, joff_h);
+    Kokkos::deep_copy(koff, koff_h);
+
+    Kokkos::deep_copy(inv_entries, inv_ent_h);
+  }
+
+  template <typename PackType>
+  KOKKOS_INLINE_FUNCTION Real MatVec(const PackType &spmat, const int imat_lo,
+                                     const int imat_hi, const PackType &v, const int iv,
+                                     const int b, const int k, const int j,
+                                     const int i) const {
+    Real matvec = 0.0;
+    for (int n = imat_lo; n <= imat_hi; n++) {
+      const int m = n - imat_lo;
+      matvec += spmat(b, n, k, j, i) * v(b, iv, k + koff(m), j + joff(m), i + ioff(m));
+    }
+    return matvec;
+  }
+
+  template <typename PackType>
+  KOKKOS_INLINE_FUNCTION Real Jacobi(const PackType &spmat, const int imat_lo,
+                                     const int imat_hi, const PackType &v, const int iv,
+                                     const int b, const int k, const int j, const int i,
+                                     const Real rhs) const {
+    const Real matvec = MatVec(spmat, imat_lo, imat_hi, v, iv, b, k, j, i);
+    return (rhs - matvec + spmat(b, imat_lo + ndiag, k, j, i) * v(b, iv, k, j, i)) /
+           spmat(b, imat_lo + ndiag, k, j, i);
+  }
+};
+
+template <typename T>
+struct Stencil {
+  ParArray1D<T> w;
+  ParArray1D<int> ioff, joff, koff;
+  const int nstencil;
+  int ndiag;
+  Stencil() : nstencil(0), ndiag(0) {}
+  Stencil(const Stencil<T> &st)
+      : w(st.w), ioff(st.ioff), joff(st.joff), koff(st.koff), nstencil(st.nstencil),
+        ndiag(st.ndiag) {}
+  Stencil(const std::string &label, const int n, std::vector<T> wgt,
+          std::vector<std::vector<int>> off)
+      : w(label + "_w", n), ioff(label + "_ioff", n), joff(label + "_joff", n),
+        koff(label + "_koff", n), nstencil(n) {
+    PARTHENON_REQUIRE_THROWS(off.size() == 3,
+                             "Offset array must have dimensions off[3][*]");
+    PARTHENON_REQUIRE_THROWS(wgt.size() >= n, "Weight array too small");
+    PARTHENON_REQUIRE_THROWS(off[0].size() >= n, "Offset array off[0][*] too small");
+    PARTHENON_REQUIRE_THROWS(off[1].size() >= n, "Offset array off[1][*] too small");
+    PARTHENON_REQUIRE_THROWS(off[2].size() >= n, "Offset array off[2][*] too small");
+    auto w_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), w);
+    auto ioff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), ioff);
+    auto joff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), joff);
+    auto koff_h = Kokkos::create_mirror_view(Kokkos::HostSpace(), koff);
+
+    for (int i = 0; i < n; i++) {
+      w_h(i) = wgt[i];
+      ioff_h(i) = off[0][i];
+      joff_h(i) = off[1][i];
+      koff_h(i) = off[2][i];
+      if (off[0][i] == 0 && off[1][i] == 0 && off[2][i] == 0) {
+        ndiag = i;
+      }
+    }
+
+    Kokkos::deep_copy(w, w_h);
+    Kokkos::deep_copy(ioff, ioff_h);
+    Kokkos::deep_copy(joff, joff_h);
+    Kokkos::deep_copy(koff, koff_h);
+  }
+
+  template <typename PackType>
+  KOKKOS_INLINE_FUNCTION Real MatVec(const PackType &v, const int iv, const int b,
+                                     const int k, const int j, const int i) const {
+    Real matvec = 0.0;
+    for (int n = 0; n < nstencil; n++) {
+      matvec += w(n) * v(b, iv, k + koff(n), j + joff(n), i + ioff(n));
+    }
+    return matvec;
+  }
+
+  template <typename PackType>
+  KOKKOS_INLINE_FUNCTION Real Jacobi(const PackType &v, const int iv, const int b,
+                                     const int k, const int j, const int i,
+                                     const Real rhs) const {
+    const Real matvec = MatVec(v, iv, b, k, j, i);
+    return (rhs - matvec + w(ndiag) * v(b, iv, k, j, i)) / w(ndiag);
+  }
+};
+
+} // namespace solvers
+
+} // namespace parthenon
+
+#endif // SOLVERS_SOLVER_UTILS_HPP_

kharma/b_cleanup/solvers.cpp

@@ -0,0 +1,22 @@
+//========================================================================================
+// (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include "bicgstab_solver.hpp"
+
+namespace parthenon {
+namespace solvers {
+
+int BiCGStabCounter::global_num_bicgstab_solvers = 0;
+
+} // namespace solvers
+} // namespace parthenon

kharma/decs.hpp

@@ -135,20 +135,24 @@ KOKKOS_INLINE_FUNCTION int dir_of(const Loci loc)
     }
 }
 
-#ifdef MPI_PARALLEL
 /**
  * Am I rank 0?  Saves typing vs comparing the global every time
  */
 inline bool MPIRank0()
 {
     return (parthenon::Globals::my_rank == 0 ? true : false);
 }
-#else
 /**
- * DUMMY version for no-MPI case: constexpr return for slight optimizations.
+ * Numbers I could just get as globals, but renamed for consistency
  */
-inline bool MPIRank0() { return true; }
-#endif // MPI_PARALLEL
+inline int MPINumRanks()
+{
+    return parthenon::Globals::nranks;
+}
+inline int MPIMyRank()
+{
+    return parthenon::Globals::my_rank;
+}
 
 // A few generic "NDArray" overloads for readability.
 // TODO torn on futures of these: they're explicitly per-block

kharma/main.cpp

@@ -174,7 +174,7 @@ int main(int argc, char *argv[])
     const int &verbose = pmesh->packages.Get("Globals")->Param<int>("verbose");
     if(MPIRank0() && verbose > 0) {
         // Print a list of variables as Parthenon used to (still does by default)
-        std::cout << "#Variables in use:\n" << *(pmesh->resolved_packages) << std::endl;
+        std::cout << "Variables in use:\n" << *(pmesh->resolved_packages) << std::endl;
 
         // Print a list of all loaded packages.  Surprisingly useful for debugging init logic
         std::cout << "Packages in use: " << std::endl;
@@ -183,6 +183,11 @@ int main(int argc, char *argv[])
         }
         std::cout << std::endl;
 
+        // Print the number of meshblocks and ranks in use
+        std::cout << "Running with " << pmesh->block_list.size() << " total meshblocks, " << MPINumRanks() << " MPI ranks." << std::endl;
+        // TODO could print entire distribution if it gets interesting
+        std::cout << "Blocks on rank " << MPIMyRank() << ": " << pmesh->GetNumMeshBlocksThisRank() << "\n" << std::endl;
+
         // Write all parameters etc. to console if we should be especially wordy
         if ((verbose > 1) && MPIRank0()) {
             // This dumps the full Kokkos config, useful for double-checking