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 return CEED_ERROR_SUCCESS; 193 } 194 195 static int CeedBasisApplyTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, 196 CeedVector v) { 197 CeedCallBackend(CeedBasisApplyTensorCore_Cuda_shared(basis, false, num_elem, t_mode, eval_mode, u, v)); 198 return CEED_ERROR_SUCCESS; 199 } 200 201 static int CeedBasisApplyAddTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 202 CeedVector u, CeedVector v) { 203 CeedCallBackend(CeedBasisApplyTensorCore_Cuda_shared(basis, true, num_elem, t_mode, eval_mode, u, v)); 204 return CEED_ERROR_SUCCESS; 205 } 206 207 //------------------------------------------------------------------------------ 208 // Basis apply - tensor AtPoints 209 //------------------------------------------------------------------------------ 210 static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, const CeedInt *num_points, 211 CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 212 Ceed ceed; 213 CeedInt Q_1d, dim, max_num_points = num_points[0]; 214 const CeedInt is_transpose = t_mode == CEED_TRANSPOSE; 215 const int max_block_size = 32; 216 const CeedScalar *d_x, *d_u; 217 CeedScalar *d_v; 218 CeedBasis_Cuda_shared *data; 219 220 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 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 // Check padded to uniform number of points per elem 232 for (CeedInt i = 1; i < num_elem; i++) max_num_points = CeedIntMax(max_num_points, num_points[i]); 233 { 234 CeedInt num_comp, q_comp; 235 CeedSize len, len_required; 236 237 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 238 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 239 CeedCallBackend(CeedVectorGetLength(is_transpose ? u : v, &len)); 240 len_required = (CeedSize)num_comp * (CeedSize)q_comp * (CeedSize)num_elem * (CeedSize)max_num_points; 241 CeedCheck(len >= len_required, ceed, CEED_ERROR_BACKEND, 242 "Vector at points must be padded to the same number of points in each element for BasisApplyAtPoints on GPU backends." 243 " Found %" CeedSize_FMT ", Required %" CeedSize_FMT, 244 len, len_required); 245 } 246 247 // Move num_points array to device 248 if (is_transpose) { 249 const CeedInt num_bytes = num_elem * sizeof(CeedInt); 250 251 if (num_elem != data->num_elem_at_points) { 252 data->num_elem_at_points = num_elem; 253 254 if (data->d_points_per_elem) CeedCallCuda(ceed, cudaFree(data->d_points_per_elem)); 255 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_points_per_elem, num_bytes)); 256 CeedCallBackend(CeedFree(&data->h_points_per_elem)); 257 CeedCallBackend(CeedCalloc(num_elem, &data->h_points_per_elem)); 258 } 259 if (memcmp(data->h_points_per_elem, num_points, num_bytes)) { 260 memcpy(data->h_points_per_elem, num_points, num_bytes); 261 CeedCallCuda(ceed, cudaMemcpy(data->d_points_per_elem, num_points, num_bytes, cudaMemcpyHostToDevice)); 262 } 263 } 264 265 // Build kernels if needed 266 if (data->num_points != max_num_points) { 267 CeedInt P_1d; 268 269 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 270 data->num_points = max_num_points; 271 272 // -- Create interp matrix to Chebyshev coefficients 273 if (!data->d_chebyshev_interp_1d) { 274 CeedSize interp_bytes; 275 CeedScalar *chebyshev_interp_1d; 276 277 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 278 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 279 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 280 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes)); 281 CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 282 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 283 } 284 285 // -- Compile kernels 286 char *basis_kernel_source; 287 const char *basis_kernel_path; 288 CeedInt num_comp; 289 290 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 291 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 292 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h", &basis_kernel_path)); 293 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n"); 294 CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source)); 295 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n"); 296 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->moduleAtPoints, 9, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN", 297 Q_1d * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim - 1), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, 298 "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS", 299 max_num_points, "POINTS_BUFF_LEN", CeedIntPow(Q_1d, dim - 1))); 300 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints)); 301 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "GradAtPoints", &data->GradAtPoints)); 302 CeedCallBackend(CeedFree(&basis_kernel_path)); 303 CeedCallBackend(CeedFree(&basis_kernel_source)); 304 } 305 306 // Get read/write access to u, v 307 CeedCallBackend(CeedVectorGetArrayRead(x_ref, CEED_MEM_DEVICE, &d_x)); 308 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 309 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 310 if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); 311 else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 312 313 // Clear v for transpose operation 314 if (is_transpose && !apply_add) { 315 CeedInt num_comp, q_comp, num_nodes; 316 CeedSize length; 317 318 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 319 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 320 CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); 321 length = 322 (CeedSize)num_elem * (CeedSize)num_comp * (t_mode == CEED_TRANSPOSE ? (CeedSize)num_nodes : ((CeedSize)max_num_points * (CeedSize)q_comp)); 323 CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar))); 324 } 325 326 // Basis action 327 switch (eval_mode) { 328 case CEED_EVAL_INTERP: { 329 void *interp_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v}; 330 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 331 332 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->InterpAtPoints, num_elem, block_size, interp_args)); 333 } break; 334 case CEED_EVAL_GRAD: { 335 void *grad_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v}; 336 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 337 338 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->GradAtPoints, num_elem, block_size, grad_args)); 339 } break; 340 case CEED_EVAL_WEIGHT: 341 case CEED_EVAL_NONE: /* handled separately below */ 342 break; 343 // LCOV_EXCL_START 344 case CEED_EVAL_DIV: 345 case CEED_EVAL_CURL: 346 return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); 347 // LCOV_EXCL_STOP 348 } 349 350 // Restore vectors, cover CEED_EVAL_NONE 351 CeedCallBackend(CeedVectorRestoreArrayRead(x_ref, &d_x)); 352 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 353 if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u)); 354 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 355 return CEED_ERROR_SUCCESS; 356 } 357 358 static int CeedBasisApplyAtPoints_Cuda_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 359 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 360 CeedCallBackend(CeedBasisApplyAtPointsCore_Cuda_shared(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 361 return CEED_ERROR_SUCCESS; 362 } 363 364 static int CeedBasisApplyAddAtPoints_Cuda_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 365 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 366 CeedCallBackend(CeedBasisApplyAtPointsCore_Cuda_shared(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 367 return CEED_ERROR_SUCCESS; 368 } 369 370 //------------------------------------------------------------------------------ 371 // Destroy basis 372 //------------------------------------------------------------------------------ 373 static int CeedBasisDestroy_Cuda_shared(CeedBasis basis) { 374 Ceed ceed; 375 CeedBasis_Cuda_shared *data; 376 377 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 378 CeedCallBackend(CeedBasisGetData(basis, &data)); 379 CeedCallCuda(ceed, cuModuleUnload(data->module)); 380 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 381 if (data->d_q_weight_1d) CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d)); 382 CeedCallBackend(CeedFree(&data->h_points_per_elem)); 383 if (data->d_points_per_elem) CeedCallCuda(ceed, cudaFree(data->d_points_per_elem)); 384 CeedCallCuda(ceed, cudaFree(data->d_interp_1d)); 385 CeedCallCuda(ceed, cudaFree(data->d_grad_1d)); 386 CeedCallCuda(ceed, cudaFree(data->d_collo_grad_1d)); 387 CeedCallCuda(ceed, cudaFree(data->d_chebyshev_interp_1d)); 388 CeedCallBackend(CeedFree(&data)); 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 char *basis_kernel_source; 399 const char *basis_kernel_path; 400 CeedInt num_comp; 401 const CeedInt q_bytes = Q_1d * sizeof(CeedScalar); 402 const CeedInt interp_bytes = q_bytes * P_1d; 403 CeedBasis_Cuda_shared *data; 404 405 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 406 CeedCallBackend(CeedCalloc(1, &data)); 407 408 // Copy basis data to GPU 409 if (q_weight_1d) { 410 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); 411 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice)); 412 } 413 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); 414 CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 415 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes)); 416 CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice)); 417 418 // Compute collocated gradient and copy to GPU 419 data->d_collo_grad_1d = NULL; 420 bool has_collocated_grad = dim == 3 && Q_1d >= P_1d; 421 422 if (has_collocated_grad) { 423 CeedScalar *collo_grad_1d; 424 425 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d)); 426 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d)); 427 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_collo_grad_1d, q_bytes * Q_1d)); 428 CeedCallCuda(ceed, cudaMemcpy(data->d_collo_grad_1d, collo_grad_1d, q_bytes * Q_1d, cudaMemcpyHostToDevice)); 429 CeedCallBackend(CeedFree(&collo_grad_1d)); 430 } 431 432 // Compile basis kernels 433 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 434 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-shared-basis-tensor.h", &basis_kernel_path)); 435 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n"); 436 CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source)); 437 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete -----\n"); 438 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "T_1D", 439 CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), 440 "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_HAS_COLLOCATED_GRAD", has_collocated_grad)); 441 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 442 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose)); 443 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTransposeAdd", &data->InterpTransposeAdd)); 444 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 445 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTranspose", &data->GradTranspose)); 446 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTransposeAdd", &data->GradTransposeAdd)); 447 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 448 CeedCallBackend(CeedFree(&basis_kernel_path)); 449 CeedCallBackend(CeedFree(&basis_kernel_source)); 450 451 CeedCallBackend(CeedBasisSetData(basis, data)); 452 453 // Register backend functions 454 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Cuda_shared)); 455 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAddTensor_Cuda_shared)); 456 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAtPoints", CeedBasisApplyAtPoints_Cuda_shared)); 457 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAddAtPoints", CeedBasisApplyAddAtPoints_Cuda_shared)); 458 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda_shared)); 459 return CEED_ERROR_SUCCESS; 460 } 461 462 //------------------------------------------------------------------------------ 463