1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-tools.h> 11 #include <cuda.h> 12 #include <cuda_runtime.h> 13 #include <stdbool.h> 14 #include <stddef.h> 15 #include <string.h> 16 17 #include "../cuda/ceed-cuda-common.h" 18 #include "../cuda/ceed-cuda-compile.h" 19 #include "ceed-cuda-shared.h" 20 21 //------------------------------------------------------------------------------ 22 // Device initalization 23 //------------------------------------------------------------------------------ 24 int CeedInit_CudaInterp(CeedScalar *d_B, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B); 25 int CeedInit_CudaGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); 26 int CeedInit_CudaCollocatedGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); 27 28 //------------------------------------------------------------------------------ 29 // Apply basis 30 //------------------------------------------------------------------------------ 31 static int CeedBasisApplyTensorCore_Cuda_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, CeedTransposeMode t_mode, 32 CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 33 Ceed ceed; 34 Ceed_Cuda *ceed_Cuda; 35 CeedInt dim, num_comp; 36 const CeedScalar *d_u; 37 CeedScalar *d_v; 38 CeedBasis_Cuda_shared *data; 39 40 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 41 CeedCallBackend(CeedGetData(ceed, &ceed_Cuda)); 42 CeedCallBackend(CeedBasisGetData(basis, &data)); 43 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 44 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 45 46 // Get read/write access to u, v 47 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 48 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 49 if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); 50 else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 51 52 // Apply basis operation 53 switch (eval_mode) { 54 case CEED_EVAL_INTERP: { 55 CeedInt P_1d, Q_1d; 56 57 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 58 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 59 CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 60 61 CeedCallBackend(CeedInit_CudaInterp(data->d_interp_1d, P_1d, Q_1d, &data->c_B)); 62 void *interp_args[] = {(void *)&num_elem, &data->c_B, &d_u, &d_v}; 63 64 if (dim == 1) { 65 CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 66 1)); // avoid >512 total threads 67 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 68 CeedInt shared_mem = elems_per_block * thread_1d * sizeof(CeedScalar); 69 70 if (t_mode == CEED_TRANSPOSE) { 71 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->InterpTransposeAdd : data->InterpTranspose, grid, thread_1d, 1, 72 elems_per_block, shared_mem, interp_args)); 73 } else { 74 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, 1, elems_per_block, shared_mem, interp_args)); 75 } 76 } else if (dim == 2) { 77 const CeedInt opt_elems[7] = {0, 32, 8, 6, 4, 2, 8}; 78 // elems_per_block must be at least 1 79 CeedInt elems_per_block = CeedIntMax(thread_1d < 7 ? opt_elems[thread_1d] / num_comp : 1, 1); 80 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 81 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 82 83 if (t_mode == CEED_TRANSPOSE) { 84 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->InterpTransposeAdd : data->InterpTranspose, grid, thread_1d, thread_1d, 85 elems_per_block, shared_mem, interp_args)); 86 } else { 87 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 88 } 89 } else if (dim == 3) { 90 CeedInt elems_per_block = 1; 91 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 92 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 93 94 if (t_mode == CEED_TRANSPOSE) { 95 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->InterpTransposeAdd : data->InterpTranspose, grid, thread_1d, thread_1d, 96 elems_per_block, shared_mem, interp_args)); 97 } else { 98 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 99 } 100 } 101 } break; 102 case CEED_EVAL_GRAD: { 103 CeedInt P_1d, Q_1d; 104 105 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 106 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 107 CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 108 109 if (data->d_collo_grad_1d) { 110 CeedCallBackend(CeedInit_CudaCollocatedGrad(data->d_interp_1d, data->d_collo_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); 111 } else { 112 CeedCallBackend(CeedInit_CudaGrad(data->d_interp_1d, data->d_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); 113 } 114 void *grad_args[] = {(void *)&num_elem, &data->c_B, &data->c_G, &d_u, &d_v}; 115 if (dim == 1) { 116 CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 117 1)); // avoid >512 total threads 118 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 119 CeedInt shared_mem = elems_per_block * thread_1d * sizeof(CeedScalar); 120 121 if (t_mode == CEED_TRANSPOSE) { 122 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->GradTransposeAdd : data->GradTranspose, grid, thread_1d, 1, 123 elems_per_block, shared_mem, grad_args)); 124 } else { 125 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, 1, elems_per_block, shared_mem, grad_args)); 126 } 127 } else if (dim == 2) { 128 const CeedInt opt_elems[7] = {0, 32, 8, 6, 4, 2, 8}; 129 // elems_per_block must be at least 1 130 CeedInt elems_per_block = CeedIntMax(thread_1d < 7 ? opt_elems[thread_1d] / num_comp : 1, 1); 131 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 132 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 133 134 if (t_mode == CEED_TRANSPOSE) { 135 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->GradTransposeAdd : data->GradTranspose, grid, thread_1d, thread_1d, 136 elems_per_block, shared_mem, grad_args)); 137 } else { 138 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 139 } 140 } else if (dim == 3) { 141 CeedInt elems_per_block = 1; 142 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 143 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 144 145 if (t_mode == CEED_TRANSPOSE) { 146 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->GradTransposeAdd : data->GradTranspose, grid, thread_1d, thread_1d, 147 elems_per_block, shared_mem, grad_args)); 148 } else { 149 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 150 } 151 } 152 } break; 153 case CEED_EVAL_WEIGHT: { 154 CeedInt Q_1d; 155 CeedInt block_size = 32; 156 157 CeedCheck(data->d_q_weight_1d, ceed, CEED_ERROR_BACKEND, "%s not supported; q_weights_1d not set", CeedEvalModes[eval_mode]); 158 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 159 void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v}; 160 if (dim == 1) { 161 const CeedInt elems_per_block = block_size / Q_1d; 162 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 163 164 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, elems_per_block, 1, weight_args)); 165 } else if (dim == 2) { 166 const CeedInt opt_elems = block_size / (Q_1d * Q_1d); 167 const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1; 168 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 169 170 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args)); 171 } else if (dim == 3) { 172 const CeedInt opt_elems = block_size / (Q_1d * Q_1d); 173 const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1; 174 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 175 176 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args)); 177 } 178 } break; 179 case CEED_EVAL_NONE: /* handled separately below */ 180 break; 181 // LCOV_EXCL_START 182 case CEED_EVAL_DIV: 183 case CEED_EVAL_CURL: 184 return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); 185 // LCOV_EXCL_STOP 186 } 187 188 // Restore vectors, cover CEED_EVAL_NONE 189 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 190 if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u)); 191 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 192 CeedCallBackend(CeedDestroy(&ceed)); 193 return CEED_ERROR_SUCCESS; 194 } 195 196 static int CeedBasisApplyTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, 197 CeedVector v) { 198 CeedCallBackend(CeedBasisApplyTensorCore_Cuda_shared(basis, false, num_elem, t_mode, eval_mode, u, v)); 199 return CEED_ERROR_SUCCESS; 200 } 201 202 static int CeedBasisApplyAddTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 203 CeedVector u, CeedVector v) { 204 CeedCallBackend(CeedBasisApplyTensorCore_Cuda_shared(basis, true, num_elem, t_mode, eval_mode, u, v)); 205 return CEED_ERROR_SUCCESS; 206 } 207 208 //------------------------------------------------------------------------------ 209 // Basis apply - tensor AtPoints 210 //------------------------------------------------------------------------------ 211 static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, const CeedInt *num_points, 212 CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 213 Ceed ceed; 214 CeedInt Q_1d, dim, max_num_points = num_points[0]; 215 const CeedInt is_transpose = t_mode == CEED_TRANSPOSE; 216 const int max_block_size = 32; 217 const CeedScalar *d_x, *d_u; 218 CeedScalar *d_v; 219 CeedBasis_Cuda_shared *data; 220 221 CeedCallBackend(CeedBasisGetData(basis, &data)); 222 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 223 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 224 225 // Weight handled separately 226 if (eval_mode == CEED_EVAL_WEIGHT) { 227 CeedCallBackend(CeedVectorSetValue(v, 1.0)); 228 return CEED_ERROR_SUCCESS; 229 } 230 231 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 232 233 // Check padded to uniform number of points per elem 234 for (CeedInt i = 1; i < num_elem; i++) max_num_points = CeedIntMax(max_num_points, num_points[i]); 235 { 236 CeedInt num_comp, q_comp; 237 CeedSize len, len_required; 238 239 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 240 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 241 CeedCallBackend(CeedVectorGetLength(is_transpose ? u : v, &len)); 242 len_required = (CeedSize)num_comp * (CeedSize)q_comp * (CeedSize)num_elem * (CeedSize)max_num_points; 243 CeedCheck(len >= len_required, ceed, CEED_ERROR_BACKEND, 244 "Vector at points must be padded to the same number of points in each element for BasisApplyAtPoints on GPU backends." 245 " Found %" CeedSize_FMT ", Required %" CeedSize_FMT, 246 len, len_required); 247 } 248 249 // Move num_points array to device 250 if (is_transpose) { 251 const CeedInt num_bytes = num_elem * sizeof(CeedInt); 252 253 if (num_elem != data->num_elem_at_points) { 254 data->num_elem_at_points = num_elem; 255 256 if (data->d_points_per_elem) CeedCallCuda(ceed, cudaFree(data->d_points_per_elem)); 257 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_points_per_elem, num_bytes)); 258 CeedCallBackend(CeedFree(&data->h_points_per_elem)); 259 CeedCallBackend(CeedCalloc(num_elem, &data->h_points_per_elem)); 260 } 261 if (memcmp(data->h_points_per_elem, num_points, num_bytes)) { 262 memcpy(data->h_points_per_elem, num_points, num_bytes); 263 CeedCallCuda(ceed, cudaMemcpy(data->d_points_per_elem, num_points, num_bytes, cudaMemcpyHostToDevice)); 264 } 265 } 266 267 // Build kernels if needed 268 if (data->num_points != max_num_points) { 269 CeedInt P_1d; 270 271 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 272 data->num_points = max_num_points; 273 274 // -- Create interp matrix to Chebyshev coefficients 275 if (!data->d_chebyshev_interp_1d) { 276 CeedSize interp_bytes; 277 CeedScalar *chebyshev_interp_1d; 278 279 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 280 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 281 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 282 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes)); 283 CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 284 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 285 } 286 287 // -- Compile kernels 288 const char basis_kernel_source[] = "// AtPoints basis source\n#include <ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h>\n"; 289 CeedInt num_comp; 290 291 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 292 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 293 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->moduleAtPoints, 9, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN", 294 Q_1d * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim - 1), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, 295 "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS", 296 max_num_points, "POINTS_BUFF_LEN", CeedIntPow(Q_1d, dim - 1))); 297 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints)); 298 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpTransposeAtPoints", &data->InterpTransposeAtPoints)); 299 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "GradAtPoints", &data->GradAtPoints)); 300 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "GradTransposeAtPoints", &data->GradTransposeAtPoints)); 301 } 302 303 // Get read/write access to u, v 304 CeedCallBackend(CeedVectorGetArrayRead(x_ref, CEED_MEM_DEVICE, &d_x)); 305 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 306 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 307 if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); 308 else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 309 310 // Clear v for transpose operation 311 if (is_transpose && !apply_add) { 312 CeedInt num_comp, q_comp, num_nodes; 313 CeedSize length; 314 315 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 316 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 317 CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); 318 length = 319 (CeedSize)num_elem * (CeedSize)num_comp * (t_mode == CEED_TRANSPOSE ? (CeedSize)num_nodes : ((CeedSize)max_num_points * (CeedSize)q_comp)); 320 CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar))); 321 } 322 323 // Basis action 324 switch (eval_mode) { 325 case CEED_EVAL_INTERP: { 326 void *interp_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v}; 327 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 328 329 CeedCallBackend( 330 CeedRunKernel_Cuda(ceed, is_transpose ? data->InterpTransposeAtPoints : data->InterpAtPoints, num_elem, block_size, interp_args)); 331 } break; 332 case CEED_EVAL_GRAD: { 333 void *grad_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v}; 334 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 335 336 CeedCallBackend(CeedRunKernel_Cuda(ceed, is_transpose ? data->GradTransposeAtPoints : data->GradAtPoints, num_elem, block_size, grad_args)); 337 } break; 338 case CEED_EVAL_WEIGHT: 339 case CEED_EVAL_NONE: /* handled separately below */ 340 break; 341 // LCOV_EXCL_START 342 case CEED_EVAL_DIV: 343 case CEED_EVAL_CURL: 344 return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); 345 // LCOV_EXCL_STOP 346 } 347 348 // Restore vectors, cover CEED_EVAL_NONE 349 CeedCallBackend(CeedVectorRestoreArrayRead(x_ref, &d_x)); 350 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 351 if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u)); 352 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 353 CeedCallBackend(CeedDestroy(&ceed)); 354 return CEED_ERROR_SUCCESS; 355 } 356 357 static int CeedBasisApplyAtPoints_Cuda_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 358 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 359 CeedCallBackend(CeedBasisApplyAtPointsCore_Cuda_shared(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 360 return CEED_ERROR_SUCCESS; 361 } 362 363 static int CeedBasisApplyAddAtPoints_Cuda_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 364 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 365 CeedCallBackend(CeedBasisApplyAtPointsCore_Cuda_shared(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 366 return CEED_ERROR_SUCCESS; 367 } 368 369 //------------------------------------------------------------------------------ 370 // Destroy basis 371 //------------------------------------------------------------------------------ 372 static int CeedBasisDestroy_Cuda_shared(CeedBasis basis) { 373 Ceed ceed; 374 CeedBasis_Cuda_shared *data; 375 376 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 377 CeedCallBackend(CeedBasisGetData(basis, &data)); 378 CeedCallCuda(ceed, cuModuleUnload(data->module)); 379 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 380 if (data->d_q_weight_1d) CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d)); 381 CeedCallBackend(CeedFree(&data->h_points_per_elem)); 382 if (data->d_points_per_elem) CeedCallCuda(ceed, cudaFree(data->d_points_per_elem)); 383 CeedCallCuda(ceed, cudaFree(data->d_interp_1d)); 384 CeedCallCuda(ceed, cudaFree(data->d_grad_1d)); 385 CeedCallCuda(ceed, cudaFree(data->d_collo_grad_1d)); 386 CeedCallCuda(ceed, cudaFree(data->d_chebyshev_interp_1d)); 387 CeedCallBackend(CeedFree(&data)); 388 CeedCallBackend(CeedDestroy(&ceed)); 389 return CEED_ERROR_SUCCESS; 390 } 391 392 //------------------------------------------------------------------------------ 393 // Create tensor basis 394 //------------------------------------------------------------------------------ 395 int CeedBasisCreateTensorH1_Cuda_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 396 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 397 Ceed ceed; 398 CeedInt num_comp; 399 const CeedInt q_bytes = Q_1d * sizeof(CeedScalar); 400 const CeedInt interp_bytes = q_bytes * P_1d; 401 CeedBasis_Cuda_shared *data; 402 403 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 404 CeedCallBackend(CeedCalloc(1, &data)); 405 406 // Copy basis data to GPU 407 if (q_weight_1d) { 408 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); 409 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice)); 410 } 411 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); 412 CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 413 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes)); 414 CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice)); 415 416 // Compute collocated gradient and copy to GPU 417 data->d_collo_grad_1d = NULL; 418 bool has_collocated_grad = dim == 3 && Q_1d >= P_1d; 419 420 if (has_collocated_grad) { 421 CeedScalar *collo_grad_1d; 422 423 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d)); 424 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d)); 425 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_collo_grad_1d, q_bytes * Q_1d)); 426 CeedCallCuda(ceed, cudaMemcpy(data->d_collo_grad_1d, collo_grad_1d, q_bytes * Q_1d, cudaMemcpyHostToDevice)); 427 CeedCallBackend(CeedFree(&collo_grad_1d)); 428 } 429 430 // Compile basis kernels 431 const char basis_kernel_source[] = "// Tensor basis source\n#include <ceed/jit-source/cuda/cuda-shared-basis-tensor.h>\n"; 432 433 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 434 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "T_1D", 435 CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), 436 "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_HAS_COLLOCATED_GRAD", has_collocated_grad)); 437 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 438 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose)); 439 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTransposeAdd", &data->InterpTransposeAdd)); 440 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 441 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTranspose", &data->GradTranspose)); 442 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTransposeAdd", &data->GradTransposeAdd)); 443 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 444 445 CeedCallBackend(CeedBasisSetData(basis, data)); 446 447 // Register backend functions 448 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Cuda_shared)); 449 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAddTensor_Cuda_shared)); 450 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAtPoints", CeedBasisApplyAtPoints_Cuda_shared)); 451 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAddAtPoints", CeedBasisApplyAddAtPoints_Cuda_shared)); 452 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda_shared)); 453 CeedCallBackend(CeedDestroy(&ceed)); 454 return CEED_ERROR_SUCCESS; 455 } 456 457 //------------------------------------------------------------------------------ 458