From 0d141db13288d3d5bc2c59d51b3baa757b45f040 Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Thu, 19 Oct 2023 12:20:34 -0600 Subject: [PATCH] Pull new Parthenon (solvers, fixes, outputs). Fork old BiCGStab --- external/parthenon | 2 +- kharma/b_cleanup/b_cleanup.cpp | 3 - kharma/b_cleanup/bicgstab_solver.hpp | 3 +- kharma/b_cleanup/solver_utils.hpp | 176 +++++++++++++++++++++++++++ kharma/b_cleanup/solvers.cpp | 22 ++++ kharma/decs.hpp | 14 ++- kharma/main.cpp | 7 +- 7 files changed, 216 insertions(+), 11 deletions(-) create mode 100644 kharma/b_cleanup/solver_utils.hpp create mode 100644 kharma/b_cleanup/solvers.cpp diff --git a/external/parthenon b/external/parthenon index 72a97564..67cbb148 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 72a975647e5548fee643952a52f12a249fc2b325 +Subproject commit 67cbb1485400051ad94d2f96735e03de76308f07 diff --git a/kharma/b_cleanup/b_cleanup.cpp b/kharma/b_cleanup/b_cleanup.cpp index 9cca1a7b..51276f26 100644 --- a/kharma/b_cleanup/b_cleanup.cpp +++ b/kharma/b_cleanup/b_cleanup.cpp @@ -62,9 +62,6 @@ void B_Cleanup::CleanupDivergence(std::shared_ptr>& 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 B_Cleanup::Initialize(ParameterInput *pin, std::shared_ptr& packages) { auto pkg = std::make_shared("B_Cleanup"); diff --git a/kharma/b_cleanup/bicgstab_solver.hpp b/kharma/b_cleanup/bicgstab_solver.hpp index dc4fe559..8721c8ff 100644 --- a/kharma/b_cleanup/bicgstab_solver.hpp +++ b/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 { diff --git a/kharma/b_cleanup/solver_utils.hpp b/kharma/b_cleanup/solver_utils.hpp new file mode 100644 index 00000000..b7185e60 --- /dev/null +++ b/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 +#include + +#include "basic_types.hpp" +#include "kokkos_abstraction.hpp" + +namespace parthenon { + +namespace solvers { + +struct SparseMatrixAccessor { + ParArray1D ioff, joff, koff; + ParArray1D ioff_inv, joff_inv, koff_inv; + ParArray1D 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> 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 + 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 + 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 +struct Stencil { + ParArray1D w; + ParArray1D ioff, joff, koff; + const int nstencil; + int ndiag; + Stencil() : nstencil(0), ndiag(0) {} + Stencil(const Stencil &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 wgt, + std::vector> 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 + 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 + 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_ diff --git a/kharma/b_cleanup/solvers.cpp b/kharma/b_cleanup/solvers.cpp new file mode 100644 index 00000000..f61284c5 --- /dev/null +++ b/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 diff --git a/kharma/decs.hpp b/kharma/decs.hpp index 8048b2f2..005d4b4e 100644 --- a/kharma/decs.hpp +++ b/kharma/decs.hpp @@ -135,7 +135,6 @@ KOKKOS_INLINE_FUNCTION int dir_of(const Loci loc) } } -#ifdef MPI_PARALLEL /** * Am I rank 0? Saves typing vs comparing the global every time */ @@ -143,12 +142,17 @@ 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 diff --git a/kharma/main.cpp b/kharma/main.cpp index 4d319fdc..1c862262 100644 --- a/kharma/main.cpp +++ b/kharma/main.cpp @@ -174,7 +174,7 @@ int main(int argc, char *argv[]) const int &verbose = pmesh->packages.Get("Globals")->Param("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