From 0b60c90a3df82a08eff8b535321094602e80c1e6 Mon Sep 17 00:00:00 2001 From: Marc Henry de Frahan Date: Thu, 11 Jun 2026 21:22:33 +0000 Subject: [PATCH] Update Kokkos policy kernels and timing --- .gitignore | 1 + .../kokkos_fortran_openacc/CMakeLists.txt | 39 ++++ .../various_policy_impls/CMakeLists.txt | 13 ++ .../shallow_mdrangepolicy.cpp | 176 +++++++++-------- .../shallow_rangepolicy.cpp | 180 ++++++++--------- .../shallow_rangepolicy_flat.cpp | 182 ++++++++--------- .../shallow_teampolicy.cpp | 187 ++++++++---------- 7 files changed, 400 insertions(+), 378 deletions(-) create mode 100644 swm_kokkos/kokkos_fortran_openacc/CMakeLists.txt create mode 100644 swm_kokkos/various_policy_impls/CMakeLists.txt diff --git a/.gitignore b/.gitignore index bd53a66..1936b51 100644 --- a/.gitignore +++ b/.gitignore @@ -38,4 +38,5 @@ swm_fortran/swm_fortran_amrex_driver *.out *.bin *.txt +!CMakeLists.txt results_* diff --git a/swm_kokkos/kokkos_fortran_openacc/CMakeLists.txt b/swm_kokkos/kokkos_fortran_openacc/CMakeLists.txt new file mode 100644 index 0000000..2ce5910 --- /dev/null +++ b/swm_kokkos/kokkos_fortran_openacc/CMakeLists.txt @@ -0,0 +1,39 @@ +set(SWM_KOKKOS_FORTRAN_OPENACC_TARGET shallow_mdrangepolicy_fortran_acc) + +if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") + if(SWM_DEVICE STREQUAL "cpu") + set(OPENACC_FORTRAN_COMPILE_FLAGS -fopenacc -foffload=disable) + elseif(SWM_DEVICE STREQUAL "gpu") + set(OPENACC_FORTRAN_COMPILE_FLAGS -fopenacc) + endif() + set(OPENACC_FORTRAN_LINK_FLAGS -fopenacc) +elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC" OR CMAKE_Fortran_COMPILER_ID STREQUAL "NVIDIA") + if(SWM_DEVICE STREQUAL "cpu") + set(OPENACC_FORTRAN_COMPILE_FLAGS -acc=multicore -Minfo=accel -Mnofma) + set(OPENACC_FORTRAN_LINK_FLAGS -acc=multicore) + elseif(SWM_DEVICE STREQUAL "gpu") + set(OPENACC_FORTRAN_COMPILE_FLAGS -acc=gpu -Minfo=accel -Mnofma) + set(OPENACC_FORTRAN_LINK_FLAGS -acc=gpu) + endif() +else() + message(WARNING "OpenACC Fortran flags are not configured for compiler ${CMAKE_Fortran_COMPILER_ID}; directives may be ignored.") + set(OPENACC_FORTRAN_COMPILE_FLAGS "") + set(OPENACC_FORTRAN_LINK_FLAGS "") +endif() + +add_executable(${SWM_KOKKOS_FORTRAN_OPENACC_TARGET} + shallow_mdrangepolicy.cpp + shallow_fortran_acc.F90 +) + +target_link_libraries(${SWM_KOKKOS_FORTRAN_OPENACC_TARGET} Kokkos::kokkos) + +foreach(flag IN LISTS OPENACC_FORTRAN_COMPILE_FLAGS) + target_compile_options(${SWM_KOKKOS_FORTRAN_OPENACC_TARGET} + PRIVATE "$<$:${flag}>" + ) +endforeach() + +target_link_options(${SWM_KOKKOS_FORTRAN_OPENACC_TARGET} + PRIVATE ${OPENACC_FORTRAN_LINK_FLAGS} +) diff --git a/swm_kokkos/various_policy_impls/CMakeLists.txt b/swm_kokkos/various_policy_impls/CMakeLists.txt new file mode 100644 index 0000000..ef1db67 --- /dev/null +++ b/swm_kokkos/various_policy_impls/CMakeLists.txt @@ -0,0 +1,13 @@ +# Kokkos various policy implementations + +set(KOKKOS_TARGETS + shallow_rangepolicy + shallow_rangepolicy_flat + shallow_mdrangepolicy + shallow_teampolicy +) + +foreach(target ${KOKKOS_TARGETS}) + add_executable(${target} ${target}.cpp) + target_link_libraries(${target} Kokkos::kokkos) +endforeach() diff --git a/swm_kokkos/various_policy_impls/shallow_mdrangepolicy.cpp b/swm_kokkos/various_policy_impls/shallow_mdrangepolicy.cpp index 81b2183..a4525da 100644 --- a/swm_kokkos/various_policy_impls/shallow_mdrangepolicy.cpp +++ b/swm_kokkos/various_policy_impls/shallow_mdrangepolicy.cpp @@ -4,6 +4,16 @@ #include #include #include +#if defined(__has_include) +#if __has_include() +#include +#define SWM_HAVE_NVTX 1 +#endif +#endif +#ifndef SWM_HAVE_NVTX +#define nvtxRangePush(name) ((void)0) +#define nvtxRangePop() ((void)0) +#endif #define MIN(x,y) ((x)>(y)?(y):(x)) #define MAX(x,y) ((x)>(y)?(x):(y)) @@ -20,7 +30,7 @@ using Layout = Kokkos::LayoutRight; using ExecSpace = Kokkos::DefaultExecutionSpace; using MemSpace = ExecSpace::memory_space; -using ViewMatrixType = Kokkos::View; +using ViewMatrixType = Kokkos::View>; using HostViewMatrixType = Kokkos::View; void write_to_file(auto array, int tM, int tN, const char *filename); @@ -83,10 +93,11 @@ int main(int argc, char **argv) { double tpi,di,dj,pcf; double tdts8,tdtsdx,tdtsdy,fsdx,fsdy; - int mnmin,ncycle; + int ncycle; // timer variables double ctime,tcyc,time,ptime; + double t100 = 0., t200 = 0., t300 = 0.; // ** Initialisations ** @@ -111,13 +122,13 @@ int main(int argc, char **argv) { pcf = pi * pi * a * a / (el * el); // Initial values of the stream function and p - Kokkos::parallel_for("init_psi_p", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { + Kokkos::parallel_for("init_psi_p", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { psi(i,j) = a * std::sin((i + .5) * di) * std::sin((j + .5) * dj); p(i,j) = pcf * (std::cos(2. * (i) * di) + std::cos(2. * (j) * dj)) + 50000.; }); // Initialize velocities - Kokkos::parallel_for("init_u_v", Kokkos::MDRangePolicy>({0,0}, {M,N}), KOKKOS_LAMBDA(const int i, const int j) { + Kokkos::parallel_for("init_u_v", Kokkos::MDRangePolicy>({0,0}, {M,N}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { u(i+1,j) = -(psi(i+1,j+1) - psi(i+1,j)) / dy; v(i,j+1) = (psi(i+1,j+1) - psi(i,j+1)) / dx; }); @@ -138,7 +149,7 @@ int main(int argc, char **argv) { v(M,0) = v(0,N); }); - Kokkos::parallel_for("init_old_arrays", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { + Kokkos::parallel_for("init_old_arrays", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { uold(i,j) = u(i,j); vold(i,j) = v(i,j); pold(i,j) = p(i,j); @@ -157,26 +168,6 @@ int main(int argc, char **argv) { printf(" grid spacing in the y direction %f\n", dy); printf(" time step %f\n", dt); printf(" time filter parameter %f\n", alpha); - - mnmin = MIN(M,N); - if constexpr (!std::is_same_v) { - Kokkos::deep_copy(u_host, u); - Kokkos::deep_copy(v_host, v); - Kokkos::deep_copy(p_host, p); - } - // printf(" initial diagonal elements of p\n"); - // for (int i=0; i>({1,0}, {M_LEN,N}), KOKKOS_LAMBDA(const int i, const int j) { - cu(i,j) = 0.5 * (p(i,j) + p(i-1,j)) * u(i,j); - }); - - // Compute cv - Kokkos::parallel_for("compute_cv", Kokkos::MDRangePolicy>({0,1}, {M,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - cv(i,j) = 0.5 * (p(i,j) + p(i,j-1)) * v(i,j); - }); - - // Compute z - Kokkos::parallel_for("compute_z", Kokkos::MDRangePolicy>({1,1}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - z(i,j) = (fsdx * (v(i,j) - v(i-1,j)) - fsdy * (u(i,j) - u(i,j-1))) / (p(i-1,j-1) + p(i,j-1) + p(i,j) + p(i-1,j)); - }); - - // Compute h - Kokkos::parallel_for("compute_h", Kokkos::MDRangePolicy>({0,0}, {M,N}), KOKKOS_LAMBDA(const int i, const int j) { - h(i,j) = p(i,j) + 0.25 * (u(i+1,j) * u(i+1,j) + u(i,j) * u(i,j) + v(i,j+1) * v(i,j+1) + v(i,j) * v(i,j)); + nvtxRangePush("UpdateIntermediateVariables"); + Kokkos::Timer timer100; + + // Compute cu, cv, z, h (fused) + Kokkos::parallel_for("compute_cu_cv_z_h", Kokkos::MDRangePolicy>({0,0}, {M,N}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { + double p_ij = p(i,j); + double p_i1j = p(i+1,j); + double p_ij1 = p(i,j+1); + double p_i1j1 = p(i+1,j+1); + double u_ij = u(i,j); + double u_i1j = u(i+1,j); + double u_i1j1 = u(i+1,j+1); + double v_ij = v(i,j); + double v_ij1 = v(i,j+1); + double v_i1j1 = v(i+1,j+1); + cu(i+1,j) = 0.5 * (p_i1j + p_ij) * u_i1j; + cv(i,j+1) = 0.5 * (p_ij1 + p_ij) * v_ij1; + z(i+1,j+1) = (fsdx * (v_i1j1 - v_ij1) - fsdy * (u_i1j1 - u_i1j)) / (p_ij + p_i1j + p_i1j1 + p_ij1); + h(i,j) = p_ij + 0.25 * (u_i1j * u_i1j + u_ij * u_ij + v_ij1 * v_ij1 + v_ij * v_ij); }); // Periodic continuation @@ -231,25 +222,37 @@ int main(int argc, char **argv) { h(M,N) = h(0,0); }); + Kokkos::fence(); + t100 += timer100.seconds(); + nvtxRangePop(); + // Compute new values u,v and p + nvtxRangePush("UpdateNewVariables"); + Kokkos::Timer timer200; tdts8 = tdt / 8.; tdtsdx = tdt / dx; tdtsdy = tdt / dy; - // Compute unew - Kokkos::parallel_for("compute_unew", Kokkos::MDRangePolicy>({1,0}, {M_LEN,N}), KOKKOS_LAMBDA(const int i, const int j) { - unew(i,j) = uold(i,j) + tdts8 * (z(i,j+1) + z(i,j)) * (cv(i,j+1) + cv(i-1,j+1) + cv(i-1,j) + cv(i,j)) - tdtsdx * (h(i,j) - h(i-1,j)); - }); - - // Compute vnew - Kokkos::parallel_for("compute_vnew", Kokkos::MDRangePolicy>({0,1}, {M,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - vnew(i,j) = vold(i,j) - tdts8 * (z(i+1,j) + z(i,j)) * (cu(i+1,j) + cu(i,j) + cu(i,j-1) + cu(i+1,j-1)) - tdtsdy * (h(i,j) - h(i,j-1)); - }); - - // Compute pnew - Kokkos::parallel_for("compute_pnew", Kokkos::MDRangePolicy>({0,0}, {M,N}), KOKKOS_LAMBDA(const int i, const int j) { - pnew(i,j) = pold(i,j) - tdtsdx * (cu(i+1,j) - cu(i,j)) - tdtsdy * (cv(i,j+1) - cv(i,j)); + // Compute unew, vnew, pnew (fused) + Kokkos::parallel_for("compute_unew_vnew_pnew", Kokkos::MDRangePolicy>({0,0}, {M,N}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { + double z_i1j1 = z(i+1,j+1); + double z_i1j = z(i+1,j); + double z_ij1 = z(i,j+1); + double cv_i1j1 = cv(i+1,j+1); + double cv_ij1 = cv(i,j+1); + double cv_ij = cv(i,j); + double cv_i1j = cv(i+1,j); + double cu_i1j1 = cu(i+1,j+1); + double cu_ij1 = cu(i,j+1); + double cu_ij = cu(i,j); + double cu_i1j = cu(i+1,j); + double h_i1j = h(i+1,j); + double h_ij = h(i,j); + double h_ij1 = h(i,j+1); + unew(i+1,j) = uold(i+1,j) + tdts8 * (z_i1j1 + z_i1j) * (cv_i1j1 + cv_ij1 + cv_ij + cv_i1j) - tdtsdx * (h_i1j - h_ij); + vnew(i,j+1) = vold(i,j+1) - tdts8 * (z_i1j1 + z_ij1) * (cu_i1j1 + cu_ij1 + cu_ij + cu_i1j) - tdtsdy * (h_ij1 - h_ij); + pnew(i,j) = pold(i,j) - tdtsdx * (cu_i1j - cu_ij) - tdtsdy * (cv_ij1 - cv_ij); }); // Periodic continuation @@ -271,41 +274,50 @@ int main(int argc, char **argv) { pnew(M,N) = pnew(0,0); }); + Kokkos::fence(); + t200 += timer200.seconds(); + nvtxRangePop(); + time = time + dt; // Time smoothing and update for next cycle + nvtxRangePush("UpdateOldVariables"); + Kokkos::Timer timer300; if ( ncycle > 1 ) { - Kokkos::parallel_for("time_smoothing_u", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - uold(i,j) = u(i,j) + alpha * (unew(i,j) - 2. * u(i,j) + uold(i,j)); - }); - - Kokkos::parallel_for("time_smoothing_v", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - vold(i,j) = v(i,j) + alpha * (vnew(i,j) - 2. * v(i,j) + vold(i,j)); - }); - - Kokkos::parallel_for("time_smoothing_p", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { - pold(i,j) = p(i,j) + alpha * (pnew(i,j) - 2. * p(i,j) + pold(i,j)); + Kokkos::parallel_for("time_smoothing_uvp", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { + double u_ij = u(i,j); + double v_ij = v(i,j); + double p_ij = p(i,j); + double uold_ij = uold(i,j); + double vold_ij = vold(i,j); + double pold_ij = pold(i,j); + uold(i,j) = u_ij + alpha * (unew(i,j) - 2. * u_ij + uold_ij); + vold(i,j) = v_ij + alpha * (vnew(i,j) - 2. * v_ij + vold_ij); + pold(i,j) = p_ij + alpha * (pnew(i,j) - 2. * p_ij + pold_ij); }); } else { tdt = tdt + tdt; - Kokkos::parallel_for("first_cycle_copy", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}), KOKKOS_LAMBDA(const int i, const int j) { + Kokkos::parallel_for("first_cycle_copy", Kokkos::MDRangePolicy>({0,0}, {M_LEN,N_LEN}, {1,256}), KOKKOS_LAMBDA(const int i, const int j) { uold(i,j) = u(i,j); vold(i,j) = v(i,j); pold(i,j) = p(i,j); }); } - // My test shows it is better to use fence and then swap for device views, rather than doing a device copy without fence. - if constexpr (!std::is_same_v) { + Kokkos::fence(); - } + t300 += timer300.seconds(); + nvtxRangePop(); + // Swap the views std::swap(u, unew); std::swap(v, vnew); std::swap(p, pnew); } // ** End of time loop ** + ctime = timer.seconds(); + // Try to use `if constexpr (!std::is_same_v)` to use swap function // for the host space, but somehow it is not compiled correctly for the device space. // Just use deep_copy for both spaces. @@ -327,23 +339,15 @@ int main(int argc, char **argv) { ptime = time / 3600.; printf(" cycle number %d model time in hours %f\n", ITMAX, ptime); - // printf(" diagonal elements of p\n"); - // for (int i=0; i #include #include +#if defined(__has_include) +#if __has_include() +#include +#define SWM_HAVE_NVTX 1 +#endif +#endif +#ifndef SWM_HAVE_NVTX +#define nvtxRangePush(name) ((void)0) +#define nvtxRangePop() ((void)0) +#endif #define MIN(x,y) ((x)>(y)?(y):(x)) #define MAX(x,y) ((x)>(y)?(x):(y)) @@ -20,7 +30,7 @@ using Layout = Kokkos::LayoutRight; using ExecSpace = Kokkos::DefaultExecutionSpace; using MemSpace = ExecSpace::memory_space; -using ViewMatrixType = Kokkos::View; +using ViewMatrixType = Kokkos::View>; using HostViewMatrixType = Kokkos::View; void write_to_file(auto array, int tM, int tN, const char *filename); @@ -83,10 +93,11 @@ int main(int argc, char **argv) { double tpi,di,dj,pcf; double tdts8,tdtsdx,tdtsdy,fsdx,fsdy; - int mnmin,ncycle; + int ncycle; // timer variables double ctime,tcyc,time,ptime; + double t100 = 0., t200 = 0., t300 = 0.; // ** Initialisations ** @@ -157,26 +168,6 @@ int main(int argc, char **argv) { printf(" grid spacing in the y direction %f\n", dy); printf(" time step %f\n", dt); printf(" time filter parameter %f\n", alpha); - - mnmin = MIN(M,N); - if constexpr (!std::is_same_v) { - Kokkos::deep_copy(u_host, u); - Kokkos::deep_copy(v_host, v); - Kokkos::deep_copy(p_host, p); - } - // printf(" initial diagonal elements of p\n"); - // for (int i=0; i(1, M_LEN), KOKKOS_LAMBDA(const int i) { + // Compute cu, cv, z, h (fused) + Kokkos::parallel_for("compute_cu_cv_z_h", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { for (int j = 0; j < N; ++j) { - cu(i,j) = 0.5 * (p(i,j) + p(i-1,j)) * u(i,j); - } - }); - - // Compute cv - Kokkos::parallel_for("compute_cv", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { - for (int j = 1; j < N_LEN; ++j) { - cv(i,j) = 0.5 * (p(i,j) + p(i,j-1)) * v(i,j); - } - }); - - // Compute z - Kokkos::parallel_for("compute_z", Kokkos::RangePolicy(1, M_LEN), KOKKOS_LAMBDA(const int i) { - for (int j = 1; j < N_LEN; ++j) { - z(i,j) = (fsdx * (v(i,j) - v(i-1,j)) - fsdy * (u(i,j) - u(i,j-1))) / (p(i-1,j-1) + p(i,j-1) + p(i,j) + p(i-1,j)); - } - }); - - // Compute h - Kokkos::parallel_for("compute_h", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { - for (int j = 0; j < N; ++j) { - h(i,j) = p(i,j) + 0.25 * ( u(i+1,j) * u(i+1,j) + u(i,j) * u(i,j) + v(i,j+1) * v(i,j+1) + v(i,j) * v(i,j) ); + double p_ij = p(i,j); + double p_i1j = p(i+1,j); + double p_ij1 = p(i,j+1); + double p_i1j1 = p(i+1,j+1); + double u_ij = u(i,j); + double u_i1j = u(i+1,j); + double u_i1j1 = u(i+1,j+1); + double v_ij = v(i,j); + double v_ij1 = v(i,j+1); + double v_i1j1 = v(i+1,j+1); + cu(i+1,j) = 0.5 * (p_i1j + p_ij) * u_i1j; + cv(i,j+1) = 0.5 * (p_ij1 + p_ij) * v_ij1; + z(i+1,j+1) = (fsdx * (v_i1j1 - v_ij1) - fsdy * (u_i1j1 - u_i1j)) / (p_ij + p_i1j + p_i1j1 + p_ij1); + h(i,j) = p_ij + 0.25 * (u_i1j * u_i1j + u_ij * u_ij + v_ij1 * v_ij1 + v_ij * v_ij); } }); @@ -239,30 +224,38 @@ int main(int argc, char **argv) { h(M,N) = h(0,0); }); + Kokkos::fence(); + t100 += timer100.seconds(); + nvtxRangePop(); + // Compute new values u,v and p + nvtxRangePush("UpdateNewVariables"); + Kokkos::Timer timer200; tdts8 = tdt / 8.; tdtsdx = tdt / dx; tdtsdy = tdt / dy; - // Compute unew - Kokkos::parallel_for("compute_unew", Kokkos::RangePolicy(1, M_LEN), KOKKOS_LAMBDA(const int i) { + // Compute unew, vnew, pnew (fused) + Kokkos::parallel_for("compute_unew_vnew_pnew", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { for (int j = 0; j < N; ++j) { - unew(i,j) = uold(i,j) + tdts8 * (z(i,j+1) + z(i,j)) * (cv(i,j+1) + cv(i-1,j+1) + cv(i-1,j) + cv(i,j)) - tdtsdx * (h(i,j) - h(i-1,j)); - } - }); - - // Compute vnew - Kokkos::parallel_for("compute_vnew", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { - for (int j = 1; j < N_LEN; ++j) { - vnew(i,j) = vold(i,j) - tdts8 * (z(i+1,j) + z(i,j)) * (cu(i+1,j) + cu(i,j) + cu(i,j-1) + cu(i+1,j-1)) - tdtsdy * (h(i,j) - h(i,j-1)); - } - }); - - // Compute pnew - Kokkos::parallel_for("compute_pnew", Kokkos::RangePolicy(0, M), KOKKOS_LAMBDA(const int i) { - for (int j = 0; j < N; ++j) { - pnew(i,j) = pold(i,j) - tdtsdx * (cu(i+1,j) - cu(i,j)) - tdtsdy * (cv(i,j+1) - cv(i,j)); + double z_i1j1 = z(i+1,j+1); + double z_i1j = z(i+1,j); + double z_ij1 = z(i,j+1); + double cv_i1j1 = cv(i+1,j+1); + double cv_ij1 = cv(i,j+1); + double cv_ij = cv(i,j); + double cv_i1j = cv(i+1,j); + double cu_i1j1 = cu(i+1,j+1); + double cu_ij1 = cu(i,j+1); + double cu_ij = cu(i,j); + double cu_i1j = cu(i+1,j); + double h_i1j = h(i+1,j); + double h_ij = h(i,j); + double h_ij1 = h(i,j+1); + unew(i+1,j) = uold(i+1,j) + tdts8 * (z_i1j1 + z_i1j) * (cv_i1j1 + cv_ij1 + cv_ij + cv_i1j) - tdtsdx * (h_i1j - h_ij); + vnew(i,j+1) = vold(i,j+1) - tdts8 * (z_i1j1 + z_ij1) * (cu_i1j1 + cu_ij1 + cu_ij + cu_i1j) - tdtsdy * (h_ij1 - h_ij); + pnew(i,j) = pold(i,j) - tdtsdx * (cu_i1j - cu_ij) - tdtsdy * (cv_ij1 - cv_ij); } }); @@ -286,26 +279,28 @@ int main(int argc, char **argv) { pnew(M,N) = pnew(0,0); }); + Kokkos::fence(); + t200 += timer200.seconds(); + nvtxRangePop(); + time = time + dt; // Time smoothing and update for next cycle + nvtxRangePush("UpdateOldVariables"); + Kokkos::Timer timer300; if ( ncycle > 1 ) { - Kokkos::parallel_for("time_smoothing_u", Kokkos::RangePolicy(0,M_LEN), KOKKOS_LAMBDA(const int i) { - for (int j = 0; j < N_LEN; ++j) { - uold(i,j) = u(i,j) + alpha * (unew(i,j) - 2. * u(i,j) + uold(i,j)); - } - }); - - Kokkos::parallel_for("time_smoothing_v", Kokkos::RangePolicy(0,M_LEN), KOKKOS_LAMBDA(const int i) { + Kokkos::parallel_for("time_smoothing_uvp", Kokkos::RangePolicy(0,M_LEN), KOKKOS_LAMBDA(const int i) { for (int j = 0; j < N_LEN; ++j) { - vold(i,j) = v(i,j) + alpha * (vnew(i,j) - 2. * v(i,j) + vold(i,j)); - } - }); - - Kokkos::parallel_for("time_smoothing_p", Kokkos::RangePolicy(0,M_LEN), KOKKOS_LAMBDA(const int i) { - for (int j = 0; j < N_LEN; ++j) { - pold(i,j) = p(i,j) + alpha * (pnew(i,j) - 2. * p(i,j) + pold(i,j)); + double u_ij = u(i,j); + double v_ij = v(i,j); + double p_ij = p(i,j); + double uold_ij = uold(i,j); + double vold_ij = vold(i,j); + double pold_ij = pold(i,j); + uold(i,j) = u_ij + alpha * (unew(i,j) - 2. * u_ij + uold_ij); + vold(i,j) = v_ij + alpha * (vnew(i,j) - 2. * v_ij + vold_ij); + pold(i,j) = p_ij + alpha * (pnew(i,j) - 2. * p_ij + pold_ij); } }); } @@ -317,16 +312,19 @@ int main(int argc, char **argv) { pold(i,j) = p(i,j); }); } - // My test shows it is better to use fence and then swap for device views, rather than doing a device copy without fence. - if constexpr (!std::is_same_v) { + Kokkos::fence(); - } + t300 += timer300.seconds(); + nvtxRangePop(); + // Swap the views std::swap(u, unew); std::swap(v, vnew); std::swap(p, pnew); } // ** End of time loop ** + ctime = timer.seconds(); + // Try to use `if constexpr (!std::is_same_v)` to use swap function // for the host space, but somehow it is not compiled correctly for the device space. // Just use deep_copy for both spaces. @@ -348,23 +346,15 @@ int main(int argc, char **argv) { ptime = time / 3600.; printf(" cycle number %d model time in hours %f\n", ITMAX, ptime); - // printf(" diagonal elements of p\n"); - // for (int i=0; i #include #include +#if defined(__has_include) +#if __has_include() +#include +#define SWM_HAVE_NVTX 1 +#endif +#endif +#ifndef SWM_HAVE_NVTX +#define nvtxRangePush(name) ((void)0) +#define nvtxRangePop() ((void)0) +#endif #define MIN(x,y) ((x)>(y)?(y):(x)) #define MAX(x,y) ((x)>(y)?(x):(y)) @@ -20,7 +30,7 @@ using Layout = Kokkos::LayoutRight; using ExecSpace = Kokkos::DefaultExecutionSpace; using MemSpace = ExecSpace::memory_space; -using ViewMatrixType = Kokkos::View; +using ViewMatrixType = Kokkos::View>; using HostViewMatrixType = Kokkos::View; void write_to_file(auto array, int tM, int tN, const char *filename); @@ -83,10 +93,11 @@ int main(int argc, char **argv) { double tpi,di,dj,pcf; double tdts8,tdtsdx,tdtsdy,fsdx,fsdy; - int mnmin,ncycle; + int ncycle; // timer variables double ctime,tcyc,time,ptime; + double t100 = 0., t200 = 0., t300 = 0.; // ** Initialisations ** @@ -157,26 +168,6 @@ int main(int argc, char **argv) { printf(" grid spacing in the y direction %f\n", dy); printf(" time step %f\n", dt); printf(" time filter parameter %f\n", alpha); - - mnmin = MIN(M,N); - if constexpr (!std::is_same_v) { - Kokkos::deep_copy(u_host, u); - Kokkos::deep_copy(v_host, v); - Kokkos::deep_copy(p_host, p); - } - // printf(" initial diagonal elements of p\n"); - // for (int i=0; i(0, M * N), KOKKOS_LAMBDA(const int idx) { - int i = idx / N + 1; - int j = idx % N; - cu(i,j) = 0.5 * (p(i,j) + p(i-1,j)) * u(i,j); - }); - - // Compute cv - Kokkos::parallel_for("compute_cv", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { - int i = idx / N; - int j = idx % N + 1; - cv(i,j) = 0.5 * (p(i,j) + p(i,j-1)) * v(i,j); - }); - - // Compute z - Kokkos::parallel_for("compute_z", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { - int i = idx / N + 1; - int j = idx % N + 1; - z(i,j) = (fsdx * (v(i,j) - v(i-1,j)) - fsdy * (u(i,j) - u(i,j-1))) / (p(i-1,j-1) + p(i,j-1) + p(i,j) + p(i-1,j)); - }); - - // Compute h - Kokkos::parallel_for("compute_h", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { + // Compute cu, cv, z, h (fused like OpenACC) + Kokkos::parallel_for("compute_cu_cv_z_h", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { int i = idx / N; int j = idx % N; - h(i,j) = p(i,j) + 0.25 * ( u(i+1,j) * u(i+1,j) + u(i,j) * u(i,j) + v(i,j+1) * v(i,j+1) + v(i,j) * v(i,j) ); + double p_ij = p(i,j); + double p_i1j = p(i+1,j); + double p_ij1 = p(i,j+1); + double p_i1j1 = p(i+1,j+1); + double u_ij = u(i,j); + double u_i1j = u(i+1,j); + double u_i1j1 = u(i+1,j+1); + double v_ij = v(i,j); + double v_ij1 = v(i,j+1); + double v_i1j1 = v(i+1,j+1); + cu(i+1,j) = 0.5 * (p_i1j + p_ij) * u_i1j; + cv(i,j+1) = 0.5 * (p_ij1 + p_ij) * v_ij1; + z(i+1,j+1) = (fsdx * (v_i1j1 - v_ij1) - fsdy * (u_i1j1 - u_i1j)) / (p_ij + p_i1j + p_i1j1 + p_ij1); + h(i,j) = p_ij + 0.25 * (u_i1j * u_i1j + u_ij * u_ij + v_ij1 * v_ij1 + v_ij * v_ij); }); // Periodic continuation @@ -239,30 +224,39 @@ int main(int argc, char **argv) { h(M,N) = h(0,0); }); + Kokkos::fence(); + t100 += timer100.seconds(); + nvtxRangePop(); + // Compute new values u,v and p + nvtxRangePush("UpdateNewVariables"); + Kokkos::Timer timer200; + tdts8 = tdt / 8.; tdtsdx = tdt / dx; tdtsdy = tdt / dy; - // Compute unew - Kokkos::parallel_for("compute_unew", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { - int i = idx / N + 1; - int j = idx % N; - unew(i,j) = uold(i,j) + tdts8 * (z(i,j+1) + z(i,j)) * (cv(i,j+1) + cv(i-1,j+1) + cv(i-1,j) + cv(i,j)) - tdtsdx * (h(i,j) - h(i-1,j)); - }); - - // Compute vnew - Kokkos::parallel_for("compute_vnew", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { - int i = idx / N; - int j = (idx % N) + 1; - vnew(i,j) = vold(i,j) - tdts8 * (z(i+1,j) + z(i,j)) * (cu(i+1,j) + cu(i,j) + cu(i,j-1) + cu(i+1,j-1)) - tdtsdy * (h(i,j) - h(i,j-1)); - }); - - // Compute pnew - Kokkos::parallel_for("compute_pnew", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { + // Compute unew, vnew, pnew (fused like OpenACC) + Kokkos::parallel_for("compute_unew_vnew_pnew", Kokkos::RangePolicy(0, M * N), KOKKOS_LAMBDA(const int idx) { int i = idx / N; int j = idx % N; - pnew(i,j) = pold(i,j) - tdtsdx * (cu(i+1,j) - cu(i,j)) - tdtsdy * (cv(i,j+1) - cv(i,j)); + double z_i1j1 = z(i+1,j+1); + double z_i1j = z(i+1,j); + double z_ij1 = z(i,j+1); + double cv_i1j1 = cv(i+1,j+1); + double cv_ij1 = cv(i,j+1); + double cv_ij = cv(i,j); + double cv_i1j = cv(i+1,j); + double cu_i1j1 = cu(i+1,j+1); + double cu_ij1 = cu(i,j+1); + double cu_ij = cu(i,j); + double cu_i1j = cu(i+1,j); + double h_i1j = h(i+1,j); + double h_ij = h(i,j); + double h_ij1 = h(i,j+1); + unew(i+1,j) = uold(i+1,j) + tdts8 * (z_i1j1 + z_i1j) * (cv_i1j1 + cv_ij1 + cv_ij + cv_i1j) - tdtsdx * (h_i1j - h_ij); + vnew(i,j+1) = vold(i,j+1) - tdts8 * (z_i1j1 + z_ij1) * (cu_i1j1 + cu_ij1 + cu_ij + cu_i1j) - tdtsdy * (h_ij1 - h_ij); + pnew(i,j) = pold(i,j) - tdtsdx * (cu_i1j - cu_ij) - tdtsdy * (cv_ij1 - cv_ij); }); // Periodic continuation @@ -284,26 +278,29 @@ int main(int argc, char **argv) { pnew(M,N) = pnew(0,0); }); + Kokkos::fence(); + t200 += timer200.seconds(); + nvtxRangePop(); + time = time + dt; // Time smoothing and update for next cycle - if ( ncycle > 1 ) { - Kokkos::parallel_for("time_smoothing_u", Kokkos::RangePolicy(0, M_LEN * N_LEN), KOKKOS_LAMBDA(const int idx) { - int i = idx / N_LEN; - int j = idx % N_LEN; - uold(i,j) = u(i,j) + alpha * (unew(i,j) - 2. * u(i,j) + uold(i,j)); - }); - - Kokkos::parallel_for("time_smoothing_v", Kokkos::RangePolicy(0, M_LEN * N_LEN), KOKKOS_LAMBDA(const int idx) { - int i = idx / N_LEN; - int j = idx % N_LEN; - vold(i,j) = v(i,j) + alpha * (vnew(i,j) - 2. * v(i,j) + vold(i,j)); - }); + nvtxRangePush("UpdateOldVariables"); + Kokkos::Timer timer300; - Kokkos::parallel_for("time_smoothing_p", Kokkos::RangePolicy(0, M_LEN * N_LEN), KOKKOS_LAMBDA(const int idx) { + if ( ncycle > 1 ) { + Kokkos::parallel_for("time_smoothing_uvp", Kokkos::RangePolicy(0, M_LEN * N_LEN), KOKKOS_LAMBDA(const int idx) { int i = idx / N_LEN; int j = idx % N_LEN; - pold(i,j) = p(i,j) + alpha * (pnew(i,j) - 2. * p(i,j) + pold(i,j)); + double u_ij = u(i,j); + double v_ij = v(i,j); + double p_ij = p(i,j); + double uold_ij = uold(i,j); + double vold_ij = vold(i,j); + double pold_ij = pold(i,j); + uold(i,j) = u_ij + alpha * (unew(i,j) - 2. * u_ij + uold_ij); + vold(i,j) = v_ij + alpha * (vnew(i,j) - 2. * v_ij + vold_ij); + pold(i,j) = p_ij + alpha * (pnew(i,j) - 2. * p_ij + pold_ij); }); } else { @@ -316,16 +313,19 @@ int main(int argc, char **argv) { pold(i,j) = p(i,j); }); } - // My test shows it is better to use fence and then swap for device views, rather than doing a device copy without fence. - if constexpr (!std::is_same_v) { + Kokkos::fence(); - } + t300 += timer300.seconds(); + nvtxRangePop(); + // Swap the views std::swap(u, unew); std::swap(v, vnew); std::swap(p, pnew); } // ** End of time loop ** + ctime = timer.seconds(); + // Try to use `if constexpr (!std::is_same_v)` to use swap function // for the host space, but somehow it is not compiled correctly for the device space. // Just use deep_copy for both spaces. @@ -347,23 +347,15 @@ int main(int argc, char **argv) { ptime = time / 3600.; printf(" cycle number %d model time in hours %f\n", ITMAX, ptime); - // printf(" diagonal elements of p\n"); - // for (int i=0; i #include #include +#if defined(__has_include) +#if __has_include() +#include +#define SWM_HAVE_NVTX 1 +#endif +#endif +#ifndef SWM_HAVE_NVTX +#define nvtxRangePush(name) ((void)0) +#define nvtxRangePop() ((void)0) +#endif #define MIN(x,y) ((x)>(y)?(y):(x)) #define MAX(x,y) ((x)>(y)?(x):(y)) @@ -20,7 +30,7 @@ using Layout = Kokkos::LayoutRight; using ExecSpace = Kokkos::DefaultExecutionSpace; using MemSpace = ExecSpace::memory_space; -using ViewMatrixType = Kokkos::View; +using ViewMatrixType = Kokkos::View>; using HostViewMatrixType = Kokkos::View; using teamPolicy = Kokkos::TeamPolicy; using memberType = teamPolicy::member_type; @@ -85,10 +95,11 @@ int main(int argc, char **argv) { double tpi,di,dj,pcf; double tdts8,tdtsdx,tdtsdy,fsdx,fsdy; - int mnmin,ncycle; + int ncycle; // timer variables double ctime,tcyc,time,ptime; + double t100 = 0., t200 = 0., t300 = 0.; // ** Initialisations ** @@ -170,26 +181,6 @@ int main(int argc, char **argv) { printf(" grid spacing in the y direction %f\n", dy); printf(" time step %f\n", dt); printf(" time filter parameter %f\n", alpha); - - mnmin = MIN(M,N); - if constexpr (!std::is_same_v) { - Kokkos::deep_copy(u_host, u); - Kokkos::deep_copy(v_host, v); - Kokkos::deep_copy(p_host, p); - } - // printf(" initial diagonal elements of p\n"); - // for (int i=0; i 1 ) { - Kokkos::parallel_for("time_smoothing_u", teamPolicy(M_LEN, Kokkos::AUTO), KOKKOS_LAMBDA(const memberType& team) { - const int i = team.league_rank(); - Kokkos::parallel_for(Kokkos::TeamThreadRange(team, N_LEN), [&](const int j) { - uold(i,j) = u(i,j) + alpha * (unew(i,j) - 2. * u(i,j) + uold(i,j)); - }); - }); - - Kokkos::parallel_for("time_smoothing_v", teamPolicy(M_LEN, Kokkos::AUTO), KOKKOS_LAMBDA(const memberType& team) { + Kokkos::parallel_for("time_smoothing_uvp", teamPolicy(M_LEN, Kokkos::AUTO), KOKKOS_LAMBDA(const memberType& team) { const int i = team.league_rank(); Kokkos::parallel_for(Kokkos::TeamThreadRange(team, N_LEN), [&](const int j) { - vold(i,j) = v(i,j) + alpha * (vnew(i,j) - 2. * v(i,j) + vold(i,j)); - }); - }); - - Kokkos::parallel_for("time_smoothing_p", teamPolicy(M_LEN, Kokkos::AUTO), KOKKOS_LAMBDA(const memberType& team) { - const int i = team.league_rank(); - Kokkos::parallel_for(Kokkos::TeamThreadRange(team, N_LEN), [&](const int j) { - pold(i,j) = p(i,j) + alpha * (pnew(i,j) - 2. * p(i,j) + pold(i,j)); + double u_ij = u(i,j); + double v_ij = v(i,j); + double p_ij = p(i,j); + double uold_ij = uold(i,j); + double vold_ij = vold(i,j); + double pold_ij = pold(i,j); + uold(i,j) = u_ij + alpha * (unew(i,j) - 2. * u_ij + uold_ij); + vold(i,j) = v_ij + alpha * (vnew(i,j) - 2. * v_ij + vold_ij); + pold(i,j) = p_ij + alpha * (pnew(i,j) - 2. * p_ij + pold_ij); }); }); } @@ -346,16 +334,19 @@ int main(int argc, char **argv) { }); }); } - // My test shows it is better to use fence and then swap for device views, rather than doing a device copy without fence. - if constexpr (!std::is_same_v) { + Kokkos::fence(); - } + t300 += timer300.seconds(); + nvtxRangePop(); + // Swap the views std::swap(u, unew); std::swap(v, vnew); std::swap(p, pnew); } // ** End of time loop ** + ctime = timer.seconds(); + // Try to use `if constexpr (!std::is_same_v)` to use swap function // for the host space, but somehow it is not compiled correctly for the device space. // Just use deep_copy for both spaces. @@ -377,23 +368,15 @@ int main(int argc, char **argv) { ptime = time / 3600.; printf(" cycle number %d model time in hours %f\n", ITMAX, ptime); - // printf(" diagonal elements of p\n"); - // for (int i=0; i