9. Computing and Transferring Data on a GPU (CUDA)
We learn the ideas for running the for_each and reduce that we wrote in Chapter 7 on an NVIDIA GPU. We distinguish where to compute from where to place arrays, and proceed through CUDA kernels, allocating and transferring memory, the CUDA versions of MISO’s for_each and reduce, and the combination with MPI.
You can read this chapter even without a GPU. The examples to run list the CPU output and the GPU output side by side. If you proceed with the CPU only, after reading this chapter go on to the serial and MPI parts of Chapter 10, and add the GPU version once you have an environment where you can use a GPU.
Contents
1. Host and Device
Even in a program that uses a GPU, it is the CPU that drives the whole program. The CPU reads the configuration file, prepares the arrays, and runs the time-step loop. In the course of this, it asks the GPU to do the computation for each grid point.
In CUDA, the CPU side is called the host and the GPU side is called the device. The two have separate memories.
Host (CPU) |
Device (GPU) |
|
|---|---|---|
Memory |
Main memory |
Memory mounted on the GPU |
Concurrent execution |
As many as the number of cores |
Thousands to tens of thousands of threads |
MISO backend |
|
|
What matters here is that an array in host memory cannot be read from the device as it is, and that an array in device memory cannot be read from the host as it is. In Fortran, once an array was declared, any loop could read and write it. On a GPU, the programmer manages which memory an array is placed in.
host (CPU) device (GPU)
+--------------------+ +--------------------+
| main memory | cudaMemcpy | GPU memory |
| a_host[0..n-1] --+-------------->+-- a_dev[0..n-1] |
| <-+---------------+-- |
+--------------------+ (PCIe etc.) +--------------------+
loop control, I/O per-grid-point computation
Transfers between the host and the device go through a path that is much slower than reading and writing memory within the device. For this reason, MISO first places the arrays of physical quantities on the device and continues the computation inside the device during the time evolution. It copies to the host only when necessary, such as for output.
2. CUDA Kernels
2.1 Function specifiers
A function that runs on the GPU is called a kernel. CUDA is a language that extends C++, and a specifier placed before a function indicates where the function runs.
Specifier |
Runs on |
Called from |
Example |
|---|---|---|---|
|
Device |
Host (launched with |
A kernel. The return type is |
|
Device |
Device (a kernel or other |
A small function called from a kernel |
|
Host |
Host |
An ordinary C++ function (the same as having no specifier) |
|
Both |
Both |
A function used on both the host and the device |
A function marked __host__ __device__ is compiled twice, once for the host and once for the device. Many of MISO’s small functions have this form. For example, util::sqrt and util::min2 in miso/include/miso/utility.hpp are marked __host__ __device__, so they can be called from both a CPU loop and a GPU kernel.
CUDA code is compiled with the compiler nvcc, which is included in NVIDIA’s CUDA Toolkit. nvcc passes the host part to the usual C++ compiler (such as g++) and converts the device part into GPU instructions by itself.
2.2 Threads and blocks
A kernel is executed simultaneously by many threads. Every thread executes the same kernel instructions, but each uses its own number to change which element it handles.
Threads are grouped into sets called blocks, and an arrangement of blocks is called a grid. When launching a kernel, you specify the number of blocks in the grid and the number of threads per block with <<<grid, block>>>. Inside the kernel, the following built-in variables tell you your own position.
Variable |
Meaning |
|---|---|
|
The number of your own block (starting from 0) |
|
The number of threads per block |
|
Your own number within the block (starting from 0) |
For a one-dimensional array, the serial index is blockIdx.x * blockDim.x + threadIdx.x. Covering 1000 elements with 256 threads per block gives the following.
number of blocks = (1000 + 256 - 1) / 256 = 4 (integer division; rounds up)
blockIdx.x : 0 1 2 3
threadIdx.x: 0 ... 255 0 ... 255 0 ... 255 0 ... 255
global i : 0 ... 255 256 ... 511 512 ... 767 768 ... 1023
~~~~~~~~~~~~
i = 1000 to 1023 is outside the array.
do nothing and exit
Threads are launched in units of blocks, so the total number of threads does not exactly equal the number of elements. To keep the extra threads from writing outside the array, the kernel checks if (i < n).
The .x is there because blocks and grids can also be arranged in three dimensions (.x, .y, .z). We will see the real three-dimensional arrangement in MISO in Section 2.4.
2.3 A vector-addition kernel
This is a CUDA-only example, without MISO, that computes c[i] = a[i] + b[i] on the GPU. The functions that handle device memory (cudaMalloc, cudaMemcpy, cudaFree) are explained in the next section. Here, look at the form of defining and launching a kernel.
#include <cstdio>
#include <vector>
// function that runs on the GPU (kernel); one thread handles one element
__global__ void add_kernel(int n, const float *a, const float *b, float *c) {
int i = blockIdx.x * blockDim.x + threadIdx.x; // global index
if (i < n) { // threads past the end do nothing
c[i] = a[i] + b[i];
}
}
int main() {
const int n = 1000;
std::vector<float> a(n), b(n), c(n); // arrays on the host (CPU)
for (int i = 0; i < n; ++i) {
a[i] = 1.0f * i;
b[i] = 2.0f * i;
}
// allocate arrays on the device (GPU)
float *a_d = nullptr, *b_d = nullptr, *c_d = nullptr;
cudaMalloc(&a_d, sizeof(float) * n);
cudaMalloc(&b_d, sizeof(float) * n);
cudaMalloc(&c_d, sizeof(float) * n);
// host → device
cudaMemcpy(a_d, a.data(), sizeof(float) * n, cudaMemcpyHostToDevice);
cudaMemcpy(b_d, b.data(), sizeof(float) * n, cudaMemcpyHostToDevice);
// launch enough blocks of 256 threads to cover n elements
const int block = 256;
const int grid = (n + block - 1) / block; // round up: 4
add_kernel<<<grid, block>>>(n, a_d, b_d, c_d);
cudaError_t err = cudaGetLastError(); // check for a launch failure
if (err != cudaSuccess) {
std::printf("launch error: %s\n", cudaGetErrorString(err));
return 1;
}
// device → host (copy after waiting for the kernel to finish)
cudaMemcpy(c.data(), c_d, sizeof(float) * n, cudaMemcpyDeviceToHost);
std::printf("grid = %d, block = %d\n", grid, block);
std::printf("c[1] = %g, c[999] = %g\n", c[1], c[999]);
cudaFree(a_d);
cudaFree(b_d);
cudaFree(c_d);
return 0;
}
Name the file vecadd.cu, and compile and run it as follows. .cu is the extension for CUDA source. -arch=sm_89 specifies which generation of GPU to generate instructions for, and you match it to the GPU you use (Section 9).
nvcc -arch=sm_89 vecadd.cu -o vecadd
./vecadd
grid = 4, block = 256
c[1] = 3, c[999] = 2997
The body of the kernel corresponds to “one iteration” of the following Fortran loop.
do i = 1, n
c(i) = a(i) + b(i)
end do
There is nothing in the kernel that corresponds to the do of the loop. “Who computes which i” is determined by the number of each thread launched with <<<grid, block>>>. Therefore, the computation for each i must be independent of the others. A recurrence that uses the new value of c(i-1) to compute c(i) cannot be handed to the GPU as it is. This is the same condition as for the lambda passed to for_each in Chapter 7.
2.4 MISO’s hand-written kernels
Most of MISO’s MHD time evolution is written directly as __global__ kernels. This is the beginning of bb_ht_vb_kernel in miso/include/miso/mhd_integrator_cuda.cuh.
template <typename Real>
__global__ void bb_ht_vb_kernel(FieldsView<const Real> qq_argm,
Array3DView<Real> pr, Array3DView<Real> bb,
Array3DView<Real> ht, Array3DView<Real> vb,
GridView<const Real> grid) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
int k = blockIdx.z * blockDim.z + threadIdx.z;
if (i < 0 || i >= grid.i_total || j < 0 || j >= grid.j_total || k < 0 ||
k >= grid.k_total)
return;
The blocks and the grid are arranged in three dimensions, and i, j, and k are built from the numbers in each direction. return ends the threads that stick out beyond the grid. The arguments are not Array3D, which owns the array, but the lightweight Array3DView and FieldsView (Chapter 4). Kernel arguments are passed to the device by value (explained in Section 3.3), so we pass Views, which hold only a pointer and sizes.
The shape of a block is determined by KernelShape3D in miso/include/miso/cuda_util.cuh as 8×8×8 = 512 threads. The launching side is in Integrator<Real, backend::CUDA>::update_sc4.
bb_ht_vb_kernel<Real><<<cu_shape.grid_dim, cu_shape.block_dim>>>(
qq_argm.view(), pr.view(), bb.view(), ht.view(), vb.view(), cgrid);
MISO_CUDA_CHECK(cudaGetLastError());
MISO_CUDA_CHECK is a macro that checks the error returned by a CUDA function (Section 7).
Exercise 1: Number of blocks and the out-of-range check
We process n = 300 elements with 128 threads per block. Find the number of blocks and the total number of launched threads, and explain what happens without if (i < n).
The number of blocks is (300 + 128 - 1) / 128 = 427 / 128 = 3 (integer division), and the total number of threads is 3 × 128 = 384. The serial index i runs from 0 to 383, and the 84 threads from 300 to 383 point outside the array.
Without if (i < n), these 84 threads write c[300] to c[383]. That is outside the 300 elements allocated with cudaMalloc, so it either corrupts other arrays or stops the program as an illegal memory access. It is the same kind of mistake as writing beyond the upper bound of an array in Fortran, but on a GPU many threads write at the same time, so it becomes hard to trace where the result shows up.
3. Allocating and Transferring Device Memory
3.1 cudaMalloc and cudaMemcpy
Device memory is allocated, transferred, and freed with CUDA functions.
Function |
Role |
Closest in Fortran and C++ |
|---|---|---|
|
Allocate memory on the device and put its address into |
|
|
Copy the contents of memory |
Array assignment |
|
Free device memory |
|
The direction of cudaMemcpy is specified by the fourth argument.
Direction |
Copy direction |
|---|---|
|
Host to device |
|
Device to host |
|
Within the device |
What cudaMalloc returns is a device address. Its type is an ordinary pointer, float *, but you cannot read it in host code like a_d[0] (because it is not an address in host memory). Conversely, if you pass a host address, such as a.data() of a std::vector, to a kernel, the device cannot read it. The type of a pointer alone does not tell which memory it points to. This is why, in the example of Section 2.3, we added _d to the variables on the device side.
3.2 MISO’s Array3D<T, backend::CUDA>
In MISO, instead of writing cudaMalloc and the like directly, the second template argument of the array class selects where the memory is (Chapter 6).
Array3D<Real, backend::Host> a_host(nx, ny, nz); // host memory
Array3D<Real, backend::CUDA> a_dev(nx, ny, nz); // device memory
These are the constructor and destructor of the CUDA version in miso/include/miso/array3d.hpp.
template <typename T> class Array3D<T, backend::CUDA> {
private:
T *data_ = nullptr;
std::array<int, 3> shape_ = {-1, -1, -1};
public:
Array3D(int nx0, int nx1, int nx2) : shape_{nx0, nx1, nx2} {
// (size check omitted)
MISO_CUDA_CHECK(cudaMalloc(&data_, sizeof(T) * array_size));
}
~Array3D() {
if (data_)
cudaFree(data_);
data_ = nullptr;
}
Where the host version uses new T[...] and delete[], this version uses cudaMalloc and cudaFree. As with the RAII of Chapter 4, the memory is released automatically when the variable’s lifetime ends.
The member function copy_from takes care of transfers. The type of the argument (which backend’s array it is) determines the direction of cudaMemcpy.
Call |
Copy direction |
Direction used inside |
|---|---|---|
|
Host to device |
|
|
Device to host |
|
|
Within the device |
|
|
Within the host |
|
This is the overload selection that we learned in Chapter 6. The host version Array3D<T, backend::Host> also has a copy_from that “copies from the CUDA version” inside #ifdef __CUDACC__. __CUDACC__ is a macro defined only while nvcc is compiling.
All of Array1D to Array4D have a CUDA version of the same form.
3.3 Why the CUDA version of an array cannot be indexed
The host version of Array3D has an operator() so that you can write a(i, j, k), but the CUDA version does not. This is because even if you wrote a_dev(i, j, k) in host code, the host cannot read device memory.
Device elements are read and written through a View inside a kernel. The member functions of Array3DView are marked __host__ __device__.
/// @brief Return reference to the element at the given indices.
__host__ __device__ T &operator()(int i0, int i1, int i2) const noexcept {
// (index range assert omitted)
return data_[(i0 * shape_[1] + i1) * shape_[2] + i2];
}
A View is a small value that holds only a pointer and the size in each direction. When you pass it by value to a kernel argument or to the lambda of the next section, only these few values are copied to the device. The elements are already in device memory, and the View just points to their address. In other words, the roles are divided as follows: the CUDA version of Array3D is used in host code for allocating, freeing, copy_from, and view(), and the elements are read and written through a View inside kernels.
Exercise 2: Where an array lives and the direction of transfer
InitialCondition::apply in demo/mhd2d_vortex/src/main.cpp writes the initial values, like qq.ro(i, j, k) = 1.0;, with a plain for loop. In a CUDA build the arrays of physical quantities are on the device, so this loop cannot write the device arrays directly. Read MHD::apply_initial_condition (miso/include/miso/mhd_solver.hpp), which is called from ModelBase::run() (miso/include/miso/mhd_model_base.hpp), and explain which array the initial condition is written to and how it reaches the device.
apply_initial_condition creates the host grid grid_h (Grid<Real, backend::Host>) and the host array qq_h (Fields<Real, backend::Host>) as local variables. With ic.apply(qq_h.view(), grid_h.const_view()), the for loop of InitialCondition::apply writes the initial values to the host array. Then qq.copy_from(qq_h) copies them to the array qq for computation (Fields<Real, Backend>).
If Backend is backend::CUDA, this copy_from is the “host to device” case of the table in Section 3.2, and cudaMemcpyHostToDevice is used inside. If Backend is backend::Host, it is a “within the host” copy. In either build, the procedure is the same: write the initial condition on the host and copy it to the array for computation. qq_h and grid_h are local variables, so they are freed when the function exits.
Output goes in the opposite direction. Checkpoint::save in miso/include/miso/mhd_checkpoint.hpp copies the device values into the host array qq with qq.copy_from(qq_) and then writes them to a file.
4. The CUDA Version of for_each
4.1 Running the same lambda on a GPU
In Chapter 7, we passed a lambda to for_each to write the computation for each grid point. The Host version of for_each merely called the lambda f as f(i) inside a for loop. In the CUDA version, the same lambda is passed to the GPU as a kernel argument, and each thread calls f(i). This is the one-dimensional CUDA version in miso/include/miso/execution.hpp.
template <class F> __global__ void for_each_kernel(Range1D range, F f) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < range.size()) {
int i = range.unravel(idx);
f(i);
}
}
/// @brief 1D parallel (concurrent) for loop in the CUDA backend.
template <class F> inline void for_each(backend::CUDA, Range1D range, F f) {
int n = range.size();
if (n == 0)
return; // Avoid launching kernel with zero threads
int block = 256;
int grid = (n + block - 1) / block;
for_each_kernel<<<grid, block>>>(range, f);
MISO_CUDA_CHECK(cudaGetLastError());
}
It has the same form as the vector addition of Section 2.3. The difference is that the body of the computation, c[i] = a[i] + b[i], is not written directly in the kernel but passed in from outside as the argument f. range.unravel(idx) converts the serial index idx (starting from 0) into the index i counted from the start of the range, range.begin.
In Chapter 7, we learned that a lambda can be thought of as “a small struct that has the captured values as its members”. Kernel arguments are passed by value, so this whole struct is copied to the device.
host-side code device (GPU)
-------------------------------------- -----------------------------------------
auto fv = f.view(); // pointer + size
const Real dx = ...;
MISO_LAMBDA(int i) { fv[i] = ... dx ...; }
│ struct generated by the compiler:
│ struct LambdaType {
│ Array1DView<Real> fv; // copy
│ Real dx; // copy
│ __device__ void operator()(int i) const;
│ };
▼
for_each(backend::CUDA{}, Range1D{0, n}, lambda)
│
└─ for_each_kernel<<<grid, 256>>>(range, f)
│ args range and f (values of fv, dx)
│ copied to device at launch ─────────▶ all threads read the same copy of f
thread idx=0 : f(range.begin + 0)
thread idx=1 : f(range.begin + 1)
...
thread idx=n-1 : f(range.begin + n-1)
threads with idx >= n do nothing
What is copied to the device is only the values the lambda captured (the pointers and sizes of Views, and numbers such as dx). The elements of the arrays that the Views point to are in device memory from the start. For that reason, the Views that a lambda captures must be made from device arrays (such as Array1D<Real, backend::CUDA>). If you pass a View of a host array to the CUDA version of for_each, the device would read a host address, which causes a run-time error.
The three-dimensional for_each(backend::CUDA, Range3D, f) also launches the threads arranged in one dimension. Range3D::unravel converts the serial index into (i0, i1, i2), assigning them in the order in which the last index changes fastest. This is so that adjacent threads read elements that are adjacent in memory.
4.2 MISO_LAMBDA and __host__ __device__
A lambda called on the device needs the __device__ specifier, just like a function. MISO_LAMBDA is a macro that adds it according to the build (miso/include/miso/execution.hpp).
#ifdef __CUDACC__
#define MISO_LAMBDA [=] __host__ __device__
#else
#define MISO_LAMBDA [=]
#endif // __CUDACC__
When compiling with nvcc it becomes [=] __host__ __device__, and when compiling with g++ or the like it becomes [=]. [=] is the specification we learned in Chapter 7 for “capturing the outer variables it uses by value”. We do not use [&], which captures by reference, because a reference points to a host variable and cannot be read from the device.
A lambda marked __host__ __device__ is, in CUDA, a feature called an extended lambda, and it can be used only when the --expt-extended-lambda option is given to nvcc. MISO adds this option in miso/CMakeLists.txt.
target_compile_options(
MISO INTERFACE $<$<COMPILE_LANGUAGE:CUDA>:--expt-extended-lambda>)
$<$<COMPILE_LANGUAGE:CUDA>:...> is CMake notation meaning “add this only to files compiled as CUDA”. When you link MISO::MISO, this option is added to your own program as well.
When building with the CPU only, __host__ and __device__ are words that g++ does not know. miso/include/miso/cuda_compat.hpp defines them as empty words in that case.
#else // __CUDACC__(when compiling with a compiler other than nvcc)
#ifndef __device__
#define __device__
#endif
#ifndef __host__
#define __host__
#endif
// Do not define __global__:
// CUDA kernel semantics and launch syntax are not portable.
Because of this, functions and Views marked __host__ __device__ can be used as ordinary C++ in a CPU build. On the other hand, __global__ is not defined. Kernels and <<< >>> launches are meaningless in a CPU build, so the CUDA versions of for_each and reduce and the hand-written kernels are placed inside #ifdef __CUDACC__.
4.3 Running the same main.cpp on a CPU and a GPU
This is an example that uses for_each and reduce to compute \(\int_0^1 x^2\,dx = 1/3\) with the midpoint rule. Just by switching Backend, the same code runs on both a CPU and a GPU. How the CUDA version of reduce works is explained in Section 6. Here, note that the way of calling it is the same as for the Host version.
#include <cstdio>
#include <miso/core.hpp>
using namespace miso;
using Real = float;
#ifdef USE_CUDA
using Backend = backend::CUDA;
#else
using Backend = backend::Host;
#endif
int main(int argc, char **argv) {
Env env(argc, argv); // initialize MPI; in a CUDA build, also select the GPU
const int n = 1000;
const Real dx = Real(1) / n;
Array1D<Real, Backend> f(n); // allocated in GPU memory if Backend is CUDA
auto fv = f.view(); // lightweight View passed to the lambda
// compute f[i] = x^2 for each element (x is the cell center)
for_each(Backend{}, Range1D{0, n}, MISO_LAMBDA(int i) {
const Real x = (i + Real(0.5)) * dx;
fv[i] = x * x;
});
// sum of f[i] * dx (integral of x^2 from 0 to 1; exact value is 1/3)
const Real integral = reduce(
Backend{}, Range1D{0, n}, Real(0),
MISO_LAMBDA(int i) { return fv[i] * dx; },
MISO_LAMBDA(Real a, Real b) { return a + b; });
// for checking, copy to a CPU array and print one element
Array1D<Real, backend::Host> f_host(n);
f_host.copy_from(f);
if (mpi::is_root()) {
std::printf("f[999] = %.7f\n", f_host[999]);
std::printf("integral = %.7f\n", integral);
}
}
Here are the key points for reading it.
fvanddxare local variables prepared outside the lambda. The lambda captures these two by value.fitself (anArray1D) is not used inside the lambda. The owning array stays on the host, and only the View is passed.The reduction operation
op(the second lambda) is also called on the device, so we write it withMISO_LAMBDA.To display the result on the host, copy it to the host with
f_host.copy_from(f).USE_CUDAis a macro that your own CMake configuration defines for CUDA builds (we gave it the same name as in MISO’s demos).
These are the outputs of the CPU build and the GPU build (the GPU was an RTX 4090).
# CPU (backend::Host)
f[999] = 0.9990003
integral = 0.3333332
# GPU (backend::CUDA)
f[999] = 0.9990003
integral = 0.3333333
f[999] matches, but integral differs in the last digit. The rounding error of float addition changes with the order of addition. The Host version adds in the order i = 0, 1, 2, ..., but the CUDA version adds by combining in a tree shape, as we will see in Section 6. Decide how to compare results by assuming that the results of a CPU and a GPU may not match bit for bit.
4.4 CMake configuration
To build this program with CUDA, add the following to the CMake configuration we wrote in Chapter 1. The CMakeLists.txt of MISO’s demos (demo/mhd2d_vortex/CMakeLists.txt) does the same thing, so we show an excerpt.
option(USE_CUDA "Enable CUDA support" OFF)
if(USE_CUDA)
# Select GPU architecture before enable_language(CUDA); otherwise CMake
# defines the default (52) and this setting is ignored.
if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
set(CMAKE_CUDA_ARCHITECTURES
native
CACHE STRING "CUDA architectures")
endif()
enable_language(CUDA)
set(CMAKE_CUDA_STANDARD 17)
set(CMAKE_CUDA_STANDARD_REQUIRED ON)
endif(USE_CUDA)
if(USE_CUDA)
find_package(CUDAToolkit REQUIRED)
target_link_libraries(${app_name} PRIVATE CUDA::cudart)
# Force main.cpp to be treated as CUDA source
set_source_files_properties(${app_main} PROPERTIES LANGUAGE CUDA)
target_compile_definitions(${app_name} PRIVATE USE_CUDA)
endif(USE_CUDA)
enable_language(CUDA) makes CMake look for nvcc. CMAKE_CUDA_ARCHITECTURES (Section 9) is decided before that. set_source_files_properties(... LANGUAGE CUDA) makes a file with the extension .cpp be compiled by nvcc instead of g++. This is so that the same file is used for the CPU build and the GPU build. target_compile_definitions(... USE_CUDA) enables the #ifdef USE_CUDA in the source. In addition, before loading MISO, the demo turns ON MISO’s own option MISO_USE_CUDA, which makes MISO attach the CUDA run-time library and --expt-extended-lambda. When you load MISO with add_subdirectory, you likewise write set(MISO_USE_CUDA ON CACHE BOOL "...") before it. Use separate build folders for the CPU and the GPU.
# CPU
cmake -B build_cpu -S .
cmake --build build_cpu -j
# GPU (nvcc must be on PATH; 89 is an example for the RTX 40 series)
cmake -B build_gpu -S . -DUSE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=89
cmake --build build_gpu -j
Exercise 3: What the lambda captures
In the example of Section 4.3, list the variables captured by the lambda passed to for_each, and estimate roughly how many bytes each takes when copied to the device. Array1DView<Real> has one pointer and one int (miso/include/miso/array1d.hpp). Also explain what goes wrong if you write f.view()[i] instead of fv[i] inside the lambda.
It captures two things: fv (Array1DView<Real>) and dx (Real = float). fv is an 8-byte pointer and a 4-byte int, about 16 bytes including the struct padding, and dx is 4 bytes. In total this is about 20 to 24 bytes, and the 1000 elements of the array (4000 bytes) are not copied. n was used only to build Range1D{0, n} and is not used inside the lambda, so it is not captured.
If you write f.view() inside the lambda, it tries to capture f (the owning Array1D<Real, backend::CUDA>). Array1D is not copyable (the = delete of Chapter 4), so it cannot be captured by value with [=], which gives a compile error (nvcc displays ... cannot be referenced -- it is a deleted function). Even if it could be captured, view() is not marked __device__, so it cannot be called from the device. Keep to the division of roles: handle owning arrays in host code, and pass only Views to lambdas.
5. The Pitfall of Capturing this
5.1 Using a member variable captures this
When you write a lambda inside a member function of a class and use a member variable in it, a mistake can occur that works on a CPU but stops on a GPU.
As we learned in Chapter 4, writing the member variable factor inside a member function is shorthand for this->factor. this is a pointer to the object on which the member function was called. By the rules of C++, when you use a member variable inside a lambda, what is captured, even with [=], is not the value of factor but the pointer this (Chapter 7).
On a GPU, this becomes a problem. this is the address of an object in host memory. The value of this (an address) is copied to the device together with the lambda, but the object that address points to is not. When a device thread tries to read this->factor, it ends up reading a host address.
host memory device (GPU)
+------------------+ lambda copy
| Scaler s | this = 0x7fff.... ──┐
| factor = 3 | ◀── 0x7fff.... │ not readable
+------------------+ │ from the device
reads this->factor ◀──┘ → illegal memory access
The following example compares two member functions that do the same thing (code for explanation).
struct Scaler {
Real factor = 3;
// wrong: use the member variable factor directly inside the lambda
void apply_wrong(Array1DView<Real> a) const {
for_each(Backend{}, Range1D{0, a.size()},
MISO_LAMBDA(int i) { a[i] = factor * i; });
}
// correct: copy to a local variable before the lambda
void apply_right(Array1DView<Real> a) const {
const Real factor_ = factor;
for_each(Backend{}, Range1D{0, a.size()},
MISO_LAMBDA(int i) { a[i] = factor_ * i; });
}
};
In main, we called apply_right and apply_wrong in turn, copied to the host with a_host.copy_from(a) after each, and displayed a[3].
# CPU (backend::Host)
right: a[3] = 9
wrong: a[3] = 9
# GPU (backend::CUDA)
right: a[3] = 9
CUDA Error: an illegal memory access was encountered .../miso/include/miso/array1d.hpp 109
On a CPU there is no distinction between host and device, so it works correctly even if this is captured. Therefore, this mistake is not found even if the CPU tests pass. On a GPU, the kernel of apply_wrong stops with an illegal memory access. The reason the reported location of the error is array1d.hpp (the cudaMemcpy inside copy_from) is explained in Section 7.
In apply_right, factor is copied to the local variable factor_ just before the lambda. The lambda captures the value of the local variable factor_ (one float), so this is not captured. a is a function argument (the value of a View), so it can be used as it is.
5.2 How to write it in MISO
MISO’s dev/notes/usage_ai.md describes this rule with examples.
const Real gm_ = gm; // copy the member variable to a local
for_each(btag, range, MISO_LAMBDA(int i) { pr[i] = (gm_ - Real(1)) * qq.ro[i] * qq.ei[i]; });
The real code is IdealEOS::gas_pressure in miso/include/miso/eos.hpp.
void gas_pressure(Backend btag, mhd::FieldsView<const Real> qq,
Array3DView<Real> pr) const {
const Real gm_ = gm; // device lambdas must not capture `this`
Range1D range{0, qq.size()};
for_each(
btag, range,
MISO_LAMBDA(int i) { pr[i] = (gm_ - Real(1)) * qq.ro[i] * qq.ei[i]; });
}
gm (the ratio of specific heats) is a member variable of IdealEOS, so it is copied to gm_ before the lambda. qq and pr are function arguments (Views), so they can be captured as they are. In Integrator<Real, backend::CUDA>::cfl (miso/include/miso/mhd_integrator_cuda.cuh) as well, the member variable cfl_number is copied to cfl_number_, and the arrays are also turned into local Views qq_v, grid_v, and cs_v before the lambda is written.
This rule comes from a bug found in MISO’s issue #170 (https://github.com/is-isee/MISO/issues/170). The lambdas in eos.hpp, cfl, table_interpolator.hpp, and demo/mhd3d_magnetosphere used member variables and captured this, and on an RTX 4090 environment mhd2d_vortex and mhd3d_magnetosphere stopped with illegal memory access. On the other hand, they happened to work in another environment. In an environment where HMM (Heterogeneous Memory Management), a mechanism that lets the device read host memory, is enabled, reading a host address does not stop the program. After we fixed them to copy the members into local variables just before the lambda, we confirmed that the magnetosphere demo runs on 2 GPUs and agrees with the CPU within the precision of float.
“It worked on some GPU environment” is not a confirmation that this is not captured, either. After writing a lambda, go back to the declarations and check whether each name used inside is a local variable or an argument.
Exercise 4: Checking the names a lambda uses
Read IdealEOS::sound_speed in miso/include/miso/eos.hpp, and classify the names used inside the lambda into “local variables”, “function arguments”, and “member variables”. Also explain what happens in the CPU build and the GPU build if you delete the line const Real gm_ = gm; and change the inside of the lambda back to gm.
The lambda of sound_speed computes cs[i] = util::sqrt(gm_ * (gm_ - Real(1)) * qq.ei[i]);.
Name |
Kind |
|---|---|
|
Local variable (a copy of the member |
|
Function argument ( |
|
Function argument ( |
|
Lambda argument |
util::sqrt is a function name, not a variable, so it is not captured. It is marked __host__ __device__, so it can be called from the device. There are no member variables inside the lambda.
If you drop gm_ and write gm, gm means this->gm, and the lambda captures this. In the CPU build, the object that this points to can be read as it is, so it works correctly. In the GPU build, it compiles, but when run the device reads a host address and stops with illegal memory access (in an environment with HMM enabled, it may happen to work).
6. The CUDA Version of reduce
6.1 Producing one value with many threads
reduce was the operation that produces a single value, such as the minimum or the sum, from a whole array (Chapter 7). The Host version simply combines the elements one at a time into a variable, acc = init, with acc = op(acc, map(i)). Doing the same on a GPU takes some care. If thousands of threads rewrite the single variable acc at the same time, they overwrite each other’s writes and the value is not correct. So MISO’s CUDA version uses a method that combines values by halves inside a block. Like a sports tournament, it pairs the values two by two, keeps the winner (the sum for a total, the smaller one for a minimum), and repeats until one remains.
6.2 Shared memory and reduction within a block
The threads in a block can share a small, fast memory called shared memory. Shared memory exists separately for each block and cannot be seen from other blocks.
This shows how the sum of 8 values is found with 8 threads per block.
shared mem s[]: s0 s1 s2 s3 s4 s5 s6 s7 each thread writes its own value
│ │ │ │ │ │ │ │
│ │ │ │ │ │ │ │ __syncthreads() (wait until everyone has written)
offset = 4: s0+s4 s1+s5 s2+s6 s3+s7 threads 0-3 add the value 4 ahead
│ │ │ │ __syncthreads()
offset = 2: (s0+s4)+(s2+s6) (s1+s5)+(s3+s7) threads 0-1
│ │ __syncthreads()
offset = 1: s0+s1+...+s7 thread 0
│
block result → out[blockIdx.x]
The 8 values become 1 in 3 stages (\(\log_2 8\)). In each stage, half as many threads work as in the previous stage.
__syncthreads() is an instruction that waits until all threads in the block reach that point. Without it, when one thread reads s[tid + offset], the thread that is supposed to write that value may not have written it yet.
The real code is block_reduce_shared in miso/include/miso/execution.hpp.
template <class T, class Op>
__device__ inline T block_reduce_shared(T *s, int tid, int n, Op op) {
__syncthreads();
for (int offset = (n + 1) >> 1; offset > 0; offset = (offset + 1) >> 1) {
if (tid < offset && tid + offset < n) {
s[tid] = op(s[tid], s[tid + offset]);
}
__syncthreads();
if (offset == 1)
break;
}
return s[0];
}
>> 1 is the operation of dividing by 2 (shifting right by one bit), and (n + 1) >> 1 is \(n/2\) rounded up. We round up so that it works even when the number of threads n is not a power of 2. With the condition tid + offset < n, a thread that has no partner does nothing.
6.3 Combining the block results further
What a one-block reduction gives is the result for that block’s share. If there are 64 blocks, 64 values remain. MISO’s reduce repeats the kernel to reduce these until one remains.
stage 1: reduce_stage_kernel compute map(i) and reduce per block
n elements ──▶ buf0[grid] (grid = n / 256, rounded up)
stage 2 onward: reduce_stage_array_kernel
buf0[grid] ──▶ buf1[g] ──▶ buf0[...] ── ... ──▶ 1 value
(swap buf0 and buf1 and repeat until one value remains)
last: copy the single value to the host with cudaMemcpyAsync and wait with cudaStreamSynchronize
In reduce_stage_kernel of stage 1 (miso/include/miso/execution.hpp), the threads inside the range make their own values with map, and the threads outside the range write init, each into the shared memory s[tid]. Then block_reduce_shared is called, and thread 0 writes the block’s result to out[blockIdx.x]. To change how map is called for one-dimensional and three-dimensional ranges, it uses the if constexpr of Chapter 6.
buf0 and buf1, which hold intermediate results, are device memory that ReduceHelper<T> allocates with cudaMalloc and keeps. Allocating and freeing them on every call of reduce takes time, so in places where it is called many times, you prepare a ReduceHelper and reuse it. Integrator<Real, backend::CUDA> has ReduceHelper<Real> cfl_helper; as a member, for cfl, which finds the time step.
const auto op = MISO_LAMBDA(Real a, Real b) { return util::min2(a, b); };
const Real dt_max = 1.e10;
const auto dt = reduce(backend::CUDA{}, range, dt_max, f, op, cfl_helper);
There is also a form that omits the last argument cfl_helper, in which case a temporary ReduceHelper is created inside reduce.
6.4 Make the initial value the identity element of op
In the CUDA version, the threads outside the range write init to shared memory, and it also takes part in the reduction. Therefore init must be a value that does not change the partner’s value when combined with op (the identity element).
Reduction |
|
|
|---|---|---|
Sum |
|
|
Minimum |
|
A sufficiently large value ( |
Maximum |
|
A sufficiently small value |
In the Host version, init takes part just once as the “first value”, so it works even if it is not the identity element. We confirmed that the same call gives different results in the Host version and the CUDA version (an example that sums 10 elements that are all 1).
init = 0: CPU 10, GPU 10
init = 5: CPU 15, GPU 1240 ← the 5 for each of the 246 out-of-range threads is also added
6.5 Order of reduction and rounding error
The Host version combines from the front in order, and the CUDA version combines in a tournament shape. For a minimum or a maximum the value is the same regardless of order, but for a floating-point sum the rounding error changes with the order. This is why the last digit of integral differed between the CPU and the GPU in the example of Section 4.3. When you compare CPU and GPU results, decide on the tolerance by referring to Testing and Comparing Results.
The size of the shared memory used in the reduction within a block (number of threads × sizeof(T)) changes with each launch. Shared memory of this kind is declared inside the kernel with extern __shared__ without writing its size, and the size is passed as the third value in <<<grid, block, bytes>>> at launch. This is why MISO’s reduce launches with four values, <<<grid, block, helper.shmem, stream>>> (the fourth is the stream of Section 7).
Earlier versions of MISO declared it as follows inside a template kernel.
extern __shared__ T s[];
It looks as if a separate array is made for each T, but an extern __shared__ array is treated as a single entity with the same name regardless of the template type. Therefore, using the reduce for int and the reduce for double in the same file produced a compile error saying “the same s was declared with different types” (MISO issue #190, https://github.com/is-isee/MISO/issues/190). (The error is error: declaration is incompatible with previous "s".) The current MISO declares the shared memory just once as an untyped sequence of bytes, and the using side reinterprets it as T *.
template <class T> __device__ inline T *reduce_shared_memory() {
extern __shared__ __align__(16) unsigned char reduce_smem[];
static_assert(alignof(T) <= 16, "Unsupported alignment for reduce");
return reinterpret_cast<T *>(reduce_smem);
}
__align__(16) aligns the starting address to a multiple of 16. This is because, to place a double and the like, the address must be a multiple that fits the size of the type. static_assert stops the compilation if a type that requires stricter alignment is used (Chapter 6). miso/tests/serial/test_execution.cu has a test that uses the reduce for int and double in the same file.
How extern __shared__ is handled is explained in the section on shared memory in NVIDIA’s CUDA C++ Programming Guide (https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html).
Exercise 5: Counting the reduction stages
In the default configuration of mhd2d_vortex (demo/mhd2d_vortex/config.yaml), i_size: 128, j_size: 128, k_size: 1, and MPI divides into two in the x direction. Find the number of elements that the reduce of cfl handles on each rank (the number of interior grid points, excluding ghost cells), and, with 256 threads per block, find the number of blocks in stage 1 and the number of kernel launches in stage 2 and after.
It is divided into two in the x direction, so the interior of each rank is 64 × 128 × 1 = 8192 points. The range Range3D of cfl covers only the interior grid points, such as from i_margin to i_total - i_margin.
The number of blocks in stage 1 is (8192 + 255) / 256 = 32. After stage 1, 32 values remain in buf0. In stage 2, for cur_n = 32 it launches with g = (32 + 255) / 256 = 1 block, which leaves 1 value. Therefore the launch in stage 2 and after happens once, and the whole reduce launches kernels twice. Finally, one value is copied to the host, and the values of the ranks are turned into the global minimum with MPI_Allreduce and MPI_MIN (Chapter 8).
7. Synchronization and Finding Errors
7.1 Kernel launches do not wait
When you launch a kernel with <<< >>>, the host moves on to the next line without waiting for the kernel to finish. This is a behavior in which you hand the work to the GPU and the CPU goes ahead. This property is called asynchronous.
Work for the GPU goes into a queue called a stream, and within a stream it is executed in the order it was put in. The fourth value of <<< >>> and the last argument of cudaMemcpyAsync specify the stream, and if omitted the default stream is used. MISO’s for_each and reduce both use the default stream by default, so the order is preserved when the next reduce reads the result of a for_each.
Before the host uses a result from the GPU, it must wait for the work to complete (synchronize).
Function |
What it waits for |
|---|---|
|
It copies after the work already put on the default stream has finished, and waits until the copy is done |
|
All work on that GPU |
|
The work on the specified stream |
|
Does not wait. It only puts the copy on the stream |
In the example of Section 4.3, we could call f_host.copy_from(f) right after for_each because the cudaMemcpy inside waits for the kernel to finish. reduce waits with cudaStreamSynchronize at the end before returning the value, so the return value can be used immediately.
copy_from also has an asynchronous version that takes a stream as an argument (using cudaMemcpyAsync inside), and in that case the caller is responsible for synchronization.
MISO calls cudaDeviceSynchronize explicitly before handing a GPU result to something other than the GPU. In the halo exchange of Section 8, it is called after the kernel that packs the data to send and before MPI_Isend. MPI does not know about CUDA streams, so we need to wait here so that a buffer that has not finished packing is not sent.
7.2 Errors are reported later
Being asynchronous also affects how errors are found. An error that occurs inside a kernel (such as an out-of-range read or write) has not yet happened at the time of the launch. It is reported as the return value when the host calls a synchronizing CUDA function after the kernel has run.
How to check |
What it tells you |
|---|---|
|
Failure of the launch itself (the block is too large, there are no instructions for the GPU, and so on) |
The return values of later calls such as |
Errors that occurred while the kernel was running |
This is why, in the example of Section 5.1, the location of the error was displayed as array1d.hpp (the cudaMemcpy in copy_from). What was wrong was the lambda of apply_wrong, but it became known only when the host next synchronized.
MISO checks the return values of CUDA functions with MISO_CUDA_CHECK (miso/include/miso/cuda_util.cuh).
#define MISO_CUDA_CHECK(ans) \
do { miso::cuda::check_error((ans), __FILE__, __LINE__); } while (0)
If there is an error, check_error displays CUDA Error: <message> <file> <line>, and if MPI has been initialized it stops all ranks with MPI_Abort. __FILE__ and __LINE__ are macros that are replaced by the file name and line number of that line. Note that what is displayed is “where the error was noticed”, which is not necessarily “where the error occurred”.
To pin down where it occurred, run with compute-sanitizer, which is included in the CUDA Toolkit. It displays which thread of which kernel made an illegal read or write. How to use it, and the display for the example of Section 5.1, are summarized in GPU Troubleshooting.
8. Selecting a GPU and MPI
8.1 Choosing a GPU for each rank
When a computer has multiple GPUs, you decide which GPU each MPI rank uses. In MISO, the Env created at the start of main does this (miso/include/miso/env.hpp).
struct Env {
mpi::Env mpi_env;
#ifdef __CUDACC__
cuda::Env cuda_env;
#endif
In a CUDA build, cuda::Env is created after MPI is initialized by mpi::Env. This is the main part of the constructor of cuda::Env.
cudaGetDeviceCount(&internal::device_count);
if (internal::device_count == 0) {
throw std::runtime_error("No CUDA-capable devices were detected.");
}
internal::device_id = mpi::local_rank() % internal::device_count;
cudaSetDevice(internal::device_id);
mpi::local_rank() is the rank number within the same computer (node). mpi::Env obtains it by creating, with MPI_Comm_split_type, a communication group of only the processes on the same node. The remainder of dividing the rank number within the node by the number of GPUs becomes the GPU number. When we ran three ranks on a computer with two GPUs and displayed cuda::device_id() of each rank, we got the following.
rank 0: local rank 0, device 0 of 2
rank 1: local rank 1, device 1 of 2
rank 2: local rank 2, device 0 of 2
Rank 2 shares the same GPU with rank 0, so normally you match the number of ranks per node to the number of GPUs.
8.2 Sending GPU memory with MPI
In the halo exchange of Chapter 8, we packed the values to send to the neighboring rank into a buffer and sent and received them with MPI_Isend and MPI_Irecv. In the CUDA version, the arrays of physical quantities are on the device, so the packing is also done by a kernel. This is a part of the x direction of HaloExchanger<Real, backend::CUDA>::apply in miso/include/miso/mhd_halo_exchange.hpp.
if (mpi_shape.x_procs_pos != MPI_PROC_NULL && grid.i_total > 1) {
pack_x_send<<<cu_shape.grid_dim_x_margin, cu_shape.block_dim>>>(
buff.view(), qq_trgt.view(), grid.view(), Face::Pos);
MISO_CUDA_CHECK(cudaGetLastError());
MISO_CUDA_CHECK(cudaDeviceSynchronize());
MPI_Isend(buff.send_x_pos,
grid.i_margin * grid.j_total * grid.k_total * n_fields,
mpi::data_type<Real>(), mpi_shape.x_procs_pos, 100,
mpi_shape.cart_comm, &reqs[req_count++]);
The flow is as follows.
The
pack_x_sendkernel packs the edge values of the nine physical quantities into the send bufferbuff.send_x_pos.cudaDeviceSynchronize()waits for the packing to finish (7.1).Pass the address of the buffer in device memory as it is to
MPI_IsendandMPI_Irecv.After waiting for the sends and receives to complete with
MPI_Waitall, theunpack_x_recvkernel writes the values of the receive buffer into the ghost cells.
The send and receive buffers (Buffers<Real, backend::CUDA> in miso/include/miso/mhd_buffers.hpp) are device memory allocated with cudaMalloc. An MPI that can send and receive directly from device addresses is called CUDA-aware MPI. MISO’s CUDA version assumes CUDA-aware MPI and has no path that goes through host buffers. To run multiple ranks on GPUs, use an MPI built with CUDA support (such as the Open MPI included in the NVIDIA HPC SDK).
Waiting for the GPU computation to finish (cudaDeviceSynchronize) and waiting for the MPI communication to finish (MPI_Waitall) are different things. Before sending, you wait for the GPU to finish; before using received values, you wait for MPI to finish.
8.3 Running
As with the CPU build, a GPU build is run with mpiexec. The app_run.sh of MISO’s demos loads the environment variables of demo/shared/openmpi_config.env and then runs with 2 ranks.
set -a
source /path/to/MISO/demo/shared/openmpi_config.env
set +a
mpiexec -np 2 ./build/mhd2d_vortex --config=config.yaml
When mhd2d_vortex is run with 2 GPUs and 2 ranks, the time and the output number are displayed in the same format as the CPU build, and the run proceeds to tend: 0.8.
9. Building MISO for a GPU
9.1 The CUDA Toolkit and the GPU generation
To build for a GPU, in addition to the CPU tools (Chapter 1), you need NVIDIA’s CUDA Toolkit (nvcc and the CUDA run-time library). A Mac has no NVIDIA GPU, so GPU builds are done on a Linux computer with an NVIDIA GPU.
nvcc generates different instructions for each GPU generation (compute capability). In CMake, you give CMAKE_CUDA_ARCHITECTURES the compute capability with the decimal point removed.
Example GPUs |
compute capability |
|
|---|---|---|
V100 |
7.0 |
|
A100 |
8.0 |
|
RTX 4090、L40 |
8.9 |
|
H100 |
9.0 |
|
You can find the value for the GPU at hand with the following command (for a reasonably recent driver).
nvidia-smi --query-gpu=name,compute_cap --format=csv
MISO’s miso/CMakeLists.txt sets CMAKE_CUDA_ARCHITECTURES to native (matched to the GPU of the computer doing the build) if you do not specify it. If the computer doing the build has no GPU (for example, when building on a login node and running on a compute node), specify the value for the GPU you will use. As the comment says, this setting must be decided before enable_language(CUDA). If it is decided afterward, CMake uses the default old value first.
If you build with a value that does not match the GPU you use, the build succeeds but the kernel launch fails at run time (GPU Troubleshooting).
9.2 Unit tests and demos
To build MISO’s unit tests for a GPU, add -DMISO_USE_CUDA=ON. The *.cu tests in miso/tests/ (such as test_execution.cu) are built only in this case.
cd /path/to/MISO
cmake -B build_gpu -S miso -DMISO_USE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=89
cmake --build build_gpu -j
(cd build_gpu && ctest --output-on-failure)
For the demos, use -DUSE_CUDA=ON (the CMake configuration of Section 4.4).
cmake -B demo/mhd2d_vortex/build -S demo/mhd2d_vortex -DUSE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=89
cmake --build demo/mhd2d_vortex/build -j
./demo/mhd2d_vortex/app_run.sh
The demo’s app_run.sh uses the executable inside build/. When switching from the CPU version to the GPU version, delete build/ and rebuild so that the previous settings do not remain in the CMake cache (app_clean.sh deletes the output data/ and figs/ along with build/).
10. Running a Program Made on a CPU on a GPU
Build and run the array computation Practice from Chapter 7 (vector addition, reduction, and one-dimensional advection) in a build folder for the GPU. Add the CMake configuration of Section 4.4 to your own CMakeLists.txt. Switch Backend with the #ifdef __CUDACC__ of the entrance to the Practice (switching with the USE_CUDA of Section 4.3 gives the same result). When you rewrite it, check the following points.
The lambdas passed to
for_eachandreduceare written withMISO_LAMBDA, and[&]is not used.The arrays that a lambda captures are Views made from arrays of
Backend. The owning arrays themselves are not captured.The names used inside a lambda are local variables or arguments. When written inside a class, member variables are not used directly (Section 5).
The
initofreduceis the identity element ofop(6.4).Before displaying the result, it has been copied with
copy_frominto an array ofbackend::Host.
Give the CPU and the GPU the same input, and compare each element of the vector, the total, and the distribution after advection. The values of the elements often match, but the last digit of the total may change because of the difference in the order of addition (6.5).
If you do not have a GPU environment yet, add the comparison later and go on to Chapter 10. Build the serial and MPI versions of heat conduction on the CPU, and add the GPU version later using what you learned in this chapter.