From 7b28eff2a950b3615d0ae0630e7c446e9a0ebdb2 Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Sun, 5 Jul 2026 12:52:44 -0700 Subject: [PATCH] Add GPU support for projection on basis in gyrokinetic app MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Introduced device-callable function pointers for evaluation and computational-to-physical (c2p) mappings in the projection on basis functionality. - Enhanced `gkyl_proj_on_basis` structure to include GPU-related fields and logic for managing device memory. - Implemented CUDA kernels for evaluating projections on the GPU, including memory management for device arrays. - Updated the `gk_species_projection` functions to handle projections on both host and device, allowing for seamless integration of GPU acceleration. - Added device-specific implementations for position and velocity mappings to ensure compatibility with GPU execution. - Refactored existing code to support device memory allocation and copying, ensuring efficient data handling between host and device. - Introduced new private headers for managing device contexts related to c2p mappings, ensuring encapsulation and clarity in the codebase. # Plan: `proj_on_basis` on GPUs and GPU projection in `gk_species_projection_calc_proj_func` **Status**: implemented, NOT tested on a GPU (no GPU available on the development machine). CPU builds compile and pass existing tests. All CUDA code follows established in-repo patterns; the "Testing on a GPU machine" section below lists exactly what to run once a GPU is available. ## Goal Port `core/zero/proj_on_basis.c` to GPUs so projections of user-specified functions (antenna sources, BGK operators, time-dependent sources, etc.) can run directly on the device, instead of projecting on the host and copying to the GPU. Then use it to perform the `proj_on_basis` update in `gk_species_projection_calc_proj_func` (gyrokinetic app, `GKYL_PROJ_FUNC` projections) on the GPU. ## Key design facts (established from the codebase) 1. **Only `.cu` files (and the generated-kernel directories) are compiled as CUDA.** With `CC=nvcc`, `NVCC_FLAGS = -x cu -dc -rdc=true` is applied to `%.cu` rules and to the kernel-directory `%.c` rules; every other `.c` file (zero/, apps/, unit/, creg/) is compiled by nvcc as *host C*. In those host TUs `__NVCC__` is still defined, so `gkyl_util.h` includes `cuda_runtime.h` (making `GKYL_HAVE_CUDA` and the host CUDA API available) and `GKYL_CU_DH`/`GKYL_CU_D` expand to attributes the host compiler ignores. Consequences: - `GKYL_CU_DH static inline` functions in *headers* are safe everywhere: plain inline functions in host TUs, device-callable when the header is included from a `.cu` TU. - `__global__` kernels, `<<<>>>` launches, and taking a *device* function address can only appear in `.cu` files. Device functions used by GPU projections must therefore live in headers (as `GKYL_CU_DH` inlines) or in `.cu` translation units — not in plain `.c` app/input files. - Helpers like `gkyl_range_idx`, `gkyl_sub_range_inv_idx`, `gkyl_rect_grid_cell_center`, `gkyl_array_(c)fetch` are header inlines, already usable in kernels. 2. **Device function pointers cannot be produced by host code.** Taking `&func` on the host yields the host address. The repo's existing pattern (see `dg_lbo_gyrokinetic_drag_cu.cu`: "Doing function pointer stuff in here avoids troublesome cudaMemcpyFromSymbol") is to assign function pointers *inside a `__global__` kernel*, where the function name resolves to its device address. For *user-supplied* functions the library cannot name the function at compile time, so we provide a macro (`GKYL_DEFINE_CU_DEV_FUNC_GETTER`, in `gkyl_util.h`, available only under `__CUDACC__`) that is instantiated in a `.cu` file next to the function; it defines a tiny setter kernel plus a host getter that returns the device address (expose it to C code with an `extern "C"` wrapper). 3. **The context must live on the GPU.** The library cannot deep-copy an opaque `void *ctx`; the user must place their context in device memory (`gkyl_cu_malloc` + `gkyl_cu_memcpy`) and pass that pointer. 4. **Updater device clones.** The standard pattern (e.g. `velocity_map_cu.cu`) is: fill the host struct, create device mirrors of its arrays, make a temporary host copy of the struct whose pointers are replaced by `array->on_dev` pointers, `gkyl_cu_memcpy` it into a `gkyl_cu_malloc`ed clone, and store the clone as `up->on_dev` (with `on_dev = up` itself on CPU). ## Changes ### 1. Core: `gkyl_proj_on_basis` GPU support - **`core/zero/gkyl_proj_on_basis.h`** - `struct gkyl_proj_on_basis_inp` gains `bool use_gpu`. When true: - `eval` (and `c2p_func`, if provided) must be **device** function pointers (obtained via `GKYL_DEFINE_CU_DEV_FUNC_GETTER`). - `ctx` / `c2p_func_ctx` must point to device-resident memory. - `c2p_func == 0` still means the identity map (set up on the device internally). - `gkyl_proj_on_basis_new()` keeps its signature (CPU only, `use_gpu = false`). - **`core/zero/gkyl_proj_on_basis_priv.h`** (new) - Holds the `struct gkyl_proj_on_basis` definition (moved out of the `.c` so the `.cu` file can see it), the `log_to_comp` helper (now `GKYL_CU_DH`), and declarations of the CUDA-side functions (`gkyl_proj_on_basis_cu_dev_new`, `gkyl_proj_on_basis_advance_cu`). - Struct gains: `bool use_gpu`, device mirrors `ordinates_cu`, `weights_cu`, `basis_at_ords_cu`, and `struct gkyl_proj_on_basis *on_dev` (self on CPU). - The host-side `ordinates`/`weights`/`basis_at_ords` arrays are retained in all cases so `gkyl_proj_on_basis_quad()` and `gkyl_proj_on_basis_fetch_ordinate()` keep working on the host. - **`core/zero/proj_on_basis.c`** - `_inew`: quadrature setup unchanged (computed on the host); when `use_gpu`, create the `*_cu` device mirrors (copy host data over) and call `gkyl_proj_on_basis_cu_dev_new()` to build the device clone. The identity c2p on GPU is installed by a `<<<1,1>>>` setter kernel that assigns the device-side identity function and points its ctx at the clone's own `grid` member. - `_advance`: dispatches to `gkyl_proj_on_basis_advance_cu()` when `use_gpu` (guarded by `#ifdef GKYL_HAVE_CUDA`). - `_release`: releases device mirrors and `gkyl_cu_free`s the clone. - **`core/zero/proj_on_basis_cu.cu`** (new) - `proj_on_basis_c2p_identity_cu` (`GKYL_CU_D`) + setter kernel. - `gkyl_proj_on_basis_cu_dev_new()`: device clone per the standard pattern. - `gkyl_proj_on_basis_advance_cu()` + kernel: **one thread per cell** of the update range (`range->nblocks, range->nthreads` launch specifiers, standard stride loop with `gkyl_sub_range_inv_idx`). Each thread mirrors the CPU per-cell algorithm exactly (same loop order, hence same FP behavior up to FMA contraction): cell center → for each quadrature node: logical→comp coords, `c2p`, `eval` → accumulate `w_q * f_q * psi_k` into the output cell. - The per-cell function values `f_q` (`num_ret_vals` doubles) are staged in a device scratch array with one slot per cell, allocated/released per advance call — this mirrors the CPU code (which also allocates `fun_at_ords` per call) and avoids any hard cap on `num_ret_vals`. - **`core/zero/gkyl_util.h`** - New macro (under `#ifdef __CUDACC__`, i.e. `.cu` files only): ```c GKYL_DEFINE_CU_DEV_FUNC_GETTER(func, func_type, getter) ``` defines `static func_type getter(void)` returning the device address of `func` by launching a `<<<1,1>>>` kernel that stores `func` (resolved on the device) into device memory and copying it back. Usage sketch — in a `.cu` file (e.g. a companion to an input file, or a library `.cu`): ```c GKYL_CU_DH static void eval_source(double t, const double *xn, double *fout, void *ctx) { ... } GKYL_DEFINE_CU_DEV_FUNC_GETTER(eval_source, evalf_t, eval_source_getter); extern "C" evalf_t eval_source_dev_ptr(void) { return eval_source_getter(); } ``` and from C code: ```c .eval = eval_source_dev_ptr(), .ctx = ctx_on_gpu, .use_gpu = true, ``` Because plain `.c` input files are host TUs even in GPU builds, a C regression test that wants GPU projection needs its device eval function in such a companion `.cu` (or in an app/library `.cu`). This is the "the user will have to ensure the context/functions are on the GPU" part of the task. ### 2. Velocity map: device-callable c2p evaluation - **`vlasov/zero/gkyl_velocity_map.h`**: new `GKYL_CU_DH static inline` function `gkyl_velocity_map_eval_c2p_dev()`. Identical logic to the host `gkyl_velocity_map_eval_c2p()`, but it uses the `vmap` array and the `vmap_basis` pointer (instead of the `vmap_ho`/`vmap_basis_ho` host copies). On the device clone (`gvm->on_dev`) those members are device pointers (`velocity_map_cu.cu` builds the clone that way), so this function is valid in kernels; on CPU builds it also works with the host object. ### 3. Position map: device copy + device-callable evaluation - **`gyrokinetic/zero/gkyl_position_map.h`** - `struct gkyl_position_map` gains `struct gkyl_array *mc2nu_dev` (device mirror of the mapping's DG coefficients) and `struct gkyl_position_map *on_dev` (self on CPU). - `gkyl_position_map_eval_mc2nu()` moves into the header as `GKYL_CU_DH static inline` (VLAs replaced with `GKYL_MAX_CDIM` arrays, since VLAs are not valid in device code). Given the host object it behaves as before; given `gpm->on_dev` inside a kernel it evaluates on the device (the clone's `basis` holds device `eval_expand` etc. pointers, and its `mc2nu` points at the device mirror). - New `gkyl_position_map_make_cu_dev()`: creates the mirror + clone. The clone's `basis` member is overwritten with a device-to-device copy from a `gkyl_cart_modal_{serendip,tensor}_cu_dev_new()` basis so its function pointers are device addresses (the temporary device basis container is released afterwards; the code addresses stay valid). Host-only members (`maps`, `ctxs`, `bmag_ctx`, `constB_ctx`, `xpt_ctx`) are nulled in the clone. No kernels are needed, so this lives in `position_map.c` under `#ifdef GKYL_HAVE_CUDA` (CPU builds get an `assert(false)` stub so the symbol always links). - `gkyl_position_map_set_mc2nu()` also syncs the device mirror when present. - `gkyl_position_map_free()` releases the mirror and clone. - **`gyrokinetic/apps/gyrokinetic.c`**: after `gkyl_position_map_set_mc2nu(...)` (post-geometry), call `gkyl_position_map_make_cu_dev(app->position_map)` when `app->use_gpu`, so every species' projection can use `app->position_map->on_dev`. ### 4. Gyrokinetic app: `GKYL_PROJ_FUNC` on the GPU - **`gyrokinetic/apps/gkyl_gyrokinetic.h`**: `struct gkyl_gyrokinetic_projection` (the `func`/`ctx_func` branch) gains: ```c void (*func_on_dev)(double t, const double *xn, double *fout, void *ctx); void *ctx_func_on_dev; ``` When running on GPUs, the user may set `func_on_dev` to the *device* pointer of their (device-callable) function and `ctx_func_on_dev` to a GPU-resident context. The projection then runs fully on the device. If `func_on_dev` is not set, the app falls back to the previous behavior (project on host, copy to device), so all existing input files keep working unchanged. - **`gyrokinetic/apps/gkyl_gyrokinetic_priv.h`**: `struct gk_proj` (FUNC branch) gains `struct gk_proj_on_basis_c2p_func_ctx *proj_on_basis_c2p_ctx_dev`, the device copy of the c2p context. - **`gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h`** (new): holds `struct gk_proj_on_basis_c2p_func_ctx` (moved out of `gkyl_gyrokinetic_priv.h`, which now includes this header) plus the declaration of the device-address getter, so the small `.cu` below does not need to compile the whole app private header as C++. - **`gyrokinetic/apps/gk_species_projection_cu.cu`** (new): the device c2p function and its address getter: ```c GKYL_CU_D static void gk_proj_on_basis_c2p_phase_func_cu(xcomp, xphys, ctx) -> gkyl_position_map_eval_mc2nu(ctx->pos_map /* on_dev */, ...) gkyl_velocity_map_eval_c2p_dev(ctx->vel_map /* on_dev */, ...) extern "C" proj_on_basis_c2p_t gk_species_projection_c2p_phase_func_cu_dev_ptr(void); ``` - **`gyrokinetic/apps/gk_species_projection.c`** - `gk_species_projection_init` (FUNC branch): if `app->use_gpu && inp.func_on_dev`, build a device-resident c2p context `{cdim, vdim, pos_map->on_dev, vel_map->on_dev}` and create the `gkyl_proj_on_basis` with `use_gpu = true`; otherwise keep the host updater + `proj_host` staging array. Per the app design rules, the conditional is resolved once at init by choosing between two `projection_calc` implementations: - `gk_species_projection_calc_proj_func` — advances directly into `f` (CPU builds, and GPU builds with a device-side updater); - `gk_species_projection_calc_proj_func_host` — advances into `proj_host` and copies to the device (GPU fallback). Both then multiply by the bmag (gyrocenter) and velocity-space Jacobians on whatever device `f` lives on (those ops already support GPUs). - `gk_species_projection_release`: release `proj_host` only if allocated; free the device c2p context if allocated. ### 5. Unit tests - **`core/unit/ctest_proj_on_basis_cu.cu`** (new): the device eval/c2p functions and `extern "C"` getters for their device addresses (unit `.cu` files are globbed into the library in GPU builds, like `ctest_array_cu.cu`). - **`core/unit/ctest_proj_on_basis.c`**: GPU tests (declared/registered under `#ifdef GKYL_HAVE_CUDA`), each comparing the device projection coefficient-by-coefficient (1e-14) against the host projection of the matching host function: - `test_1_cu`: 1x p1 projection of `x^2`, no ctx, identity c2p (set on the device internally). - `test_2d_2c_cu`: 2x p2, `num_ret_vals = 2`, device ctx allocated with `gkyl_cu_malloc` — exercises multi-component projection and device contexts. - `test_c2p_1d_cu`: projection with a nontrivial device c2p mapping vs. the CPU projection with the equivalent host c2p — the same machinery the gyrokinetic app uses. ## Out of scope (follow-ups) - **`eval_on_nodes`** (the nodal evaluator): the same treatment applies — mirror `basis_at_nodes` on the device, clone the updater, evaluate the device `eval` fptr at nodes and do the nodal→modal transform per cell in a kernel. Not needed for `gk_species_projection_calc_proj_func`; deferred to keep this change reviewable. - **Neutral species** (`gk_neut_species_projection.c`): same recipe as `gk_species_projection.c` once this lands. - **Other proj_on_basis app call sites** (Maxwellian moments projections, etc.): these are init-only and can migrate incrementally by supplying device fptrs/ctxs. - **Lua layer**: Lua-defined functions cannot run on the device; Lua apps keep the host+copy path. ## Risks / review notes - **`cudaFree` after kernel launch** (scratch array in `advance_cu`): safe — `cudaFree` synchronizes with prior work on the default stream (and the repo frees kernel-visible memory this way elsewhere). - **Device basis embedded by value** (position map clone): the D2D copy of `struct gkyl_basis` copies device *code* addresses; releasing the temporary container does not invalidate them. - **FP differences**: GPU FMA contraction can make GPU results differ from CPU at the last-ulp level; the unit tests compare to 1e-14, matching other GPU tests in the repo. - **`num_quad`**: quadrature tables are precomputed on the host and copied; no device-side Gauss-Legendre generation needed. ## Testing on a GPU machine ```bash # configure with a machine file that sets CC=nvcc, then: make core-unit -j 8 && ./cuda-build/core/unit/ctest_proj_on_basis make gyrokinetic -j 8 # any gyrokinetic C regression test with GKYL_PROJ_FUNC ICs, adapted to set # .func_on_dev/.ctx_func_on_dev, run with -g; compare frame 0 to the CPU run: ./cuda-build/gyrokinetic/creg/rt_gk_sheath_2x2v_p1 -g -s1 compute-sanitizer --tool memcheck --leak-check full ./cuda-build/core/unit/ctest_proj_on_basis ``` On CPU-only machines: `make core-unit && ./build/core/unit/ctest_proj_on_basis` plus valgrind must remain clean (GPU code paths compile out). --- core/unit/ctest_proj_on_basis.c | 220 +++++++++++++++++- core/unit/ctest_proj_on_basis_cu.cu | 64 +++++ core/zero/gkyl_proj_on_basis.h | 11 +- core/zero/gkyl_proj_on_basis_priv.h | 68 ++++++ core/zero/gkyl_util.h | 31 +++ core/zero/proj_on_basis.c | 85 ++++--- core/zero/proj_on_basis_cu.cu | 123 ++++++++++ gyrokinetic/apps/gk_species_projection.c | 111 ++++++--- gyrokinetic/apps/gk_species_projection_cu.cu | 40 ++++ .../apps/gkyl_gk_proj_on_basis_c2p_priv.h | 29 +++ gyrokinetic/apps/gkyl_gyrokinetic.h | 14 +- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 11 +- gyrokinetic/apps/gyrokinetic.c | 8 + gyrokinetic/zero/gkyl_position_map.h | 57 ++++- gyrokinetic/zero/position_map.c | 77 ++++-- vlasov/zero/gkyl_velocity_map.h | 44 ++++ 16 files changed, 890 insertions(+), 103 deletions(-) create mode 100644 core/unit/ctest_proj_on_basis_cu.cu create mode 100644 core/zero/gkyl_proj_on_basis_priv.h create mode 100644 core/zero/proj_on_basis_cu.cu create mode 100644 gyrokinetic/apps/gk_species_projection_cu.cu create mode 100644 gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h diff --git a/core/unit/ctest_proj_on_basis.c b/core/unit/ctest_proj_on_basis.c index ae780a95af..4ba8550589 100644 --- a/core/unit/ctest_proj_on_basis.c +++ b/core/unit/ctest_proj_on_basis.c @@ -1,9 +1,11 @@ #include +#include #include #include #include #include +#include #include void evalFunc(double t, const double *xn, double* restrict fout, void *ctx) @@ -380,11 +382,221 @@ test_3_3d() gkyl_array_release(distf); } +// Cuda specific tests +#ifdef GKYL_HAVE_CUDA + +// Getters for device function addresses, defined in ctest_proj_on_basis_cu.cu. +// The device functions must match their host counterparts in this file. +evalf_t ctest_proj_on_basis_f_1d_cu_dev_ptr(void); +evalf_t ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void); +proj_on_basis_c2p_t ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void); + +// Context for the 2d two-component function; must match the definition in +// ctest_proj_on_basis_cu.cu. +struct ctest_proj_on_basis_2d_ctx { + double c0, c1; +}; + +// Host counterpart of the device function ctest_pob_f_2d_2c. +void evalFunc_2d_2c(double t, const double *xn, double* restrict fout, void *ctx) +{ + struct ctest_proj_on_basis_2d_ctx *tctx = ctx; + double x = xn[0], y = xn[1]; + fout[0] = tctx->c0 + x*y; + fout[1] = tctx->c1*x*x + y; +} + +// Host counterpart of the device function ctest_pob_c2p_1d. +void c2pFunc_1d(const double *xcomp, double *xphys, void *ctx) +{ + xphys[0] = 0.5*xcomp[0] + 0.1; +} + +// Compare a device array against a host reference, coefficient by coefficient. +static void +check_same_dev_ho(const struct gkyl_array *ref_ho, struct gkyl_array *arr_cu) +{ + struct gkyl_array *arr_ho = gkyl_array_new(GKYL_DOUBLE, arr_cu->ncomp, arr_cu->size); + gkyl_array_copy(arr_ho, arr_cu); + + for (size_t i=0; isize; ++i) { + const double *ref_c = gkyl_array_cfetch(ref_ho, i); + const double *arr_c = gkyl_array_cfetch(arr_ho, i); + for (size_t k=0; kncomp; ++k) { + TEST_CHECK( gkyl_compare(ref_c[k], arr_c[k], 1e-14) ); + TEST_MSG("cell %zu coeff %zu | Expected: %.13e | Produced: %.13e", i, k, ref_c[k], arr_c[k]); + } + } + + gkyl_array_release(arr_ho); +} + +void +test_1_cu() +{ + int poly_order = 1; + double lower[] = {-2.0}, upper[] = {2.0}; + int cells[] = {2}; + struct gkyl_rect_grid grid; + gkyl_rect_grid_init(&grid, 1, lower, upper, cells); + + struct gkyl_basis basis; + gkyl_cart_modal_serendip(&basis, 1, poly_order); + + int nghost[GKYL_MAX_DIM] = { 0 }; + struct gkyl_range arr_range, arr_ext_range; + gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range); + + // Reference: project on the host. + gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis, + poly_order+1, 1, evalFunc, NULL); + struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf); + + // Project the same function on the device. + gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &grid, + .basis = &basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = poly_order+1, + .num_ret_vals = 1, + .eval = ctest_proj_on_basis_f_1d_cu_dev_ptr(), + .ctx = NULL, + .use_gpu = true, + } + ); + struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu); + + check_same_dev_ho(distf, distf_cu); + + gkyl_proj_on_basis_release(projDistf); + gkyl_proj_on_basis_release(projDistf_cu); + gkyl_array_release(distf); + gkyl_array_release(distf_cu); +} + +void +test_2d_2c_cu() +{ + int poly_order = 2; + double lower[] = {-2.0,-1.0}, upper[] = {2.0,3.0}; + int cells[] = {4, 3}; + int ndim = sizeof(cells)/sizeof(cells[0]); + struct gkyl_rect_grid grid; + gkyl_rect_grid_init(&grid, ndim, lower, upper, cells); + + struct gkyl_basis basis; + gkyl_cart_modal_serendip(&basis, ndim, poly_order); + + int nghost[GKYL_MAX_DIM] = { 0 }; + struct gkyl_range arr_range, arr_ext_range; + gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range); + + int num_ret_vals = 2; + struct ctest_proj_on_basis_2d_ctx tctx = { .c0 = 0.5, .c1 = 2.0 }; + + // Reference: project on the host. + gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis, + poly_order+1, num_ret_vals, evalFunc_2d_2c, &tctx); + struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE, + num_ret_vals*basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf); + + // Project the same function on the device; the context must be placed + // in device memory by the user. + struct ctest_proj_on_basis_2d_ctx *tctx_cu = gkyl_cu_malloc(sizeof(struct ctest_proj_on_basis_2d_ctx)); + gkyl_cu_memcpy(tctx_cu, &tctx, sizeof(struct ctest_proj_on_basis_2d_ctx), GKYL_CU_MEMCPY_H2D); + + gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &grid, + .basis = &basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = poly_order+1, + .num_ret_vals = num_ret_vals, + .eval = ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(), + .ctx = tctx_cu, + .use_gpu = true, + } + ); + struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, + num_ret_vals*basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu); + + check_same_dev_ho(distf, distf_cu); + + gkyl_proj_on_basis_release(projDistf); + gkyl_proj_on_basis_release(projDistf_cu); + gkyl_array_release(distf); + gkyl_array_release(distf_cu); + gkyl_cu_free(tctx_cu); +} + +void +test_c2p_1d_cu() +{ + int poly_order = 1; + double lower[] = {-2.0}, upper[] = {2.0}; + int cells[] = {8}; + struct gkyl_rect_grid grid; + gkyl_rect_grid_init(&grid, 1, lower, upper, cells); + + struct gkyl_basis basis; + gkyl_cart_modal_serendip(&basis, 1, poly_order); + + int nghost[GKYL_MAX_DIM] = { 0 }; + struct gkyl_range arr_range, arr_ext_range; + gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range); + + // Reference: project on the host with the host c2p mapping. + gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &grid, + .basis = &basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = poly_order+1, + .num_ret_vals = 1, + .eval = evalFunc, + .c2p_func = c2pFunc_1d, + } + ); + struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf); + + // Project on the device with the equivalent device c2p mapping. + gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &grid, + .basis = &basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = poly_order+1, + .num_ret_vals = 1, + .eval = ctest_proj_on_basis_f_1d_cu_dev_ptr(), + .c2p_func = ctest_proj_on_basis_c2p_1d_cu_dev_ptr(), + .use_gpu = true, + } + ); + struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume); + gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu); + + check_same_dev_ho(distf, distf_cu); + + gkyl_proj_on_basis_release(projDistf); + gkyl_proj_on_basis_release(projDistf_cu); + gkyl_array_release(distf); + gkyl_array_release(distf_cu); +} + +#endif + TEST_LIST = { { "test_1", test_1 }, - { "test_2", test_2 }, - { "test_2_2d", test_2_2d }, - { "test_2_3d", test_2_3d }, - { "test_3_3d", test_3_3d }, + { "test_2", test_2 }, + { "test_2_2d", test_2_2d }, + { "test_2_3d", test_2_3d }, + { "test_3_3d", test_3_3d }, +#ifdef GKYL_HAVE_CUDA + { "test_1_cu", test_1_cu }, + { "test_2d_2c_cu", test_2d_2c_cu }, + { "test_c2p_1d_cu", test_c2p_1d_cu }, +#endif { NULL, NULL }, }; diff --git a/core/unit/ctest_proj_on_basis_cu.cu b/core/unit/ctest_proj_on_basis_cu.cu new file mode 100644 index 0000000000..476bad9590 --- /dev/null +++ b/core/unit/ctest_proj_on_basis_cu.cu @@ -0,0 +1,64 @@ +/* -*- c++ -*- */ + +// Device-side functions (and getters for their device addresses) used by the +// GPU tests in ctest_proj_on_basis.c. Each function must match its host +// counterpart in that file exactly. + +extern "C" { +#include +#include + +evalf_t ctest_proj_on_basis_f_1d_cu_dev_ptr(void); +evalf_t ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void); +proj_on_basis_c2p_t ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void); +} + +// Context for the 2d two-component function; must match the definition in +// ctest_proj_on_basis.c. +struct ctest_proj_on_basis_2d_ctx { + double c0, c1; +}; + +GKYL_CU_DH static void +ctest_pob_f_1d(double t, const double *xn, double *fout, void *ctx) +{ + double x = xn[0]; + fout[0] = x*x; +} + +GKYL_CU_DH static void +ctest_pob_f_2d_2c(double t, const double *xn, double *fout, void *ctx) +{ + struct ctest_proj_on_basis_2d_ctx *tctx = (struct ctest_proj_on_basis_2d_ctx *) ctx; + double x = xn[0], y = xn[1]; + fout[0] = tctx->c0 + x*y; + fout[1] = tctx->c1*x*x + y; +} + +GKYL_CU_DH static void +ctest_pob_c2p_1d(const double *xcomp, double *xphys, void *ctx) +{ + xphys[0] = 0.5*xcomp[0] + 0.1; +} + +GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_f_1d, evalf_t, ctest_pob_f_1d_getter); +GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_f_2d_2c, evalf_t, ctest_pob_f_2d_2c_getter); +GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_c2p_1d, proj_on_basis_c2p_t, ctest_pob_c2p_1d_getter); + +extern "C" evalf_t +ctest_proj_on_basis_f_1d_cu_dev_ptr(void) +{ + return ctest_pob_f_1d_getter(); +} + +extern "C" evalf_t +ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void) +{ + return ctest_pob_f_2d_2c_getter(); +} + +extern "C" proj_on_basis_c2p_t +ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void) +{ + return ctest_pob_c2p_1d_getter(); +} diff --git a/core/zero/gkyl_proj_on_basis.h b/core/zero/gkyl_proj_on_basis.h index 878bc6efe6..993995ed4b 100644 --- a/core/zero/gkyl_proj_on_basis.h +++ b/core/zero/gkyl_proj_on_basis.h @@ -1,5 +1,7 @@ #pragma once +#include + #include #include #include @@ -28,6 +30,13 @@ struct gkyl_proj_on_basis_inp { proj_on_basis_c2p_t c2p_func; // Function that transforms a set of ndim // computational coordinates to physical ones. void *c2p_func_ctx; // Context for c2p_func. + + bool use_gpu; // Whether to run the projection on the GPU. If true, 'eval' + // (and 'c2p_func', if provided) must be device function + // pointers (see GKYL_DEFINE_CU_DEV_FUNC_GETTER in gkyl_util.h) + // and 'ctx'/'c2p_func_ctx' must point to GPU-resident memory; + // the user must place their context on the GPU themselves as + // this updater has no way to perform that copy. }; /** @@ -63,7 +72,7 @@ gkyl_proj_on_basis *gkyl_proj_on_basis_new(const struct gkyl_rect_grid *grid, * @param pob Project on basis updater to run * @param tm Time at which projection must be computed * @param update_rng Range on which to run projection. - * @param out Output array + * @param out Output array (a device array if the updater was created with use_gpu=true). */ void gkyl_proj_on_basis_advance(const gkyl_proj_on_basis *pob, double tm, const struct gkyl_range *update_rng, struct gkyl_array *out); diff --git a/core/zero/gkyl_proj_on_basis_priv.h b/core/zero/gkyl_proj_on_basis_priv.h new file mode 100644 index 0000000000..5bedc7a37d --- /dev/null +++ b/core/zero/gkyl_proj_on_basis_priv.h @@ -0,0 +1,68 @@ +// Private header for the proj_on_basis updater. Not for direct use in user code. +#pragma once + +#include +#include +#include +#include + +struct gkyl_proj_on_basis { + struct gkyl_rect_grid grid; + int num_quad; // number of quadrature points to use in each direction + int num_ret_vals; // number of values returned by eval function + evalf_t eval; // function to project + void *ctx; // evaluation context + + int num_basis; // number of basis functions + int tot_quad; // total number of quadrature points + struct gkyl_array *ordinates; // ordinates for quadrature + struct gkyl_array *weights; // weights for quadrature + struct gkyl_array *basis_at_ords; // basis functions at ordinates + + proj_on_basis_c2p_t c2p; // Function transformin comp to phys coords. + void *c2p_ctx; // Context for the c2p mapping. + + bool use_gpu; // Whether the projection runs on the GPU. + // Device copies of the quadrature data (only allocated when use_gpu=true). + struct gkyl_array *ordinates_cu; + struct gkyl_array *weights_cu; + struct gkyl_array *basis_at_ords_cu; + struct gkyl_proj_on_basis *on_dev; // Device clone of this updater (points to itself on CPU). +}; + +GKYL_CU_DH +static inline void +proj_on_basis_log_to_comp(int ndim, const double *eta, + const double * GKYL_RESTRICT dx, const double * GKYL_RESTRICT xc, + double* GKYL_RESTRICT xout) +{ + // Convert logical to computational coordinates. + for (int d=0; d>>(fptr_d); \ + checkCuda(cudaMemcpy(&fptr_ho, fptr_d, sizeof(func_type), \ + cudaMemcpyDeviceToHost)); \ + checkCuda(cudaFree(fptr_d)); \ + return fptr_ho; \ + } +#endif // __CUDACC__ + #else #undef GKYL_HAVE_CUDA diff --git a/core/zero/proj_on_basis.c b/core/zero/proj_on_basis.c index 2d23038175..2f9bcb3938 100644 --- a/core/zero/proj_on_basis.c +++ b/core/zero/proj_on_basis.c @@ -5,25 +5,9 @@ #include #include #include +#include #include -struct gkyl_proj_on_basis { - struct gkyl_rect_grid grid; - int num_quad; // number of quadrature points to use in each direction - int num_ret_vals; // number of values returned by eval function - evalf_t eval; // function to project - void *ctx; // evaluation context - - int num_basis; // number of basis functions - int tot_quad; // total number of quadrature points - struct gkyl_array *ordinates; // ordinates for quadrature - struct gkyl_array *weights; // weights for quadrature - struct gkyl_array *basis_at_ords; // basis functions at ordinates - - proj_on_basis_c2p_t c2p; // Function transformin comp to phys coords. - void *c2p_ctx; // Context for the c2p mapping. -}; - // Identity comp to phys coord mapping, for when user doesn't provide a map. static inline void c2p_identity(const double *xcomp, double *xphys, void *ctx) @@ -47,6 +31,7 @@ gkyl_proj_on_basis_new(const struct gkyl_rect_grid *grid, const struct gkyl_basi .ctx = ctx, .c2p_func = 0, .c2p_func_ctx = NULL, + .use_gpu = false, } ); } @@ -62,8 +47,15 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp) up->eval = inp->eval; up->ctx = inp->ctx; up->num_basis = inp->basis->num_basis; + up->use_gpu = inp->use_gpu; - if (inp->c2p_func == 0) { + if (up->use_gpu) { + // On the GPU eval and c2p_func must be device function pointers; a NULL + // c2p_func selects an identity mapping assigned on the device. + up->c2p = inp->c2p_func; + up->c2p_ctx = inp->c2p_func_ctx; + } + else if (inp->c2p_func == 0) { up->c2p = c2p_identity; up->c2p_ctx = &up->grid; // Use grid as the context since all we need is ndim. } @@ -87,7 +79,7 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp) } else if (inp->qtype == GKYL_GAUSS_LOBATTO_QUAD) { assert( (num_quad > 1) && (num_quad <= gkyl_gauss_max) ); - + // Gauss-Lobatto quadrature memcpy(ordinates1, gkyl_gauss_lobatto_ordinates[num_quad], sizeof(double[num_quad])); memcpy(weights1, gkyl_gauss_lobatto_weights[num_quad], sizeof(double[num_quad])); @@ -105,7 +97,7 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp) int tot_quad = up->tot_quad = qrange.volume; - // create ordinates and weights for multi-D quadrature + // create ordinates and weights for multi-D quadrature up->ordinates = gkyl_array_new(GKYL_DOUBLE, inp->grid->ndim, tot_quad); up->weights = gkyl_array_new(GKYL_DOUBLE, 1, tot_quad); @@ -114,12 +106,12 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp) while (gkyl_range_iter_next(&iter)) { long node = gkyl_range_idx(&qrange, iter.idx); - + // set ordinates double *ord = gkyl_array_fetch(up->ordinates, node); for (int i=0; igrid->ndim; ++i) ord[i] = ordinates1[iter.idx[i]-qrange.lower[i]]; - + // set weights double *wgt = gkyl_array_fetch(up->weights, node); wgt[0] = 1.0; @@ -133,6 +125,25 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp) inp->basis->eval(gkyl_array_fetch(up->ordinates, n), gkyl_array_fetch(up->basis_at_ords, n)); + up->ordinates_cu = 0; + up->weights_cu = 0; + up->basis_at_ords_cu = 0; + up->on_dev = up; // On the CPU the updater points to itself. +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + // Mirror the quadrature data on the device and create the device clone. + up->ordinates_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, inp->grid->ndim, tot_quad); + up->weights_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, 1, tot_quad); + up->basis_at_ords_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, inp->basis->num_basis, tot_quad); + gkyl_array_copy(up->ordinates_cu, up->ordinates); + gkyl_array_copy(up->weights_cu, up->weights); + gkyl_array_copy(up->basis_at_ords_cu, up->basis_at_ords); + up->on_dev = gkyl_proj_on_basis_cu_dev_new(up); + } +#else + assert(!up->use_gpu); +#endif + return up; } @@ -146,15 +157,6 @@ double* gkyl_proj_on_basis_fetch_ordinate(const struct gkyl_proj_on_basis *up, l return gkyl_array_fetch(up->ordinates, node); } -static inline void -log_to_comp(int ndim, const double *eta, - const double * GKYL_RESTRICT dx, const double * GKYL_RESTRICT xc, - double* GKYL_RESTRICT xout) -{ - // Convert logical to computational coordinates. - for (int d=0; duse_gpu) { + gkyl_proj_on_basis_advance_cu(up, tm, update_range, arr); + return; + } +#endif + double xc[GKYL_MAX_DIM], xmu[GKYL_MAX_DIM]; int num_ret_vals = up->num_ret_vals; int tot_quad = up->tot_quad; struct gkyl_array *fun_at_ords = gkyl_array_new(GKYL_DOUBLE, num_ret_vals, tot_quad); - + struct gkyl_range_iter iter; gkyl_range_iter_init(&iter, update_range); - + while (gkyl_range_iter_next(&iter)) { gkyl_rect_grid_cell_center(&up->grid, iter.idx, xc); for (int i=0; igrid.ndim, gkyl_array_cfetch(up->ordinates, i), + proj_on_basis_log_to_comp(up->grid.ndim, gkyl_array_cfetch(up->ordinates, i), up->grid.dx, xc, xmu); up->c2p(xmu, xmu, up->c2p_ctx); up->eval(tm, xmu, gkyl_array_fetch(fun_at_ords, i), up->ctx); @@ -218,5 +227,13 @@ gkyl_proj_on_basis_release(struct gkyl_proj_on_basis* up) gkyl_array_release(up->ordinates); gkyl_array_release(up->weights); gkyl_array_release(up->basis_at_ords); +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + gkyl_array_release(up->ordinates_cu); + gkyl_array_release(up->weights_cu); + gkyl_array_release(up->basis_at_ords_cu); + gkyl_cu_free(up->on_dev); + } +#endif gkyl_free(up); } diff --git a/core/zero/proj_on_basis_cu.cu b/core/zero/proj_on_basis_cu.cu new file mode 100644 index 0000000000..31b0fcc6b2 --- /dev/null +++ b/core/zero/proj_on_basis_cu.cu @@ -0,0 +1,123 @@ +/* -*- c++ -*- */ + +extern "C" { +#include + +#include +#include +#include +#include +#include +#include +#include +} + +// Identity comp to phys coord mapping, for when user doesn't provide a map. +GKYL_CU_D static void +proj_on_basis_c2p_identity_cu(const double *xcomp, double *xphys, void *ctx) +{ + struct gkyl_rect_grid *grid = (struct gkyl_rect_grid *) ctx; + int ndim = grid->ndim; + for (int d=0; dc2p = proj_on_basis_c2p_identity_cu; + up->c2p_ctx = &up->grid; // Use grid as the context since all we need is ndim. +} + +struct gkyl_proj_on_basis* +gkyl_proj_on_basis_cu_dev_new(struct gkyl_proj_on_basis *up) +{ + struct gkyl_proj_on_basis *up_cu = (struct gkyl_proj_on_basis*) + gkyl_cu_malloc(sizeof(struct gkyl_proj_on_basis)); + + // Clone with device pointers to the quadrature data. The eval/c2p members + // and their contexts are device pointers provided by the user. + struct gkyl_proj_on_basis up_ho = *up; + up_ho.ordinates = up->ordinates_cu->on_dev; + up_ho.weights = up->weights_cu->on_dev; + up_ho.basis_at_ords = up->basis_at_ords_cu->on_dev; + up_ho.on_dev = up_cu; + + gkyl_cu_memcpy(up_cu, &up_ho, sizeof(struct gkyl_proj_on_basis), GKYL_CU_MEMCPY_H2D); + + // When the user does not provide a c2p mapping use the identity mapping, + // which must be assigned on the device to obtain a device function pointer. + if (up->c2p == 0) + proj_on_basis_set_c2p_identity_cu_ker<<<1,1>>>(up_cu); + + return up_cu; +} + +__global__ static void +gkyl_proj_on_basis_advance_cu_ker(const struct gkyl_proj_on_basis *up, double tm, + struct gkyl_range update_range, struct gkyl_array* GKYL_RESTRICT fun_at_ords, + struct gkyl_array* GKYL_RESTRICT arr) +{ + int idx[GKYL_MAX_DIM]; + double xc[GKYL_MAX_DIM], xmu[GKYL_MAX_DIM]; + + int num_basis = up->num_basis; + int tot_quad = up->tot_quad; + int num_ret_vals = up->num_ret_vals; + int ndim = up->grid.ndim; + + const double* GKYL_RESTRICT weights = (const double*) up->weights->data; + const double* GKYL_RESTRICT basis_at_ords = (const double*) up->basis_at_ords->data; + + for (unsigned long linc1 = threadIdx.x + blockIdx.x*blockDim.x; + linc1 < update_range.volume; + linc1 += gridDim.x*blockDim.x) + { + gkyl_sub_range_inv_idx(&update_range, linc1, idx); + gkyl_rect_grid_cell_center(&up->grid, idx, xc); + + // Scratch space for the function evaluated at one quadrature node, + // one slot per cell (i.e. per linc1). + double *fq = (double*) gkyl_array_fetch(fun_at_ords, linc1); + + long lidx = gkyl_range_idx(&update_range, idx); + double *f = (double*) gkyl_array_fetch(arr, lidx); + + // Arrangement of f is as: + // c0[0], c0[1], ... c1[0], c1[1], .... + // where c0, c1, ... are components of f (num_ret_vals). + for (int k=0; kordinates, imu), + up->grid.dx, xc, xmu); + up->c2p(xmu, xmu, up->c2p_ctx); + up->eval(tm, xmu, fq, up->ctx); + + long offset = 0; + for (int n=0; nnum_ret_vals, update_range->volume); + + gkyl_proj_on_basis_advance_cu_ker<<nblocks, update_range->nthreads>>>( + up->on_dev, tm, *update_range, fun_at_ords->on_dev, arr->on_dev); + + gkyl_array_release(fun_at_ords); +} diff --git a/gyrokinetic/apps/gk_species_projection.c b/gyrokinetic/apps/gk_species_projection.c index 843e0601e8..fcc6763dd5 100644 --- a/gyrokinetic/apps/gk_species_projection.c +++ b/gyrokinetic/apps/gk_species_projection.c @@ -67,23 +67,36 @@ func_gaussian(double t, const double* xn, double* GKYL_RESTRICT fout, void *ctx) } static void -gk_species_projection_calc_proj_func(gkyl_gyrokinetic_app *app, struct gk_species *s, - struct gk_proj *proj, struct gkyl_array *f, double tm) +gk_species_projection_multiply_jacobians(gkyl_gyrokinetic_app *app, struct gk_species *s, + struct gkyl_array *f) { - if (app->use_gpu) { - gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, proj->proj_host); - gkyl_array_copy(f, proj->proj_host); - } - else { - gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, f); - } // Multiply by the gyrocenter coord jacobian (bmag). - gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, f, - app->gk_geom->geo_corn.bmag, f, &app->local, &s->local); - // Multiply by the velocity-space jacobian. + gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, f, + app->gk_geom->geo_corn.bmag, f, &app->local, &s->local); + // Multiply by the velocity-space jacobian. gkyl_array_scale_by_cell(f, s->vel_map->jacobvel); } +static void +gk_species_projection_calc_proj_func(gkyl_gyrokinetic_app *app, struct gk_species *s, + struct gk_proj *proj, struct gkyl_array *f, double tm) +{ + // The projection runs on the host or on the device (directly into f) + // depending on how proj_func was created. + gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, f); + gk_species_projection_multiply_jacobians(app, s, f); +} + +static void +gk_species_projection_calc_proj_func_host(gkyl_gyrokinetic_app *app, struct gk_species *s, + struct gk_proj *proj, struct gkyl_array *f, double tm) +{ + // Project on the host and copy to the device. + gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, proj->proj_host); + gkyl_array_copy(f, proj->proj_host); + gk_species_projection_multiply_jacobians(app, s, f); +} + static void project_moment_if_needed(bool from_file, struct gkyl_proj_on_basis *proj_op, double tm, const struct gkyl_range *conf_range, struct gkyl_array *arr, double scale_fac) @@ -430,22 +443,59 @@ gk_species_projection_init(struct gkyl_gyrokinetic_app *app, struct gk_species * proj->proj_on_basis_c2p_ctx.vel_map = s->vel_map; proj->proj_on_basis_c2p_ctx.pos_map = app->position_map; if (proj->proj_id == GKYL_PROJ_FUNC) { - proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { - .grid = &s->grid, - .basis = &s->basis, - .qtype = GKYL_GAUSS_QUAD, - .num_quad = app->basis.poly_order+1, - .num_ret_vals = 1, - .eval = inp.func, - .ctx = inp.ctx_func, - .c2p_func = proj_on_basis_c2p_phase_func, - .c2p_func_ctx = &proj->proj_on_basis_c2p_ctx, - } - ); - if (app->use_gpu) - proj->proj_host = mkarr(false, s->basis.num_basis, s->local_ext.volume); + proj->proj_host = 0; + proj->proj_on_basis_c2p_ctx_dev = 0; + + bool proj_on_gpu = false; +#ifdef GKYL_HAVE_CUDA + proj_on_gpu = app->use_gpu && (inp.func_on_dev != 0); + if (proj_on_gpu) { + // Project on the device with the user's device function. The c2p + // context must live on the device and hold the maps' device objects. + struct gk_proj_on_basis_c2p_func_ctx c2p_ctx_dev_ho = { + .cdim = app->cdim, + .vdim = s->local_vel.ndim, + .vel_map = s->vel_map->on_dev, + .pos_map = app->position_map->on_dev, + }; + proj->proj_on_basis_c2p_ctx_dev = gkyl_cu_malloc(sizeof(struct gk_proj_on_basis_c2p_func_ctx)); + gkyl_cu_memcpy(proj->proj_on_basis_c2p_ctx_dev, &c2p_ctx_dev_ho, + sizeof(struct gk_proj_on_basis_c2p_func_ctx), GKYL_CU_MEMCPY_H2D); + + proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &s->grid, + .basis = &s->basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = app->basis.poly_order+1, + .num_ret_vals = 1, + .eval = inp.func_on_dev, + .ctx = inp.ctx_func_on_dev, + .c2p_func = gk_species_projection_c2p_phase_func_cu_dev_ptr(), + .c2p_func_ctx = proj->proj_on_basis_c2p_ctx_dev, + .use_gpu = true, + } + ); + } +#endif + if (!proj_on_gpu) { + proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) { + .grid = &s->grid, + .basis = &s->basis, + .qtype = GKYL_GAUSS_QUAD, + .num_quad = app->basis.poly_order+1, + .num_ret_vals = 1, + .eval = inp.func, + .ctx = inp.ctx_func, + .c2p_func = proj_on_basis_c2p_phase_func, + .c2p_func_ctx = &proj->proj_on_basis_c2p_ctx, + } + ); + if (app->use_gpu) + proj->proj_host = mkarr(false, s->basis.num_basis, s->local_ext.volume); + } - proj->projection_calc = gk_species_projection_calc_proj_func; + proj->projection_calc = (app->use_gpu && !proj_on_gpu) ? + gk_species_projection_calc_proj_func_host : gk_species_projection_calc_proj_func; proj->moms_correct = gk_species_projection_correct_all_moms_none; } else { @@ -478,9 +528,14 @@ gk_species_projection_release(const struct gkyl_gyrokinetic_app *app, const stru { if (proj->proj_id == GKYL_PROJ_FUNC) { gkyl_proj_on_basis_release(proj->proj_func); - if (app->use_gpu) { + if (proj->proj_host) { gkyl_array_release(proj->proj_host); } +#ifdef GKYL_HAVE_CUDA + if (proj->proj_on_basis_c2p_ctx_dev) { + gkyl_cu_free(proj->proj_on_basis_c2p_ctx_dev); + } +#endif } else if (proj->proj_id == GKYL_PROJ_MAXWELLIAN_PRIM || proj->proj_id == GKYL_PROJ_BIMAXWELLIAN) { gkyl_array_release(proj->dens); diff --git a/gyrokinetic/apps/gk_species_projection_cu.cu b/gyrokinetic/apps/gk_species_projection_cu.cu new file mode 100644 index 0000000000..a493517a43 --- /dev/null +++ b/gyrokinetic/apps/gk_species_projection_cu.cu @@ -0,0 +1,40 @@ +/* -*- c++ -*- */ + +extern "C" { +#include +#include +#include +#include +#include +} + +// Comp. to phys. mapping for phase-space projections (position map in +// configuration space, velocity map in velocity space), callable on the +// device. The context must be a GPU-resident gk_proj_on_basis_c2p_func_ctx +// whose pos_map/vel_map members are the maps' device (on_dev) objects. +GKYL_CU_D static void +gk_proj_on_basis_c2p_phase_func_cu(const double *xcomp, double *xphys, void *ctx) +{ + struct gk_proj_on_basis_c2p_func_ctx *c2p_ctx = (struct gk_proj_on_basis_c2p_func_ctx *) ctx; + int cdim = c2p_ctx->cdim; // Assumes update range is a phase range. + gkyl_position_map_eval_mc2nu(c2p_ctx->pos_map, xcomp, xphys); + gkyl_velocity_map_eval_c2p_dev(c2p_ctx->vel_map, &xcomp[cdim], &xphys[cdim]); +} + +__global__ static void +gk_proj_on_basis_c2p_phase_func_set_cu_ker(proj_on_basis_c2p_t *fptr_d) +{ + // Assigned in device code so the function pointer is a device address. + *fptr_d = gk_proj_on_basis_c2p_phase_func_cu; +} + +extern "C" proj_on_basis_c2p_t +gk_species_projection_c2p_phase_func_cu_dev_ptr(void) +{ + proj_on_basis_c2p_t *fptr_d, fptr_ho; + checkCuda(cudaMalloc(&fptr_d, sizeof(proj_on_basis_c2p_t))); + gk_proj_on_basis_c2p_phase_func_set_cu_ker<<<1,1>>>(fptr_d); + checkCuda(cudaMemcpy(&fptr_ho, fptr_d, sizeof(proj_on_basis_c2p_t), cudaMemcpyDeviceToHost)); + checkCuda(cudaFree(fptr_d)); + return fptr_ho; +} diff --git a/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h b/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h new file mode 100644 index 0000000000..c2faa069ef --- /dev/null +++ b/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h @@ -0,0 +1,29 @@ +// Private header for use in the Gyrokinetic app: do not include in +// user-facing header files! Context (and device helpers) for the +// computational-to-physical coordinate mappings passed to proj_on_basis. +#pragma once + +#include +#include +#include + +// Context for c2p function passed to proj_on_basis. +struct gk_proj_on_basis_c2p_func_ctx { + int cdim, vdim; + struct gkyl_position_map *pos_map; + struct gkyl_velocity_map *vel_map; +}; + +#ifdef GKYL_HAVE_CUDA + +/** + * Return the device address of the phase-space c2p function (position map in + * configuration space, velocity map in velocity space), for projecting with + * proj_on_basis on the GPU. The context passed along with it must be a + * GPU-resident gk_proj_on_basis_c2p_func_ctx whose pos_map/vel_map members + * are the maps' device (on_dev) objects. Defined in + * gk_species_projection_cu.cu. + */ +proj_on_basis_c2p_t gk_species_projection_c2p_phase_func_cu_dev_ptr(void); + +#endif diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h index e04e1510ba..a148d5c5d0 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic.h @@ -38,8 +38,18 @@ struct gkyl_gyrokinetic_projection { // Distribution function to project and its context. For fluid neutrals // this function returns mass density, momentum density, and total // energy density. - void (*func)(double t, const double *xn, double *fout, void *ctx); - void *ctx_func; + void (*func)(double t, const double *xn, double *fout, void *ctx); + void *ctx_func; + + // Optionally, when running on GPUs, the device address of a + // device-callable version of 'func' (defined in a CUDA-compiled + // translation unit and obtained with GKYL_DEFINE_CU_DEV_FUNC_GETTER) + // and a device-resident context. When provided, the projection is + // performed on the GPU; the user must place the context on the GPU + // themselves as the app has no way to perform that copy. If absent, + // the function is projected on the host and copied to the GPU. + void (*func_on_dev)(double t, const double *xn, double *fout, void *ctx); + void *ctx_func_on_dev; }; struct { // For Maxwellians (or BiMaxwellians), specify density, parallel speed diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 5615160973..41ebabeeb7 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -715,12 +715,9 @@ struct gk_scaling { struct gk_scaling *sca, int ridx, double tm, int frame); }; -// Context for c2p function passed to proj_on_basis. -struct gk_proj_on_basis_c2p_func_ctx { - int cdim, vdim; - struct gkyl_position_map *pos_map; - struct gkyl_velocity_map *vel_map; -}; +// Context for c2p function passed to proj_on_basis is defined in +// gkyl_gk_proj_on_basis_c2p_priv.h (shared with gk_species_projection_cu.cu). +#include struct gk_proj { enum gkyl_projection_id proj_id; // Type of projection. @@ -730,6 +727,8 @@ struct gk_proj { struct { struct gkyl_proj_on_basis *proj_func; // Projection operator for specified function. struct gkyl_array *proj_host; // Array for projection on host-side if running on GPUs. + struct gk_proj_on_basis_c2p_func_ctx *proj_on_basis_c2p_ctx_dev; // Device copy of the + // c2p context (GPU projection). }; // Maxwellian and Bi-Maxwellian projection from primitive moments. struct { diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c index 4b907cb034..b5cb4f7134 100644 --- a/gyrokinetic/apps/gyrokinetic.c +++ b/gyrokinetic/apps/gyrokinetic.c @@ -416,6 +416,14 @@ gkyl_gyrokinetic_app_new_geom(struct gkyl_gk *gk) gkyl_position_map_set_mc2nu(app->position_map, app->gk_geom->geo_corn.mc2nu_pos); +#ifdef GKYL_HAVE_CUDA + if (app->use_gpu) { + // Device copy of the position map, used e.g. by projections that + // evaluate the c2p mapping inside device kernels. + gkyl_position_map_make_cu_dev(app->position_map); + } +#endif + const struct gkyl_dg_geom_inp dg_geom_inp = { .grid = &app->grid, .range = &app->local_ext, diff --git a/gyrokinetic/zero/gkyl_position_map.h b/gyrokinetic/zero/gkyl_position_map.h index 8d88742e85..e3493566d5 100644 --- a/gyrokinetic/zero/gkyl_position_map.h +++ b/gyrokinetic/zero/gkyl_position_map.h @@ -6,6 +6,9 @@ #include #include #include +#include + +#include enum gkyl_position_map_id { GKYL_PMAP_USER_INPUT = 0, // Function projection. User specified. Default @@ -62,7 +65,11 @@ struct gkyl_position_map { // Stuff for constant B mapping struct gkyl_bmag_ctx *bmag_ctx; // Context for magnetic field calculation struct gkyl_position_map_const_B_ctx *constB_ctx; // Context for constant B mapping - struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping + struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping + + struct gkyl_array *mc2nu_dev; // Device mirror of mc2nu (only with a device copy, see + // gkyl_position_map_make_cu_dev; kept in sync by set_mc2nu). + struct gkyl_position_map *on_dev; // Device copy of itself (points to itself on CPU). }; struct gkyl_position_map_const_B_ctx { @@ -175,16 +182,54 @@ void gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, double zcenter, double w, double psisep); +/** + * Create a device (GPU) copy of the position map object, stored in + * gpm->on_dev, for use inside device kernels with + * gkyl_position_map_eval_mc2nu. Call it after the geometry has set the + * mc2nu array (gkyl_position_map_set_mc2nu keeps the device mirror in sync + * afterwards). Only available in GPU builds. + * + * @param gpm Position map object. + */ +void +gkyl_position_map_make_cu_dev(struct gkyl_position_map* gpm); + /** * Evaluate the position mapping at a specific computational (position) coordinate. - * NOTE: done on the host. + * Callable on the host with the host object, or inside device kernels with the + * device object (gpm->on_dev). On GPU builds do not call it on the host with + * the device object or vice-versa. * * @param gpm Gkyl position map object. - * @param xc Computational position coordinates. - * @param xnu Resulting non-uniform position coordinates. + * @param x_comp Computational position coordinates. + * @param x_fa Resulting non-uniform position coordinates. */ -void -gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *xc, double *xnu); +GKYL_CU_DH static inline void +gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *x_comp, double *x_fa) +{ + int cidx[GKYL_MAX_CDIM]; + for (int i=0; igrid.ndim; i++) { + int idxtemp = gpm->global.lower[i] + (int) floor((x_comp[i] - (gpm->grid.lower[i]) )/gpm->grid.dx[i]); + idxtemp = GKYL_MAX2(GKYL_MIN2(idxtemp, gpm->local.upper[i]), gpm->local.lower[i]); + cidx[i] = idxtemp; + } + long lidx = gkyl_range_idx(&gpm->local, cidx); + const double *pmap_coeffs = (const double *) gkyl_array_cfetch(gpm->mc2nu, lidx); + double cxc[GKYL_MAX_CDIM]; + double x_log[GKYL_MAX_CDIM]; + gkyl_rect_grid_cell_center(&gpm->grid, cidx, cxc); + for (int i=0; igrid.ndim; i++) { + x_log[i] = (x_comp[i]-cxc[i])/(gpm->grid.dx[i]*0.5); + } + double xyz_fa[3]; + for (int i=0; i<3; i++) { + xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i*gpm->basis.num_basis]); + } + for (int i=0; igrid.ndim; i++) { + x_fa[i] = xyz_fa[i]; + } + x_fa[gpm->grid.ndim-1] = xyz_fa[2]; +} /** * Evaluate the slope of the position mapping at a specific computational (position) coordinate. diff --git a/gyrokinetic/zero/position_map.c b/gyrokinetic/zero/position_map.c index 47b0e95230..9a2e70a6c1 100644 --- a/gyrokinetic/zero/position_map.c +++ b/gyrokinetic/zero/position_map.c @@ -9,6 +9,7 @@ #include #include +#include // Remove with the print statements at the bottom #include @@ -39,6 +40,8 @@ gkyl_position_map_null_new() gpm->xpt_ctx = gkyl_malloc(sizeof(struct gkyl_position_map_xpt_ctx)); gpm->bmag_ctx = gkyl_malloc(sizeof(struct gkyl_bmag_ctx)); gpm->bmag_ctx->bmag = gkyl_array_new(GKYL_DOUBLE, 1, 1); + gpm->mc2nu_dev = 0; + gpm->on_dev = gpm; // On the CPU the object points to itself. gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free); for (int i = 0; i < 3; i++){ @@ -161,6 +164,8 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g gpm->basis = basis; gpm->cdim = grid.ndim; gpm->mc2nu = gkyl_array_new(GKYL_DOUBLE, 3*gpm->basis.num_basis, gpm->local_ext.volume); + gpm->mc2nu_dev = 0; + gpm->on_dev = gpm; // On the CPU the object points to itself. gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free); struct gkyl_position_map *gpm_out = gpm; @@ -171,6 +176,8 @@ void gkyl_position_map_set_mc2nu(struct gkyl_position_map* gpm, struct gkyl_array* mc2nu) { gkyl_array_copy(gpm->mc2nu, mc2nu); + if (gpm->mc2nu_dev) + gkyl_array_copy(gpm->mc2nu_dev, gpm->mc2nu); } void @@ -225,31 +232,53 @@ gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, do } } -void -gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *x_comp, double *x_fa) +void +gkyl_position_map_make_cu_dev(struct gkyl_position_map* gpm) { - int cidx[GKYL_MAX_CDIM]; - for(int i = 0; i < gpm->grid.ndim; i++){ - int idxtemp = gpm->global.lower[i] + (int) floor((x_comp[i] - (gpm->grid.lower[i]) )/gpm->grid.dx[i]); - idxtemp = GKYL_MAX2(GKYL_MIN2(idxtemp, gpm->local.upper[i]), gpm->local.lower[i]); - cidx[i] = idxtemp; - } - long lidx = gkyl_range_idx(&gpm->local, cidx); - const double *pmap_coeffs = gkyl_array_cfetch(gpm->mc2nu, lidx); - double cxc[gpm->grid.ndim]; - double x_log[gpm->grid.ndim]; - gkyl_rect_grid_cell_center(&gpm->grid, cidx, cxc); - for(int i = 0; i < gpm->grid.ndim; i++){ - x_log[i] = (x_comp[i]-cxc[i])/(gpm->grid.dx[i]*0.5); - } - double xyz_fa[3]; - for(int i = 0; i < 3; i++){ - xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i*gpm->basis.num_basis]); +#ifdef GKYL_HAVE_CUDA + // Device mirror of the mapping's DG coefficients, kept in sync by set_mc2nu. + gpm->mc2nu_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, gpm->mc2nu->ncomp, gpm->mc2nu->size); + gkyl_array_copy(gpm->mc2nu_dev, gpm->mc2nu); + + // Clone with device pointers; host-only members are nulled. + struct gkyl_position_map gpm_ho = *gpm; + gpm_ho.mc2nu = gpm->mc2nu_dev->on_dev; + for (int i=0; i<3; i++) { + gpm_ho.maps[i] = 0; + gpm_ho.map_derivs[i] = 0; + gpm_ho.ctxs[i] = 0; } - for (int i=0; igrid.ndim; i++) { - x_fa[i] = xyz_fa[i]; + gpm_ho.bmag_ctx = 0; + gpm_ho.constB_ctx = 0; + gpm_ho.xpt_ctx = 0; + + struct gkyl_position_map *gpm_cu = gkyl_cu_malloc(sizeof(struct gkyl_position_map)); + gpm_ho.on_dev = gpm_cu; + gkyl_cu_memcpy(gpm_cu, &gpm_ho, sizeof(struct gkyl_position_map), GKYL_CU_MEMCPY_H2D); + + // Overwrite the basis in the clone with one whose function pointers are + // device addresses (the temporary device container can be released after + // copying its contents; the code addresses it holds remain valid). + struct gkyl_basis *basis_cu; + switch (gpm->basis.b_type) { + case GKYL_BASIS_MODAL_SERENDIPITY: + basis_cu = gkyl_cart_modal_serendip_cu_dev_new(gpm->basis.ndim, gpm->basis.poly_order); + break; + case GKYL_BASIS_MODAL_TENSOR: + basis_cu = gkyl_cart_modal_tensor_cu_dev_new(gpm->basis.ndim, gpm->basis.poly_order); + break; + default: + assert(false); + break; } - x_fa[gpm->grid.ndim-1] = xyz_fa[2]; + gkyl_cu_memcpy((char *)gpm_cu + offsetof(struct gkyl_position_map, basis), basis_cu, + sizeof(struct gkyl_basis), GKYL_CU_MEMCPY_D2D); + gkyl_cart_modal_basis_release_cu(basis_cu); + + gpm->on_dev = gpm_cu; +#else + assert(false); +#endif } void @@ -364,6 +393,10 @@ gkyl_position_map_free(const struct gkyl_ref_count *ref) struct gkyl_position_map *gpm = container_of(ref, struct gkyl_position_map, ref_count); gkyl_array_release(gpm->mc2nu); gkyl_array_release(gpm->bmag_ctx->bmag); + if (gpm->mc2nu_dev) + gkyl_array_release(gpm->mc2nu_dev); + if (gpm->on_dev != gpm) + gkyl_cu_free(gpm->on_dev); if (gpm->to_optimize == true) { gkyl_free(gpm->constB_ctx->theta_extrema); diff --git a/vlasov/zero/gkyl_velocity_map.h b/vlasov/zero/gkyl_velocity_map.h index b7054bfa45..17ce622f9a 100644 --- a/vlasov/zero/gkyl_velocity_map.h +++ b/vlasov/zero/gkyl_velocity_map.h @@ -6,6 +6,9 @@ #include #include #include +#include + +#include typedef void (*mapc2p_t)(double t, const double *zc, double *vp, void *ctx); @@ -115,6 +118,47 @@ gkyl_velocity_map_reduce_dv_range(const struct gkyl_velocity_map* gvm, enum gkyl void gkyl_velocity_map_eval_c2p(const struct gkyl_velocity_map* gvm, const double *zc, double *vp); +/** + * Evaluate the velocity mapping at a specific computational (velocity) + * coordinate. Unlike gkyl_velocity_map_eval_c2p, this variant uses the vmap + * and vmap_basis members (not their _ho host copies), so it is callable + * inside device kernels when given the device object (gvm->on_dev). On GPU + * builds do not call it on the host with the host object, whose vmap data + * lives on the device. + * + * @param gvm Velocity map object (on_dev object inside kernels). + * @param zc Computational velocity coordinates. + * @param vp Resulting physical velocity coordinates. + */ +GKYL_CU_DH static inline void +gkyl_velocity_map_eval_c2p_dev(const struct gkyl_velocity_map* gvm, const double *zc, double *vp) +{ + // Find the index of the cell containing zc. + int idx_zc[GKYL_MAX_VDIM]; + for (int d=0; dlocal_ext_vel.ndim; d++) { + int idx = gvm->local_ext_vel.lower[d] + (int) floor((zc[d] - (gvm->grid_vel.lower[d]) )/gvm->grid_vel.dx[d]); + // Bound idx to the range in the grid. If it falls outside of that is due + // to floating point arithmetic, or due to an error in the code. + idx = GKYL_MIN2(idx, gvm->local_ext_vel.upper[d]); + idx = GKYL_MAX2(idx, gvm->local_ext_vel.lower[d]); + idx_zc[d] = idx; + } + + // Fetch DG coefficients of the velocity map in idx_zc. + long lidx_zc = gkyl_range_idx(&gvm->local_ext_vel, idx_zc); + const double *vmap_c = (const double *) gkyl_array_cfetch(gvm->vmap, lidx_zc); + + double zc_cc[GKYL_MAX_VDIM]; + gkyl_rect_grid_cell_center(&gvm->grid_vel, idx_zc, zc_cc); + + for (int d=0; dlocal_ext_vel.ndim; d++) { + // Convert computational to logical coord. + double zlog[] = {(zc[d] - zc_cc[d]) / (0.5*gvm->grid_vel.dx[d])}; + // Evaluate vmap expansion at logical coord. + vp[d] = gvm->vmap_basis->eval_expand(zlog, &vmap_c[d*gvm->vmap_basis->num_basis]); + } +} + /** * Indicate if this velocity map object is allocated on the GPU. *