Merge pull request #21 from AFD-Illinois/dev

KHARMA 2022.3

.gitignore

@@ -12,8 +12,9 @@ frames_*/
 *.rhdf
 *.xdmf
 *.hst
-# Archival parsed parameters file
+# Archival files
 kharma_parsed_*.par
+log_*.txt
 
 # Editor documents
 .project

.gitlab-ci.yml

@@ -34,6 +34,7 @@ bondi:
     paths:
       - tests/bondi/*.png
       - tests/bondi/*.hst
+      - tests/bondi/*.txt
 
 mhdmodes:
   stage: tests
@@ -47,6 +48,19 @@ mhdmodes:
     paths:
       - tests/mhdmodes/*.png
       - tests/mhdmodes/*.hst
+      - tests/mhdmodes/*.txt
+
+emhdmodes:
+  stage: tests
+  before_script:
+    - cd tests/emhdmodes/
+  script:
+    - bash run.sh
+    - bash check.sh
+  artifacts:
+    when: always
+    paths:
+      - tests/emhdmodes/*.png
 
 noh:
   stage: tests

.gitmodules

@@ -5,3 +5,6 @@
 [submodule "external/variant"]
 	path = external/variant
 	url = https://github.com/mpark/variant.git
+[submodule "external/kokkos-kernels"]
+	path = external/kokkos-kernels
+	url = https://github.com/kokkos/kokkos-kernels

CMakeLists.txt

@@ -13,40 +13,56 @@ set(CMAKE_CXX_STANDARD 14)
 #set(CMAKE_CXX17_STANDARD_COMPILE_OPTION "-std=c++17")
 #set(PARTHENON_ENABLE_CPP17 ON CACHE BOOL "KHARMA Override")
 
-#foreach(path ${CMAKE_PREFIX_PATH})
-#  include_directories(${path})
-#endforeach()
-
 # Parthenon options
 set(PARTHENON_DISABLE_EXAMPLES ON CACHE BOOL "KHARMA Override")
 set(PARTHENON_LINT_DEFAULT OFF CACHE BOOL "KHARMA Override")
-# To use old Summit built-in HDF5
-set(PARTHENON_DISABLE_HDF5_COMPRESSION ON CACHE BOOL "KHARMA Override")
+# Attempt HDF5 compression, requires recent/standard HDF5. YMMV
+set(PARTHENON_DISABLE_HDF5_COMPRESSION OFF CACHE BOOL "KHARMA Override")
 
 # Parthenon internal build options
 set(BUILD_TESTING OFF CACHE BOOL "KHARMA Override")
 set(ENABLE_COMPILER_WARNINGS OFF CACHE BOOL "KHARMA Override")
+# TODO set this here when I upstream or otherwise unfork
+#set(COORDINATE_TYPE GRCoordinates)
 
 # Kokkos options
-set(Kokkos_ENABLE_OPENMP ON)
-set(Kokkos_ENABLE_CUDA_LAMBDA ON)
-set(Kokkos_ENABLE_CUDA_CONSTEXPR ON)
-set(Kokkos_ENABLE_HWLOC OFF) # Possible speed improvement?
-set(Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION ON)
+set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "KHARMA Override")
+set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "KHARMA Override")
+set(Kokkos_ENABLE_CUDA_CONSTEXPR ON CACHE BOOL "KHARMA Override")
+set(Kokkos_ENABLE_HWLOC OFF CACHE BOOL "KHARMA Override") # Possible speed improvement?
+set(Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION ON CACHE BOOL "KHARMA Override")
+
+# Build only what we need of KokkosKernels
+set(KokkosKernels_ENABLE_TPL_CUSPARSE OFF CACHE BOOL "KHARMA Override")
+set(KokkosKernels_ENABLE_TPL_CUBLAS OFF CACHE BOOL "KHARMA Override")
 
-# TODO set this here when I upstream or otherwise unfork
-#set(COORDINATE_TYPE GRCoordinates)
 
 # Parthenon says it doesn't need MPI.  It just *strongly prefers* it, and so do we.
+# Builds without MPI have pretty limited support, you can usually find distribution 
+# packages or other ways to install it on personal machines without too much work.
+# Check out oneAPI or NVHPC for software distributions that include easily-usable,
+# fast MPI modules
+# If you really want to disable MPI, set this to ON and comment the next two lines
+set(PARTHENON_DISABLE_MPI OFF CACHE BOOL "KHARMA Override")
 find_package(MPI REQUIRED)
 include_directories(SYSTEM ${MPI_INCLUDE_PATH})
+
+# OpenMP is strictly required
 find_package(OpenMP REQUIRED)
 
-# TODO don't build parthenon unit tests etc just the library
+# Build Parthenon
 add_subdirectory(external/parthenon)
 include_directories(external/parthenon/src)
 # mpark::variant is header only, don't build anything
 include_directories(external/variant/include)
+# Kokkos kernels: don't build them (very slow), just import all headers
+# Requires KokkosKernels_config.h shipped with KHARMA, YMMV
+# In case of issues, uncomment the following line to build them
+#add_subdirectory(external/kokkos-kernels)
+include_directories(external/kokkos-kernels/src)
+include_directories(external/kokkos-kernels/src/batched)
+include_directories(external/kokkos-kernels/src/batched/dense)
+include_directories(external/kokkos-kernels/src/batched/dense/impl)
 
-# KHARMA folder
+# Finally, build KHARMA
 add_subdirectory(kharma)

kharma/CMakeLists.txt

@@ -13,37 +13,47 @@ AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR} EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/prob EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/coordinates EXE_NAME_SRC)
 
-AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_cd EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/current EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/electrons EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/emhd EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/floors EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/grmhd EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/implicit EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/reductions EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/emhd EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/wind EXE_NAME_SRC)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/prob)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/coordinates)
 
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_cd)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/current)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/electrons)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/emhd)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/floors)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/grmhd)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/implicit)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/reductions)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/emhd)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/wind)
 
 add_executable(${EXE_NAME} ${EXE_NAME_SRC})
 
+target_link_libraries(${EXE_NAME} PUBLIC kokkos)
+# We actually only need the header
+#target_link_libraries(${EXE_NAME} PUBLIC kokkoskernels)
 target_link_libraries(${EXE_NAME} PUBLIC parthenon)
 
 # OPTIONS
 # These are almost universally performance trade-offs
 # TODO is there any way to make compile options less painful in CMake?
 option(FUSE_FLUX_KERNELS "Bundle the usual four flux calculation kernels (floors,R,L,apply) into one" ON)
-option(FUSE_EMF_KERNELS "Bundle the three emf direction kernels into one. Likely won't affect much" ON)
 option(FUSE_FLOOR_KERNELS "Bundle applying the floors and ceilings into one kernel" ON)
 option(FAST_CARTESIAN "Break operation in curved spacetimes to make Cartesian Minkowski space computations faster" OFF)
 if(FUSE_FLUX_KERNELS)

kharma/b_cd/b_cd.cpp

@@ -59,7 +59,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Packages_t pack
     Real damping = pin->GetOrAddReal("b_field", "damping", 0.1);
     params.Add("damping", damping);
 
-    std::vector<int> s_vector({3});
+    std::vector<int> s_vector({NVEC});
 
     MetadataFlag isPrimitive = packages.Get("GRMHD")->Param<MetadataFlag>("PrimitiveFlag");
 

kharma/b_cd/b_cd.hpp

@@ -91,37 +91,4 @@ TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData<Real> *rc);
  */
 void FillOutput(MeshBlock *pmb, ParameterInput *pin);
 
-/**
- * Turn the primitive B field into the local conserved flux
- */
-KOKKOS_INLINE_FUNCTION void prim_to_flux(const GRCoordinates& G, const ScratchPad2D<Real> &P, const VarMap& m_p, const FourVectors D,
-                                           const int& k, const int& j, const int& i, const int& dir,
-                                           ScratchPad2D<Real>& flux, const VarMap& m_u, const Loci& loc=Loci::center)
-{
-    Real gdet = G.gdet(loc, j, i);
-    if (dir == 0) { // Parent is templated on dir, so we should get the speed here still
-        VLOOP flux(m_u.B1 + v, i) = P(m_p.B1 + v, i) * gdet;
-        flux(m_u.PSI, i) = P(m_p.PSI, i) * gdet;
-    } else {
-        // Dual of Maxwell tensor
-        // Dedner would have e.g. P(m.psip, i) * gdet,
-        // but for us this is in the source term
-        VLOOP flux(m_u.B1 + v, i) = (D.bcon[v+1] * D.ucon[dir] - D.bcon[dir] * D.ucon[v+1]) * gdet;
-        // Psi field update as in Mosta et al (IllinoisGRMHD), alternate explanation Jesse et al (2020)
-        //Real alpha = 1. / sqrt(-G.gcon(Loci::center, j, i, 0, 0));
-        //Real beta_dir = G.gcon(Loci::center, j, i, 0, dir) * alpha * alpha;
-        flux(m_u.PSI, i) = (D.bcon[dir] - G.gcon(Loci::center, j, i, 0, dir) * P(m_p.PSI, i)) * gdet;
-    }
-}
-
-KOKKOS_INLINE_FUNCTION void p_to_u(const GRCoordinates& G, const VariablePack<Real>& P, const VarMap& m_p,
-                                           const int& k, const int& j, const int& i,
-                                           const VariablePack<Real>& U, const VarMap& m_u, const Loci& loc=Loci::center)
-{
-    Real gdet = G.gdet(loc, j, i);
-    VLOOP U(m_u.B1 + v, k, j, i) = P(m_p.B1 + v, k, j, i) * gdet;
-    U(m_u.PSI, k, j, i) = P(m_p.PSI, k, j, i) * gdet;
-
-}
-
 }

kharma/b_cd/seed_B_cd.cpp

@@ -39,7 +39,7 @@
 #include "b_field_tools.hpp"
 
 #include "b_flux_ct.hpp"
-#include "mhd_functions.hpp"
+#include "grmhd_functions.hpp"
 
 TaskStatus B_CD::SeedBField(MeshBlockData<Real> *rc, ParameterInput *pin)
 {
@@ -165,12 +165,7 @@ TaskStatus B_CD::SeedBField(MeshBlockData<Real> *rc, ParameterInput *pin)
             B_P(2, k, j, i) = 0.;
         }
     );
-    pmb->par_for("first_U_B", ks, ke, js, je, is, ie,
-        KOKKOS_LAMBDA_3D {
-            // Use the "other" P to U, because we're content that psi = 0 to begin
-            B_FluxCT::p_to_u(G, B_P, k, j, i, B_U);
-        }
-    );
+    B_FluxCT::PtoU(rc);
 
     return TaskStatus::complete;
 }