feat(cpu): add CUDA execution-model emulation layer for OpenMP backend - 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 e9a992a493dfb106030c3d6e4c65b985cd562f20
DIR parent 2d9fa8ec01c60f97db3c4d4431cb0fb5f6e518c8
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Sat, 4 Jul 2026 23:41:03 +0200
feat(cpu): add CUDA execution-model emulation layer for OpenMP backend
Diffstat:
A src/cpu_backend.cpp | 45 +++++++++++++++++++++++++++++++
A src/cpu_backend.h | 254 +++++++++++++++++++++++++++++++
A src/device_cpu.cpp | 3 +++
A src/launch.h | 18 ++++++++++++++++++
4 files changed, 320 insertions(+), 0 deletions(-)
---
DIR diff --git a/src/cpu_backend.cpp b/src/cpu_backend.cpp
@@ -0,0 +1,45 @@
+// cpu_backend.cpp -- CPU-backend replacement for utility.cu plus the
+// definitions of the CUDA execution-model emulation globals.
+#include <iostream>
+#include <omp.h>
+
+#include "cpu_backend.h"
+#include "sphere.h"
+
+thread_local uint3 threadIdx;
+thread_local uint3 blockIdx;
+dim3 blockDim;
+dim3 gridDim;
+
+void DEM::diagnostics()
+{
+ // Retrieve information from device to host and run diagnostic tests
+ transferFromGlobalDeviceMemory();
+ checkValues();
+
+ // Clean up memory before exiting
+ if (fluid == 1 && cfd_solver == 0) {
+ freeNSmemDev();
+ freeNSmem();
+ }
+ if (fluid == 1 && cfd_solver == 1) {
+ freeDarcyMemDev();
+ freeDarcyMem();
+ }
+ freeGlobalDeviceMemory();
+ // CPU memory freed upon object destruction
+}
+
+// The CPU shims cannot fail asynchronously, so the checks are no-ops
+void DEM::checkForCudaErrors(const char*, const int) {}
+
+void DEM::checkForCudaErrorsIter(const char*, const unsigned int, const int) {}
+
+void DEM::initializeGPU()
+{
+ ndevices = 1;
+ if (verbose == 1)
+ std::cout << " CPU backend (OpenMP), max threads: "
+ << omp_get_max_threads() << std::endl;
+}
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
DIR diff --git a/src/cpu_backend.h b/src/cpu_backend.h
@@ -0,0 +1,254 @@
+// cpu_backend.h -- CUDA execution-model emulation for the OpenMP CPU backend.
+// Supplies the vector types, qualifier macros, thread-index globals and CUDA
+// runtime API shims needed to compile the kernel sources as plain C++.
+#ifndef CPU_BACKEND_H_
+#define CPU_BACKEND_H_
+
+#ifdef SPHERE_GPU
+#error "cpu_backend.h must not be included in CUDA builds"
+#endif
+
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <chrono>
+
+//// Vector types normally provided by the CUDA toolkit headers
+
+struct double2 { double x, y; };
+struct double3 { double x, y, z; };
+struct double4 { double x, y, z, w; };
+struct float3 { float x, y, z; };
+struct float4 { float x, y, z, w; };
+struct int3 { int x, y, z; };
+struct uint2 { unsigned int x, y; };
+struct uint3 { unsigned int x, y, z; };
+
+struct dim3 {
+ unsigned int x, y, z;
+ dim3(unsigned int x_ = 1, unsigned int y_ = 1, unsigned int z_ = 1)
+ : x(x_), y(y_), z(z_) {}
+};
+
+inline double2 make_double2(double x, double y)
+{
+ double2 v; v.x = x; v.y = y; return v;
+}
+
+inline double3 make_double3(double x, double y, double z)
+{
+ double3 v; v.x = x; v.y = y; v.z = z; return v;
+}
+
+inline double4 make_double4(double x, double y, double z, double w)
+{
+ double4 v; v.x = x; v.y = y; v.z = z; v.w = w; return v;
+}
+
+inline float3 make_float3(float x, float y, float z)
+{
+ float3 v; v.x = x; v.y = y; v.z = z; return v;
+}
+
+inline float4 make_float4(float x, float y, float z, float w)
+{
+ float4 v; v.x = x; v.y = y; v.z = z; v.w = w; return v;
+}
+
+inline int3 make_int3(int x, int y, int z)
+{
+ int3 v; v.x = x; v.y = y; v.z = z; return v;
+}
+
+inline uint2 make_uint2(unsigned int x, unsigned int y)
+{
+ uint2 v; v.x = x; v.y = y; return v;
+}
+
+inline uint3 make_uint3(unsigned int x, unsigned int y, unsigned int z)
+{
+ uint3 v; v.x = x; v.y = y; v.z = z; return v;
+}
+
+// Subset of helper_math.h operators used by the kernels
+inline int3 operator+(int3 a, int3 b)
+{
+ return make_int3(a.x + b.x, a.y + b.y, a.z + b.z);
+}
+
+// float3/float4 operators used by the raytracer
+inline float4 make_float4(float3 v, float w)
+{
+ return make_float4(v.x, v.y, v.z, w);
+}
+
+inline float3 operator-(const float3& a)
+{
+ return make_float3(-a.x, -a.y, -a.z);
+}
+
+inline float3 operator+(float3 a, float3 b)
+{
+ return make_float3(a.x + b.x, a.y + b.y, a.z + b.z);
+}
+
+inline float3 operator-(float3 a, float3 b)
+{
+ return make_float3(a.x - b.x, a.y - b.y, a.z - b.z);
+}
+
+inline float3 operator*(float3 a, float3 b)
+{
+ return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
+}
+
+inline float3 operator*(float3 a, float b)
+{
+ return make_float3(a.x*b, a.y*b, a.z*b);
+}
+
+inline float3 operator*(float b, float3 a)
+{
+ return make_float3(b*a.x, b*a.y, b*a.z);
+}
+
+inline float3 operator/(float3 a, float b)
+{
+ return make_float3(a.x/b, a.y/b, a.z/b);
+}
+
+inline float dot(float3 a, float3 b)
+{
+ return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+inline float3 cross(float3 a, float3 b)
+{
+ return make_float3(a.y*b.z - a.z*b.y,
+ a.z*b.x - a.x*b.z,
+ a.x*b.y - a.y*b.x);
+}
+
+inline float length(float3 v)
+{
+ return sqrtf(dot(v, v));
+}
+
+inline float3 normalize(float3 v)
+{
+ return v * (1.0f/sqrtf(dot(v, v)));
+}
+
+//// Qualifier macros and thread-index emulation
+
+#define __global__
+#define __device__
+#define __host__
+#define __constant__
+#define __forceinline__ inline
+#define __inline__ inline
+#define __shared__ static
+
+extern thread_local uint3 threadIdx;
+extern thread_local uint3 blockIdx;
+extern dim3 blockDim; // uniform per launch
+extern dim3 gridDim;
+
+// Safe as a no-op: the only kernels relying on intra-block cooperation
+// (reorderArrays, summation) have dedicated CPU implementations.
+inline void __syncthreads() {}
+
+//// CUDA runtime API shims
+
+enum cudaMemcpyKind {
+ cudaMemcpyHostToHost,
+ cudaMemcpyHostToDevice,
+ cudaMemcpyDeviceToHost,
+ cudaMemcpyDeviceToDevice
+};
+
+typedef int cudaError_t;
+const cudaError_t cudaSuccess = 0;
+
+// calloc, not malloc: some kernels (e.g. the Darcy solvers reading
+// dev_darcy_div_v_p) consume arrays no kernel ever writes, and rely on
+// fresh CUDA allocations being zeroed pages.
+inline cudaError_t cudaMalloc(void** p, size_t n)
+{
+ *p = calloc(1, n);
+ return cudaSuccess;
+}
+
+inline cudaError_t cudaFree(void* p)
+{
+ free(p);
+ return cudaSuccess;
+}
+
+inline cudaError_t cudaMemcpy(void* dst, const void* src, size_t n,
+ cudaMemcpyKind)
+{
+ memcpy(dst, src, n);
+ return cudaSuccess;
+}
+
+inline cudaError_t cudaMemset(void* p, int v, size_t n)
+{
+ memset(p, v, n);
+ return cudaSuccess;
+}
+
+inline cudaError_t cudaDeviceSynchronize() { return cudaSuccess; }
+inline cudaError_t cudaThreadSynchronize() { return cudaSuccess; }
+inline cudaError_t cudaGetLastError() { return cudaSuccess; }
+inline cudaError_t cudaSetDevice(int) { return cudaSuccess; }
+inline cudaError_t cudaDeviceReset() { return cudaSuccess; }
+
+#define cudaMemcpyToSymbol(sym, src, ...) memcpy(&(sym), (src), sizeof(sym))
+#define cudaMemcpyFromSymbol(dst, sym, ...) memcpy((dst), &(sym), sizeof(sym))
+
+// Event shims for the PROFILING timers
+typedef std::chrono::high_resolution_clock::time_point cudaEvent_t;
+
+inline void cudaEventCreate(cudaEvent_t*) {}
+
+inline void cudaEventRecord(cudaEvent_t& e, int = 0)
+{
+ e = std::chrono::high_resolution_clock::now();
+}
+
+inline void cudaEventSynchronize(cudaEvent_t&) {}
+
+inline void cudaEventElapsedTime(float* ms, cudaEvent_t& tic, cudaEvent_t& toc)
+{
+ *ms = std::chrono::duration<float, std::milli>(toc - tic).count();
+}
+
+inline void cudaEventDestroy(cudaEvent_t&) {}
+
+//// CPU launch driver: one OpenMP-parallel loop over emulated blocks*threads
+
+template<typename F>
+void cpuKernelLaunch(const dim3 grid, const dim3 block, size_t /*smem*/, F f)
+{
+ gridDim = grid;
+ blockDim = block;
+ const long nb = (long)grid.x * grid.y * grid.z;
+ const long nt = (long)block.x * block.y * block.z;
+ // The if clause keeps tiny launches serial and cheap
+#pragma omp parallel for schedule(static) collapse(2) if (nb*nt > 512)
+ for (long b = 0; b < nb; ++b) {
+ for (long t = 0; t < nt; ++t) {
+ blockIdx.x = b % grid.x;
+ blockIdx.y = (b / grid.x) % grid.y;
+ blockIdx.z = b / ((long)grid.x * grid.y);
+ threadIdx.x = t % block.x;
+ threadIdx.y = (t / block.x) % block.y;
+ threadIdx.z = t / ((long)block.x * block.y);
+ f();
+ }
+ }
+}
+
+#endif
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
DIR diff --git a/src/device_cpu.cpp b/src/device_cpu.cpp
@@ -0,0 +1,3 @@
+// CPU-backend translation unit: compiles device.cu as plain C++.
+#include "device.cu"
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
DIR diff --git a/src/launch.h b/src/launch.h
@@ -0,0 +1,18 @@
+// launch.h -- backend-agnostic kernel launch macro.
+// The argument list is passed parenthesized, e.g.
+// KERNEL_LAUNCH(interact, dimGrid, dimBlock, 0, (dev_x, dev_vel));
+// so commas inside it never split macro arguments, and templated kernels
+// like setNSghostNodes<Float3> work as the kernel argument.
+#ifndef LAUNCH_H_
+#define LAUNCH_H_
+
+#ifdef SPHERE_GPU
+#define KERNEL_LAUNCH(kernel, grid, block, smem, args) \
+ kernel<<<(grid), (block), (smem)>>>args
+#else
+#define KERNEL_LAUNCH(kernel, grid, block, smem, args) \
+ cpuKernelLaunch((grid), (block), (smem), [&] { kernel args; })
+#endif
+
+#endif
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4