feat(sphere): add OpenMP CPU backend sharing kernel sources with CUDA - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit 07cbc0187d7275aa422010407afa0f346d005309
DIR parent a59d99c34192952f05e3485f16adce80bcc2aa63
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Sat, 4 Jul 2026 23:46:42 +0200
feat(sphere): add OpenMP CPU backend sharing kernel sources with CUDA
Diffstat:
M CMakeLists.txt | 46 +++++++++++++++++++------------
M src/CMakeLists.txt | 105 +++++++++++++++++++------------
M src/darcy.cuh | 2 ++
M src/datatypes.h | 2 ++
M src/device.cu | 491 ++++++++++++++++---------------
M src/integration.cuh | 20 ++++++++++++++++++++
M src/navierstokes.cuh | 4 +++-
M src/raytracer.cuh | 21 +++++++++++----------
M src/sorting.cuh | 40 ++++++++++++++++++++++++++++++-
M src/typedefs.h | 4 ++++
M src/vector_arithmetic.h | 2 ++
M tests/CMakeLists.txt | 13 ++++++++++---
12 files changed, 441 insertions(+), 309 deletions(-)
---
DIR diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,26 +1,38 @@
+# CMake minimum version required
+cmake_minimum_required(VERSION 3.12)
+
+# The name of the project.
+project(sphere_CUDA)
+
# Create input/output folders
file(MAKE_DIRECTORY input)
file(MAKE_DIRECTORY output)
file(MAKE_DIRECTORY img_out)
file(MAKE_DIRECTORY gnuplot/data)
-# The name of the project.
-project(sphere_CUDA)
-
-# CMake minimum version required
-cmake_minimum_required(VERSION 3.0)
-
-# Find CUDA
-find_package(CUDA REQUIRED)
-
-# Set CUDA GPU generation
-# 0: Fermi
-# 1: Kepler
-# 2: Ampere
-set(GPU_GENERATION 2)
-
-# Find OpenMP
-find_package(OpenMP)
+# Backend selection: CUDA (default) or OpenMP CPU
+option(SPHERE_GPU "Build the CUDA backend (requires CUDA toolkit)" ON)
+
+if (SPHERE_GPU)
+ # Find CUDA
+ find_package(CUDA REQUIRED)
+
+ # Set CUDA GPU generation
+ # 0: Fermi
+ # 1: Kepler
+ # 2: Ampere
+ set(GPU_GENERATION 2)
+
+ # Find OpenMP
+ find_package(OpenMP)
+else()
+ # On macOS with Apple clang, OpenMP comes from Homebrew libomp
+ if (APPLE AND NOT DEFINED OpenMP_ROOT
+ AND EXISTS /opt/homebrew/opt/libomp)
+ set(OpenMP_ROOT /opt/homebrew/opt/libomp)
+ endif()
+ find_package(OpenMP REQUIRED COMPONENTS CXX)
+endif()
# Find Boost components
#find_package(Boost COMPONENTS system filesystem unit_test_framework REQUIRED)
DIR diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -1,48 +1,71 @@
-# Link with libcutil (expected to be in ${CUDA_SDK_ROOT_DIR}/C/lib)
-#LINK_LIBRARIES("-L${CUDA_SDK_ROOT_DIR}/lib -lcutil") # For 32 bit systems
-#LINK_LIBRARIES("-L${CUDA_SDK_ROOT_DIR}/lib -lcutil_x86_64") # For 64 bit systems
-
-# Ohter folders to include
-#SET(CUDA_SDK_ROOT_DIR "/usr/local/cuda-5.0/samples")
-INCLUDE_DIRECTORIES("${CUDA_SDK_ROOT_DIR}/common/inc")
INCLUDE_DIRECTORIES("${CMAKE_BINARY_DIR}/src")
SET(EXECUTABLE_OUTPUT_PATH "../")
-# Include FindCUDA script
-INCLUDE(FindCUDA)
-
-# Additional NVCC command line arguments
-# NOTE: Multiple arguments must be semi-colon selimited
-IF (GPU_GENERATION EQUAL 2) # Ampere
- SET(CUDA_NVCC_FLAGS
- "--use_fast_math;-O3;-gencode=arch=compute_80,code=\"sm_80,compute_80\";--fmad=false -ccbin gcc")
-ELSEIF (GPU_GENERATION EQUAL 1) # Kepler
- SET(CUDA_NVCC_FLAGS
- "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc")
- #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-8.3")
- #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-4.6")
- #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
- #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false;-ccbin clang-3.8")
-ELSE() # Fermi
- SET(CUDA_NVCC_FLAGS
- "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc")
- #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc-8.3")
- #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
- #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false;-ccbin clang-3.8")
-ENDIF ()
+IF (SPHERE_GPU)
+
+ # Link with libcutil (expected to be in ${CUDA_SDK_ROOT_DIR}/C/lib)
+ #LINK_LIBRARIES("-L${CUDA_SDK_ROOT_DIR}/lib -lcutil") # For 32 bit systems
+ #LINK_LIBRARIES("-L${CUDA_SDK_ROOT_DIR}/lib -lcutil_x86_64") # For 64 bit systems
+
+ # Ohter folders to include
+ #SET(CUDA_SDK_ROOT_DIR "/usr/local/cuda-5.0/samples")
+ INCLUDE_DIRECTORIES("${CUDA_SDK_ROOT_DIR}/common/inc")
+
+ # Include FindCUDA script
+ INCLUDE(FindCUDA)
+
+ # Additional NVCC command line arguments
+ # NOTE: Multiple arguments must be semi-colon selimited
+ IF (GPU_GENERATION EQUAL 2) # Ampere
+ SET(CUDA_NVCC_FLAGS
+ "--use_fast_math;-O3;-gencode=arch=compute_80,code=\"sm_80,compute_80\";--fmad=false -ccbin gcc")
+ ELSEIF (GPU_GENERATION EQUAL 1) # Kepler
+ SET(CUDA_NVCC_FLAGS
+ "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc")
+ ELSE() # Fermi
+ SET(CUDA_NVCC_FLAGS
+ "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc")
+ ENDIF ()
+
+ # Backend selection, visible to both nvcc and the host compiler
+ SET(CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS};-DSPHERE_GPU")
+ ADD_DEFINITIONS(-DSPHERE_GPU)
+
+ SET(CMAKE_CXX_FLAGS "-fPIC ${CMAKE_CXX_FLAGS}")
-SET(CMAKE_CXX_FLAGS "-fPIC ${CMAKE_CXX_FLAGS}")
-
-# Rule to build executable program
-CUDA_ADD_EXECUTABLE(sphere
- main.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
- navierstokes.cpp darcy.cpp)
-CUDA_ADD_EXECUTABLE(porosity
- porosity.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
- navierstokes.cpp darcy.cpp)
-CUDA_ADD_EXECUTABLE(forcechains
- forcechains.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
- navierstokes.cpp darcy.cpp)
+ # Rule to build executable program
+ CUDA_ADD_EXECUTABLE(sphere
+ main.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
+ navierstokes.cpp darcy.cpp)
+ CUDA_ADD_EXECUTABLE(porosity
+ porosity.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
+ navierstokes.cpp darcy.cpp)
+ CUDA_ADD_EXECUTABLE(forcechains
+ forcechains.cpp file_io.cpp sphere.cpp device.cu utility.cu utility.cpp
+ navierstokes.cpp darcy.cpp)
+
+ELSE()
+
+ SET(CMAKE_CXX_STANDARD 11)
+ SET(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+ # device_cpu.cpp compiles device.cu as plain C++ through the shims in
+ # cpu_backend.h; cpu_backend.cpp replaces utility.cu
+ SET(SPHERE_CPU_SRC
+ file_io.cpp sphere.cpp device_cpu.cpp cpu_backend.cpp utility.cpp
+ navierstokes.cpp darcy.cpp)
+
+ ADD_EXECUTABLE(sphere main.cpp ${SPHERE_CPU_SRC})
+ ADD_EXECUTABLE(porosity porosity.cpp ${SPHERE_CPU_SRC})
+ ADD_EXECUTABLE(forcechains forcechains.cpp ${SPHERE_CPU_SRC})
+
+ FOREACH(TARGET sphere porosity forcechains)
+ TARGET_LINK_LIBRARIES(${TARGET} OpenMP::OpenMP_CXX)
+ # -ffp-contract=off mirrors the nvcc --fmad=false choice
+ TARGET_COMPILE_OPTIONS(${TARGET} PRIVATE -O3 -ffp-contract=off)
+ ENDFOREACH()
+
+ENDIF ()
#ADD_EXECUTABLE(unittests boost-unit-tests.cpp sphere.cpp)
#TARGET_LINK_LIBRARIES(unittests
DIR diff --git a/src/darcy.cuh b/src/darcy.cuh
@@ -2,9 +2,11 @@
// CUDA implementation of Darcy porous flow
#include <iostream>
+#ifdef SPHERE_GPU
#include <cuda.h>
//#include <cutil_math.h>
#include <helper_math.h>
+#endif
#include "vector_arithmetic.h" // for arbitrary precision vectors
#include "sphere.h"
DIR diff --git a/src/datatypes.h b/src/datatypes.h
@@ -3,7 +3,9 @@
#define DATATYPES_H_
#include <math.h>
+#ifdef SPHERE_GPU
#include "vector_functions.h"
+#endif
#include "typedefs.h"
#include "constants.h"
DIR diff --git a/src/device.cu b/src/device.cu
@@ -3,15 +3,24 @@
#include <fstream>
#include <string>
#include <cstdio>
-#include <cuda.h>
-#include <helper_math.h>
#include <iomanip>
#include <time.h>
-#include "vector_arithmetic.h" // for arbitrary prec. vectors
-//#include <vector_functions.h> // for single prec. vectors
+#ifdef SPHERE_GPU
+#include <cuda.h>
+#include <helper_math.h>
#include "thrust/device_ptr.h"
#include "thrust/sort.h"
+#else
+#include <algorithm>
+#include <utility>
+#include <vector>
+#include "cpu_backend.h"
+#endif
+
+#include "vector_arithmetic.h" // for arbitrary prec. vectors
+//#include <vector_functions.h> // for single prec. vectors
+#include "launch.h"
#include "sphere.h"
#include "datatypes.h"
@@ -29,6 +38,7 @@
#include "navierstokes.cuh"
#include "darcy.cuh"
+#ifdef SPHERE_GPU
// Returns the number of cores per streaming multiprocessor, which is
// a function of the device compute capability
// https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities
@@ -155,6 +165,7 @@ void DEM::initializeGPU(void)
checkForCudaErrors("While initializing CUDA device");
}
+#endif // SPHERE_GPU
// Start timer for kernel profiling
void startTimer(cudaEvent_t* kernel_tic)
@@ -306,7 +317,8 @@ void DEM::checkConstantMemory()
*equal = 0;
int* dev_equal;
cudaMalloc((void**)&dev_equal, sizeof(int));
- checkConstantValues<<<1,1>>>(dev_equal, dev_grid, dev_params);
+ KERNEL_LAUNCH(checkConstantValues, 1, 1, 0,
+ (dev_equal, dev_grid, dev_params));
checkForCudaErrors("After constant memory check");
// Copy result to host
@@ -1049,7 +1061,8 @@ void DEM::startTime()
// check if particle positions have finite values
#ifdef CHECK_PARTICLES_FINITE
- checkParticlePositions<<<dimGrid, dimBlock>>>(dev_x);
+ KERNEL_LAUNCH(checkParticlePositions, dimGrid, dimBlock, 0,
+ (dev_x));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post checkParticlePositions", iter);
#endif
@@ -1065,9 +1078,10 @@ void DEM::startTime()
// in the fine, uniform and homogenous grid.
if (PROFILING == 1)
startTimer(&kernel_tic);
- calcParticleCellID<<<dimGrid, dimBlock>>>(dev_gridParticleCellID,
+ KERNEL_LAUNCH(calcParticleCellID, dimGrid, dimBlock, 0,
+ (dev_gridParticleCellID,
dev_gridParticleIndex,
- dev_x);
+ dev_x));
// Synchronization point
cudaDeviceSynchronize();
@@ -1081,10 +1095,27 @@ void DEM::startTime()
// radix sort
if (PROFILING == 1)
startTimer(&kernel_tic);
+#ifdef SPHERE_GPU
thrust::sort_by_key(
thrust::device_ptr<uint>(dev_gridParticleCellID),
thrust::device_ptr<uint>(dev_gridParticleCellID + np),
thrust::device_ptr<uint>(dev_gridParticleIndex));
+#else
+ { // stable sort (key, index) pairs to match thrust radix order
+ std::vector<std::pair<unsigned int, unsigned int> > kv(np);
+ for (unsigned int i = 0; i < np; ++i)
+ kv[i] = std::make_pair(dev_gridParticleCellID[i],
+ dev_gridParticleIndex[i]);
+ std::stable_sort(kv.begin(), kv.end(),
+ [](const std::pair<unsigned int, unsigned int>& a,
+ const std::pair<unsigned int, unsigned int>& b)
+ { return a.first < b.first; });
+ for (unsigned int i = 0; i < np; ++i) {
+ dev_gridParticleCellID[i] = kv[i].first;
+ dev_gridParticleIndex[i] = kv[i].second;
+ }
+ }
+#endif
cudaDeviceSynchronize(); // Maybe Thrust synchronizes implicitly?
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1105,7 +1136,8 @@ void DEM::startTime()
// configurations in new arrays (*_sorted).
if (PROFILING == 1)
startTimer(&kernel_tic);
- reorderArrays<<<dimGrid, dimBlock, smemSize>>>(dev_cellStart,
+ KERNEL_LAUNCH(reorderArrays, dimGrid, dimBlock, smemSize,
+ (dev_cellStart,
dev_cellEnd,
dev_gridParticleCellID,
dev_gridParticleIndex,
@@ -1113,7 +1145,7 @@ void DEM::startTime()
dev_angvel,
dev_x_sorted,
dev_vel_sorted,
- dev_angvel_sorted);
+ dev_angvel_sorted));
// Synchronization point
cudaDeviceSynchronize();
@@ -1129,12 +1161,13 @@ void DEM::startTime()
// For each particle: Search contacts in neighbor cells
if (PROFILING == 1)
startTimer(&kernel_tic);
- topology<<<dimGrid, dimBlock>>>(dev_cellStart,
+ KERNEL_LAUNCH(topology, dimGrid, dimBlock, 0,
+ (dev_cellStart,
dev_cellEnd,
dev_gridParticleIndex,
dev_x_sorted,
dev_contacts,
- dev_distmod);
+ dev_distmod));
// Synchronization point
cudaDeviceSynchronize();
@@ -1152,7 +1185,8 @@ void DEM::startTime()
//cudaPrintfInit();
if (PROFILING == 1)
startTimer(&kernel_tic);
- interact<<<dimGrid, dimBlock>>>(dev_gridParticleIndex,
+ KERNEL_LAUNCH(interact, dimGrid, dimBlock, 0,
+ (dev_gridParticleIndex,
dev_cellStart,
dev_cellEnd,
dev_x,
@@ -1173,7 +1207,7 @@ void DEM::startTime()
dev_walls_force_pp,
dev_contacts,
dev_distmod,
- dev_delta_t);
+ dev_delta_t));
// Synchronization point
cudaDeviceSynchronize();
@@ -1189,15 +1223,15 @@ void DEM::startTime()
if (params.nb0 > 0) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- bondsLinear<<<dimGridBonds, dimBlock>>>(
- dev_bonds,
+ KERNEL_LAUNCH(bondsLinear, dimGridBonds, dimBlock, 0,
+ (dev_bonds,
dev_bonds_delta,
dev_bonds_omega,
dev_x,
dev_vel,
dev_angvel,
dev_force,
- dev_torque);
+ dev_torque));
// Synchronization point
cudaDeviceSynchronize();
//cudaPrintfDisplay(stdout, true);
@@ -1220,10 +1254,9 @@ void DEM::startTime()
// the fluid velocities
if (PROFILING == 1)
startTimer(&kernel_tic);
- findPorositiesVelocitiesDiametersSpherical
//findPorositiesVelocitiesDiametersSphericalGradient
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_cellStart,
+ KERNEL_LAUNCH(findPorositiesVelocitiesDiametersSpherical, dimGridFluid, dimBlockFluid, 0,
+ (dev_cellStart,
dev_cellEnd,
dev_x_sorted,
dev_vel_sorted,
@@ -1233,7 +1266,7 @@ void DEM::startTime()
dev_ns_d_avg,
iter,
np,
- ns.c_phi);
+ ns.c_phi));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1249,59 +1282,56 @@ void DEM::startTime()
}*/
if (iter == 0) {
// set cell center ghost nodes
- setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v, ns.bc_bot, ns.bc_top));
// find cell face velocities
- interpolateCenterToFace
- <<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v,
+ KERNEL_LAUNCH(interpolateCenterToFace, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v,
dev_ns_v_x,
dev_ns_v_y,
- dev_ns_v_z);
+ dev_ns_v_z));
cudaDeviceSynchronize();
checkForCudaErrors("Post interpolateCenterToFace");
}
- setNSghostNodesFace<Float>
- <<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_x,
+ KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
ns.bc_bot,
- ns.bc_top);
+ ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodesFace", iter);
- findFaceDivTau<<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_x,
+ KERNEL_LAUNCH(findFaceDivTau, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
ns.mu,
dev_ns_div_tau_x,
dev_ns_div_tau_y,
- dev_ns_div_tau_z);
+ dev_ns_div_tau_z));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post findFaceDivTau", iter);
- setNSghostNodesFace<Float>
- <<<dimGridFluidFace, dimBlockFluid>>>(
- dev_ns_div_tau_x,
+ KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluid, 0,
+ (dev_ns_div_tau_x,
dev_ns_div_tau_y,
dev_ns_div_tau_z,
ns.bc_bot,
- ns.bc_top);
+ ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_div_tau)",
iter);
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_p, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_p, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)", iter);
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_phi, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_phi, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)", iter);
@@ -1310,8 +1340,8 @@ void DEM::startTime()
// Per particle, find the fluid-particle interaction force
// f_pf and apply it to the particle
- findInteractionForce<<<dimGrid, dimBlock>>>(
- dev_x,
+ KERNEL_LAUNCH(findInteractionForce, dimGrid, dimBlock, 0,
+ (dev_x,
dev_vel,
dev_ns_phi,
dev_ns_p,
@@ -1329,24 +1359,23 @@ void DEM::startTime()
dev_ns_f_d,
dev_ns_f_p,
dev_ns_f_v,
- dev_ns_f_sum);
+ dev_ns_f_sum));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post findInteractionForce", iter);
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_p, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_p, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)",
iter);
// Apply fluid-particle interaction force to the fluid
- applyInteractionForceToFluid
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_gridParticleIndex,
+ KERNEL_LAUNCH(applyInteractionForceToFluid, dimGridFluid, dimBlockFluid, 0,
+ (dev_gridParticleIndex,
dev_cellStart,
dev_cellEnd,
dev_ns_f_pf,
- dev_ns_F_pf);
+ dev_ns_F_pf));
//dev_ns_F_pf_x,
//dev_ns_F_pf_y,
//dev_ns_F_pf_z);
@@ -1354,8 +1383,8 @@ void DEM::startTime()
checkForCudaErrorsIter("Post applyInteractionForceToFluid",
iter);
- setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_F_pf, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_F_pf, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(F_pf)", iter);
}
@@ -1368,10 +1397,10 @@ void DEM::startTime()
Float pressure = ns.p[idx(0,0,ns.nz-1)];
Float pressure_new = pressure; // Dirichlet
Float epsilon_value = pressure_new - ns.beta*pressure;
- setNSepsilonTop<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(setNSepsilonTop, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
dev_ns_epsilon_new,
- epsilon_value);
+ epsilon_value));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSepsilonTop", iter);
@@ -1385,12 +1414,12 @@ void DEM::startTime()
if (walls.nw > 0 &&
(walls.wmode[0] == 1 || walls.wmode[0] == 3)) {
wall0_iz = walls.nx->w/(grid.L[2]/grid.num[2]);
- setNSepsilonAtTopWall<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(setNSepsilonAtTopWall, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
dev_ns_epsilon_new,
wall0_iz,
epsilon_value,
- dp_dz);
+ dp_dz));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSepsilonAtTopWall",
iter);
@@ -1411,12 +1440,12 @@ void DEM::startTime()
Float new_pressure = ns.p[idx(0,0,ns.nz-1)]
+ ns.p_mod_A*sin(2.0*M_PI*ns.p_mod_f*time.current
+ ns.p_mod_phi);
- setUpperPressureNS<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_p,
+ KERNEL_LAUNCH(setUpperPressureNS, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_p,
dev_ns_epsilon,
dev_ns_epsilon_new,
ns.beta,
- new_pressure);
+ new_pressure));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setUpperPressureNS", iter);
@@ -1434,24 +1463,23 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_p, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_p, ns.bc_bot, ns.bc_top));
//setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
//dev_ns_v, ns.bc_bot, ns.bc_top);
- setNSghostNodesFace<Float>
- <<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_p_x,
+ KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z,
- ns.bc_bot, ns.bc_top);
+ ns.bc_bot, ns.bc_top));
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_phi, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_phi, ns.bc_bot, ns.bc_top));
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_dphi, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_dphi, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
@@ -1466,18 +1494,18 @@ void DEM::startTime()
// interpolate velocities to cell centers which makes
// velocity prediction easier
- interpolateFaceToCenter<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v_x,
+ KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
- dev_ns_v);
+ dev_ns_v));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"Post interpolateFaceToCenter", iter);
// Set cell-center velocity ghost nodes
- setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(v)", iter);
@@ -1485,10 +1513,10 @@ void DEM::startTime()
// the fluid velocities
if (PROFILING == 1)
startTimer(&kernel_tic);
- findNSdivphiviv<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_phi,
+ KERNEL_LAUNCH(findNSdivphiviv, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_phi,
dev_ns_v,
- dev_ns_div_phi_vi_v);
+ dev_ns_div_phi_vi_v));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1496,8 +1524,8 @@ void DEM::startTime()
checkForCudaErrorsIter("Post findNSdivphiviv", iter);
// Set cell-center ghost nodes
- setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_div_phi_vi_v, ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_div_phi_vi_v, ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodes(div_phi_vi_v)",
iter);
@@ -1507,8 +1535,8 @@ void DEM::startTime()
// constraint
if (PROFILING == 1)
startTimer(&kernel_tic);
- findPredNSvelocities<<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_p,
+ KERNEL_LAUNCH(findPredNSvelocities, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_p,
dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
@@ -1529,28 +1557,27 @@ void DEM::startTime()
ns.rho_f,
dev_ns_v_p_x,
dev_ns_v_p_y,
- dev_ns_v_p_z);
+ dev_ns_v_p_z));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
&t_findPredNSvelocities);
checkForCudaErrorsIter("Post findPredNSvelocities", iter);
- setNSghostNodesFace<Float>
- <<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_p_x,
+ KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z,
- ns.bc_bot, ns.bc_top);
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"Post setNSghostNodesFace(dev_ns_v_p)", iter);
- interpolateFaceToCenter<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v_p_x,
+ KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z,
- dev_ns_v_p);
+ dev_ns_v_p));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"Post interpolateFaceToCenter", iter);
@@ -1568,12 +1595,12 @@ void DEM::startTime()
Float epsilon_value = pressure_new - ns.beta*pressure;
if (PROFILING == 1)
startTimer(&kernel_tic);
- setNSepsilonInterior<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon, epsilon_value);
+ KERNEL_LAUNCH(setNSepsilonInterior, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon, epsilon_value));
cudaDeviceSynchronize();
- setNSnormZero<<<dimGridFluid, dimBlockFluid>>>
- (dev_ns_norm);
+ KERNEL_LAUNCH(setNSnormZero, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_norm));
cudaDeviceSynchronize();
if (PROFILING == 1)
@@ -1596,10 +1623,10 @@ void DEM::startTime()
epsilon_value = pressure_new - ns.beta*pressure;
if (PROFILING == 1)
startTimer(&kernel_tic);
- setNSepsilonBottom<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(setNSepsilonBottom, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
dev_ns_epsilon_new,
- epsilon_value);
+ epsilon_value));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1621,9 +1648,9 @@ void DEM::startTime()
checkForCudaErrors(
"Post setNSghostNodesFloat(dev_ns_epsilon)",
iter);*/
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
- ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodesEpsilon(1)",
iter);
@@ -1666,8 +1693,8 @@ void DEM::startTime()
// iteration.
if (PROFILING == 1)
startTimer(&kernel_tic);
- findNSforcing<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(findNSforcing, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
dev_ns_phi,
dev_ns_dphi,
dev_ns_v_p,
@@ -1680,7 +1707,7 @@ void DEM::startTime()
ns.rho_f,
dev_ns_f1,
dev_ns_f2,
- dev_ns_f);
+ dev_ns_f));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1696,9 +1723,9 @@ void DEM::startTime()
checkForCudaErrors("Post setNSghostNodesForcing",
iter);*/
- setNSghostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
- ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setNSghostNodesEpsilon(2)",
iter);
@@ -1715,8 +1742,8 @@ void DEM::startTime()
// Perform a single Jacobi iteration
if (PROFILING == 1)
startTimer(&kernel_tic);
- jacobiIterationNS<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(jacobiIterationNS, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
dev_ns_epsilon_new,
dev_ns_norm,
dev_ns_f,
@@ -1724,7 +1751,7 @@ void DEM::startTime()
ns.bc_top,
ns.theta,
wall0_iz,
- dp_dz);
+ dp_dz));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1747,9 +1774,9 @@ void DEM::startTime()
}*/
// Copy new values to current values
- copyValues<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon_new,
- dev_ns_epsilon);
+ KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon_new,
+ dev_ns_epsilon));
cudaDeviceSynchronize();
checkForCudaErrorsIter
("Post copyValues (epsilon_new->epsilon)", iter);
@@ -1787,10 +1814,9 @@ void DEM::startTime()
&& iter % conv_log_interval == 0)
convlog << iter+1 << '\t' << nijac << std::endl;
- setNSghostNodes<Float>
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
- ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter
("Post setNSghostNodesEpsilon(4)", iter);
@@ -1798,17 +1824,16 @@ void DEM::startTime()
// Apply smoothing if requested
if (ns.gamma > 0.0) {
- smoothing<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(smoothing, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
ns.gamma,
- ns.bc_bot, ns.bc_top);
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post smoothing", iter);
- setNSghostNodes<Float>
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
- ns.bc_bot, ns.bc_top);
+ KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter
("Post setNSghostNodesEpsilon(4)", iter);
@@ -1848,15 +1873,15 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- updateNSpressure<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_epsilon,
+ KERNEL_LAUNCH(updateNSpressure, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_epsilon,
ns.beta,
- dev_ns_p);
+ dev_ns_p));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post updateNSpressure", iter);
- updateNSvelocity<<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_p_x,
+ KERNEL_LAUNCH(updateNSvelocity, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z,
dev_ns_phi,
@@ -1871,28 +1896,27 @@ void DEM::startTime()
iter,
dev_ns_v_x,
dev_ns_v_y,
- dev_ns_v_z);
+ dev_ns_v_z));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
&t_updateNSvelocityPressure);
checkForCudaErrorsIter("Post updateNSvelocity", iter);
- setNSghostNodesFace<Float>
- <<<dimGridFluidFace, dimBlockFluidFace>>>(
- dev_ns_v_p_x,
+ KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
+ (dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z,
- ns.bc_bot, ns.bc_top);
+ ns.bc_bot, ns.bc_top));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"Post setNSghostNodesFace(dev_ns_v)", iter);
- interpolateFaceToCenter<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v_x,
+ KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
- dev_ns_v);
+ dev_ns_v));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post interpolateFaceToCenter",
iter);
@@ -1916,11 +1940,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1930,18 +1954,18 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- findDarcyPressureGradient<<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
- dev_darcy_grad_p);
+ KERNEL_LAUNCH(findDarcyPressureGradient, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
+ dev_darcy_grad_p));
cudaDeviceSynchronize();
checkForCudaErrorsIter("After findDarcyPressureGradient",
iter);
- setDarcyGhostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_grad_p,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_grad_p,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -1976,15 +2000,15 @@ void DEM::startTime()
darcy.rho_f,
dev_force,
dev_darcy_f_p);*/
- findDarcyPressureForceLinear<<<dimGrid, dimBlock>>>(
- dev_x,
+ KERNEL_LAUNCH(findDarcyPressureForceLinear, dimGrid, dimBlock, 0,
+ (dev_x,
dev_darcy_grad_p,
dev_darcy_phi,
wall0_iz,
darcy.rho_f,
darcy.bc_top,
dev_force,
- dev_darcy_f_p);
+ dev_darcy_f_p));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2008,8 +2032,8 @@ void DEM::startTime()
darcy.c_phi,
dev_darcy_phi,
dev_darcy_dphi);*/
- findDarcyPorositiesLinear<<<dimGridFluid, dimBlockFluid>>>(
- dev_cellStart,
+ KERNEL_LAUNCH(findDarcyPorositiesLinear, dimGridFluid, dimBlockFluid, 0,
+ (dev_cellStart,
dev_cellEnd,
dev_x_sorted,
dev_vel_sorted,
@@ -2020,7 +2044,7 @@ void DEM::startTime()
dev_darcy_phi,
dev_darcy_dphi,
dev_darcy_div_v_p,
- dev_darcy_vp_avg);
+ dev_darcy_vp_avg));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2029,23 +2053,21 @@ void DEM::startTime()
// copy porosities to the frictionless Y boundaries
if (grid.periodic == 2) {
- copyDarcyPorositiesToEdges<<<dimGridFluid,
- dimBlockFluid>>>(
- dev_darcy_phi,
+ KERNEL_LAUNCH(copyDarcyPorositiesToEdges, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_phi,
dev_darcy_dphi,
dev_darcy_div_v_p,
- dev_darcy_vp_avg);
+ dev_darcy_vp_avg));
cudaDeviceSynchronize();
}
// copy porosities to the frictionless lower Z boundary
if (grid.periodic == 2) {
- copyDarcyPorositiesToBottom<<<dimGridFluid,
- dimBlockFluid>>>(
- dev_darcy_phi,
+ KERNEL_LAUNCH(copyDarcyPorositiesToBottom, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_phi,
dev_darcy_dphi,
dev_darcy_div_v_p,
- dev_darcy_vp_avg);
+ dev_darcy_vp_avg));
cudaDeviceSynchronize();
}
@@ -2059,19 +2081,18 @@ void DEM::startTime()
+ darcy.p_mod_phi);
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyTopPressure<<<dimGridFluid, dimBlockFluid>>>(
- new_pressure,
+ KERNEL_LAUNCH(setDarcyTopPressure, dimGridFluid, dimBlockFluid, 0,
+ (new_pressure,
dev_darcy_p,
- wall0_iz);
+ wall0_iz));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setUpperPressureNS", iter);
// Modulate the pressures at the top wall
- setDarcyTopWallPressure
- <<<dimGridFluid, dimBlockFluid>>>(
- new_pressure,
+ KERNEL_LAUNCH(setDarcyTopWallPressure, dimGridFluid, dimBlockFluid, 0,
+ (new_pressure,
wall0_iz,
- dev_darcy_p);
+ dev_darcy_p));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setDarcyTopWallPressure",
iter);
@@ -2083,8 +2104,8 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- findDarcyPermeabilities<<<dimGridFluid, dimBlockFluid>>>(
- darcy.k_c, dev_darcy_phi, dev_darcy_k);
+ KERNEL_LAUNCH(findDarcyPermeabilities, dimGridFluid, dimBlockFluid, 0,
+ (darcy.k_c, dev_darcy_phi, dev_darcy_k));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2094,11 +2115,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_phi,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_phi,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2108,11 +2129,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodes<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_k,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_k,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2122,9 +2143,8 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- findDarcyPermeabilityGradients
- <<<dimGridFluid, dimBlockFluid>>>
- (dev_darcy_k, dev_darcy_grad_k);
+ KERNEL_LAUNCH(findDarcyPermeabilityGradients, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_k, dev_darcy_grad_k));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2133,16 +2153,16 @@ void DEM::startTime()
"Post findDarcyPermeabilityGradients", iter);
if (iter == 0) {
- setDarcyNormZero<<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_norm);
+ KERNEL_LAUNCH(setDarcyNormZero, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_norm));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post setDarcyNormZero", iter);
if (PROFILING == 1)
startTimer(&kernel_tic);
- copyValues<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
- dev_darcy_p_old);
+ KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
+ dev_darcy_p_old));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2154,15 +2174,14 @@ void DEM::startTime()
if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodesFlux<Float>
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(setDarcyGhostNodesFlux<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
darcy.bc_bot,
darcy.bc_top,
darcy.bc_bot_flux,
darcy.bc_top_flux,
dev_darcy_k,
- darcy.mu);
+ darcy.mu));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2192,9 +2211,9 @@ void DEM::startTime()
if (nijac == 0) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- copyValues<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
- dev_darcy_p_old);
+ KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
+ dev_darcy_p_old));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc,
@@ -2205,12 +2224,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodes<Float>
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2221,8 +2239,8 @@ void DEM::startTime()
if (nijac == 0) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- firstDarcySolution<<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(firstDarcySolution, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
dev_darcy_k,
dev_darcy_phi,
dev_darcy_dphi,
@@ -2240,7 +2258,7 @@ void DEM::startTime()
darcy.ndem,
wall0_iz,
dev_darcy_p_constant,
- dev_darcy_dp_expl);
+ dev_darcy_dp_expl));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc,
@@ -2252,8 +2270,8 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- updateDarcySolution<<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p_old,
+ KERNEL_LAUNCH(updateDarcySolution, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p_old,
//dev_darcy_dpdt,
dev_darcy_dp_expl,
dev_darcy_p,
@@ -2275,7 +2293,7 @@ void DEM::startTime()
wall0_iz,
dev_darcy_p_constant,
dev_darcy_p_new,
- dev_darcy_norm);
+ dev_darcy_norm));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2286,9 +2304,8 @@ void DEM::startTime()
if (darcy.bc_top == 1) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyTopWallFixedFlow
- <<<dimGridFluid, dimBlockFluid>>>
- (wall0_iz, dev_darcy_p);
+ KERNEL_LAUNCH(setDarcyTopWallFixedFlow, dimGridFluid, dimBlockFluid, 0,
+ (wall0_iz, dev_darcy_p));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc,
@@ -2301,15 +2318,14 @@ void DEM::startTime()
if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodesFlux<Float>
- <<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(setDarcyGhostNodesFlux<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
darcy.bc_bot,
darcy.bc_top,
darcy.bc_bot_flux,
darcy.bc_top_flux,
dev_darcy_k,
- darcy.mu);
+ darcy.mu));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc,
@@ -2322,9 +2338,9 @@ void DEM::startTime()
// Copy new values to current values
if (PROFILING == 1)
startTimer(&kernel_tic);
- copyValues<Float><<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p_new,
- dev_darcy_p);
+ KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p_new,
+ dev_darcy_p));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2390,8 +2406,8 @@ void DEM::startTime()
// Zero all dphi values right after they are used in fluid
// solution, unless a file is written in this step.
if (filetimeclock + time.dt < time.file_dt) {
- setDarcyZeros<Float> <<<dimGridFluid, dimBlockFluid>>>
- (dev_darcy_dphi);
+ KERNEL_LAUNCH(setDarcyZeros<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_dphi));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"After setDarcyZeros(dev_darcy_dphi)", iter);
@@ -2399,11 +2415,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- setDarcyGhostNodes<Float> <<<dimGridFluid, dimBlockFluid>>>
- (dev_darcy_p,
+ KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
darcy.bc_xn, darcy.bc_xp,
darcy.bc_yn, darcy.bc_yp,
- darcy.bc_bot, darcy.bc_top);
+ darcy.bc_bot, darcy.bc_top));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2413,12 +2429,12 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- findDarcyVelocities<<<dimGridFluid, dimBlockFluid>>>(
- dev_darcy_p,
+ KERNEL_LAUNCH(findDarcyVelocities, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_p,
dev_darcy_phi,
dev_darcy_k,
darcy.mu,
- dev_darcy_v);
+ dev_darcy_v));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2437,11 +2453,11 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- findShearStressOnFixedMovingParticles<<<dimGrid, dimBlock>>>
- (dev_x,
+ KERNEL_LAUNCH(findShearStressOnFixedMovingParticles, dimGrid, dimBlock, 0,
+ (dev_x,
dev_vel,
dev_force,
- dev_walls_tau_eff_x_pp);
+ dev_walls_tau_eff_x_pp));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2451,8 +2467,9 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
- summation<<<dimGrid, dimBlock>>>(dev_walls_tau_eff_x_pp,
- dev_walls_tau_eff_x_partial);
+ KERNEL_LAUNCH(summation, dimGrid, dimBlock, 0,
+ (dev_walls_tau_eff_x_pp,
+ dev_walls_tau_eff_x_partial));
cudaDeviceSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
@@ -2473,7 +2490,8 @@ void DEM::startTime()
// Update particle kinematics
if (PROFILING == 1)
startTimer(&kernel_tic);
- integrate<<<dimGrid, dimBlock>>>(dev_x_sorted,
+ KERNEL_LAUNCH(integrate, dimGrid, dimBlock, 0,
+ (dev_x_sorted,
dev_vel_sorted,
dev_angvel_sorted,
dev_x,
@@ -2496,7 +2514,7 @@ void DEM::startTime()
walls.tau_x[0],
change_velocity_state,
velocity_factor,
- blocksPerGrid);
+ blocksPerGrid));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post integrate", iter);
if (PROFILING == 1)
@@ -2510,8 +2528,9 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
if (walls.nw > 0) {
- summation<<<dimGrid, dimBlock>>>(dev_walls_force_pp,
- dev_walls_force_partial);
+ KERNEL_LAUNCH(summation, dimGrid, dimBlock, 0,
+ (dev_walls_force_pp,
+ dev_walls_force_partial));
}
cudaDeviceSynchronize();
if (PROFILING == 1)
@@ -2523,15 +2542,15 @@ void DEM::startTime()
if (PROFILING == 1)
startTimer(&kernel_tic);
if (walls.nw > 0) {
- integrateWalls<<< 1, walls.nw>>>(
- dev_walls_nx,
+ KERNEL_LAUNCH(integrateWalls, 1, walls.nw, 0,
+ (dev_walls_nx,
dev_walls_mvfd,
dev_walls_wmode,
dev_walls_force_partial,
dev_walls_acc,
blocksPerGrid,
time.current,
- iter);
+ iter));
}
cudaDeviceSynchronize();
if (PROFILING == 1)
@@ -2584,11 +2603,11 @@ void DEM::startTime()
// v_x, v_y, v_z -> v
if (fluid == 1 && cfd_solver == 0) {
- interpolateFaceToCenter<<<dimGridFluid, dimBlockFluid>>>(
- dev_ns_v_x,
+ KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
+ (dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z,
- dev_ns_v);
+ dev_ns_v));
cudaDeviceSynchronize();
checkForCudaErrorsIter("Post interpolateFaceToCenter", iter);
}
@@ -2601,8 +2620,8 @@ void DEM::startTime()
// Empty the dphi values after device to host transfer
if (fluid == 1) {
if (cfd_solver == 1) {
- setDarcyZeros<Float> <<<dimGridFluid, dimBlockFluid>>>
- (dev_darcy_dphi);
+ KERNEL_LAUNCH(setDarcyZeros<Float>, dimGridFluid, dimBlockFluid, 0,
+ (dev_darcy_dphi));
cudaDeviceSynchronize();
checkForCudaErrorsIter(
"After setDarcyZeros(dev_darcy_dphi) after transfer",
DIR diff --git a/src/integration.cuh b/src/integration.cuh
@@ -310,6 +310,7 @@ __global__ void integrate(
// Reduce wall force contributions from particles to a single value per wall
+#ifdef SPHERE_GPU
__global__ void summation(
const Float* __restrict__ in,
Float* __restrict__ out)
@@ -343,6 +344,25 @@ __global__ void summation(
if (cacheIdx == 0)
out[blockIdx.x] = cache[0];
}
+#else
+// CPU variant: thread 0 of each emulated block sums its chunk, preserving
+// the one-partial-per-block contract that integrate and integrateWalls rely
+// on. All launch sites use a grid covering devC_np, so the GPU grid-stride
+// loop reduces to the same per-block chunks.
+__global__ void summation(
+ const Float* __restrict__ in,
+ Float* __restrict__ out)
+{
+ if (threadIdx.x != 0)
+ return;
+ const unsigned int b = blockIdx.x;
+ Float sum = 0.0;
+ for (unsigned int i = b*blockDim.x;
+ i < (b+1)*blockDim.x && i < devC_np; ++i)
+ sum += in[i];
+ out[b] = sum;
+}
+#endif // SPHERE_GPU
// Update wall positions
__global__ void integrateWalls(
DIR diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
@@ -2,9 +2,11 @@
// CUDA implementation of porous flow
#include <iostream>
+#ifdef SPHERE_GPU
#include <cuda.h>
//#include <cutil_math.h>
#include <helper_math.h>
+#endif
#include "vector_arithmetic.h" // for arbitrary precision vectors
#include "sphere.h"
@@ -2880,7 +2882,7 @@ __global__ void findNormalizedResiduals(
// Computes the new velocity and pressure using the corrector
__global__ void updateNSpressure(
const Float* __restrict__ dev_ns_epsilon, // in
- const Float __restrict__ beta, // in
+ const Float beta, // in
Float* __restrict__ dev_ns_p) // out
{
// 3D thread index
DIR diff --git a/src/raytracer.cuh b/src/raytracer.cuh
@@ -404,7 +404,8 @@ __host__ void DEM::render(
const unsigned int blocksPerGrid = iDivUp(pixels, threadsPerBlock);
// Initialize image to background color
- imageInit<<< blocksPerGrid, threadsPerBlock >>>(dev_img, pixels);
+ KERNEL_LAUNCH(imageInit, blocksPerGrid, threadsPerBlock, 0,
+ (dev_img, pixels));
// Initialize camera values and transfer to constant memory
float imgw = grid.L[0]*1.35f; // Screen width in world coordinates
@@ -418,9 +419,9 @@ __host__ void DEM::render(
cameraInit(eye, lookat, imgw, focalLength);
// Construct rays for perspective projection
- rayInitPerspective<<< blocksPerGrid, threadsPerBlock >>>(
- dev_ray_origo, dev_ray_direction,
- make_float4(eye.x, eye.y, eye.z, 0.0f), width, height);
+ KERNEL_LAUNCH(rayInitPerspective, blocksPerGrid, threadsPerBlock, 0,
+ (dev_ray_origo, dev_ray_direction,
+ make_float4(eye.x, eye.y, eye.z, 0.0f), width, height));
cudaThreadSynchronize();
Float* linarr; // Linear array to use for color visualization
@@ -433,9 +434,9 @@ __host__ void DEM::render(
// Visualize spheres without color scale overlay
if (method == 0) {
- rayIntersectSpheres<<< blocksPerGrid, threadsPerBlock >>>(
- dev_ray_origo, dev_ray_direction,
- dev_x, dev_img);
+ KERNEL_LAUNCH(rayIntersectSpheres, blocksPerGrid, threadsPerBlock, 0,
+ (dev_ray_origo, dev_ray_direction,
+ dev_x, dev_img));
} else {
if (method == 1) { // Visualize pressure
@@ -511,11 +512,11 @@ __host__ void DEM::render(
}
// Start raytracing kernel
- rayIntersectSpheresColormap<<< blocksPerGrid, threadsPerBlock >>>(
- dev_ray_origo, dev_ray_direction,
+ KERNEL_LAUNCH(rayIntersectSpheresColormap, blocksPerGrid, threadsPerBlock, 0,
+ (dev_ray_origo, dev_ray_direction,
dev_x, dev_vel,
dev_linarr, maxval, lower_cutoff,
- dev_img);
+ dev_img));
}
DIR diff --git a/src/sorting.cuh b/src/sorting.cuh
@@ -53,7 +53,8 @@ __global__ void calcParticleCellID(unsigned int* dev_gridParticleCellID,
// Reorder particle data into sorted order, and find the start and end particle
// indexes of each cell in the sorted hash array.
-__global__ void reorderArrays(unsigned int* dev_cellStart,
+#ifdef SPHERE_GPU
+__global__ void reorderArrays(unsigned int* dev_cellStart,
unsigned int* dev_cellEnd,
unsigned int* dev_gridParticleCellID,
unsigned int* dev_gridParticleIndex,
@@ -131,5 +132,42 @@ __global__ void reorderArrays(unsigned int* dev_cellStart,
}
} // End of reorderArrays(...)
+#else
+// CPU variant: identical result, but reads the neighbor hash from global
+// memory instead of the shared-memory staging used on the GPU.
+__global__ void reorderArrays(unsigned int* dev_cellStart,
+ unsigned int* dev_cellEnd,
+ unsigned int* dev_gridParticleCellID,
+ unsigned int* dev_gridParticleIndex,
+ Float4* dev_x,
+ Float4* dev_vel,
+ Float4* dev_angvel,
+ Float4* dev_x_sorted,
+ Float4* dev_vel_sorted,
+ Float4* dev_angvel_sorted)
+{
+ unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;
+
+ if (idx < devC_np) { // Condition prevents block size error
+ const unsigned int cellID = dev_gridParticleCellID[idx];
+
+ // Find lowest and highest particle index in each cell
+ if (idx == 0 || cellID != dev_gridParticleCellID[idx-1]) {
+ dev_cellStart[cellID] = idx;
+ if (idx > 0)
+ dev_cellEnd[dev_gridParticleCellID[idx-1]] = idx;
+ }
+ if (idx == devC_np - 1)
+ dev_cellEnd[cellID] = idx + 1;
+
+ // Use the sorted index to reorder the position and velocity data
+ const unsigned int sortedIndex = dev_gridParticleIndex[idx];
+ dev_x_sorted[idx] = dev_x[sortedIndex];
+ dev_vel_sorted[idx] = dev_vel[sortedIndex];
+ dev_angvel_sorted[idx] = dev_angvel[sortedIndex];
+ }
+} // End of reorderArrays(...)
+#endif // SPHERE_GPU
+
#endif
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
DIR diff --git a/src/typedefs.h b/src/typedefs.h
@@ -1,7 +1,11 @@
#ifndef TYPEDEFS_H_
#define TYPEDEFS_H_
+#ifdef SPHERE_GPU
#include "vector_functions.h"
+#else
+#include "cpu_backend.h"
+#endif
//////////////////////
// TYPE DEFINITIONS //
DIR diff --git a/src/vector_arithmetic.h b/src/vector_arithmetic.h
@@ -1,7 +1,9 @@
#ifndef VECTOR_ARITHMETIC_H_
#define VECTOR_ARITHMETIC_H_
+#ifdef SPHERE_GPU
#include "cuda_runtime.h"
+#endif
#include "datatypes.h"
typedef unsigned int uint;
DIR diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -1,4 +1,9 @@
-find_package(PythonInterp REQUIRED)
+if (CMAKE_VERSION VERSION_LESS 3.12)
+ find_package(PythonInterp REQUIRED)
+else()
+ find_package(Python3 COMPONENTS Interpreter REQUIRED)
+ set(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
+endif()
add_test(io_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/io_tests.py)
@@ -24,8 +29,10 @@ add_test(io_tests_fluid ${PYTHON_EXECUTABLE}
add_test(porosity_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/porosity_tests.py)
-add_test(memory_tests ${PYTHON_EXECUTABLE}
- ${CMAKE_CURRENT_BINARY_DIR}/memcheck_tests.py)
+if (SPHERE_GPU)
+ add_test(memory_tests ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/memcheck_tests.py)
+endif()
add_test(bond_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/bond_tests.py)