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 const char basis_kernel_source[] = "// AtPoints basis source\n#include <ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h>\n"; 287 CeedInt num_comp; 288 289 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 290 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 291 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->moduleAtPoints, 9, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN", 292 Q_1d * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim - 1), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, 293 "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS", 294 max_num_points, "POINTS_BUFF_LEN", CeedIntPow(Q_1d, dim - 1))); 295 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints)); 296 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "GradAtPoints", &data->GradAtPoints)); 297 } 298 299 // Get read/write access to u, v 300 CeedCallBackend(CeedVectorGetArrayRead(x_ref, CEED_MEM_DEVICE, &d_x)); 301 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 302 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 303 if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); 304 else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 305 306 // Clear v for transpose operation 307 if (is_transpose && !apply_add) { 308 CeedInt num_comp, q_comp, num_nodes; 309 CeedSize length; 310 311 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 312 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 313 CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); 314 length = 315 (CeedSize)num_elem * (CeedSize)num_comp * (t_mode == CEED_TRANSPOSE ? (CeedSize)num_nodes : ((CeedSize)max_num_points * (CeedSize)q_comp)); 316 CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar))); 317 } 318 319 // Basis action 320 switch (eval_mode) { 321 case CEED_EVAL_INTERP: { 322 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}; 323 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 324 325 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->InterpAtPoints, num_elem, block_size, interp_args)); 326 } break; 327 case CEED_EVAL_GRAD: { 328 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}; 329 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 330 331 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->GradAtPoints, num_elem, block_size, grad_args)); 332 } break; 333 case CEED_EVAL_WEIGHT: 334 case CEED_EVAL_NONE: /* handled separately below */ 335 break; 336 // LCOV_EXCL_START 337 case CEED_EVAL_DIV: 338 case CEED_EVAL_CURL: 339 return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); 340 // LCOV_EXCL_STOP 341 } 342 343 // Restore vectors, cover CEED_EVAL_NONE 344 CeedCallBackend(CeedVectorRestoreArrayRead(x_ref, &d_x)); 345 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 346 if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u)); 347 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 348 return CEED_ERROR_SUCCESS; 349 } 350 351 static int CeedBasisApplyAtPoints_Cuda_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 352 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 353 CeedCallBackend(CeedBasisApplyAtPointsCore_Cuda_shared(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 354 return CEED_ERROR_SUCCESS; 355 } 356 357 static int CeedBasisApplyAddAtPoints_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, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 360 return CEED_ERROR_SUCCESS; 361 } 362 363 //------------------------------------------------------------------------------ 364 // Destroy basis 365 //------------------------------------------------------------------------------ 366 static int CeedBasisDestroy_Cuda_shared(CeedBasis basis) { 367 Ceed ceed; 368 CeedBasis_Cuda_shared *data; 369 370 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 371 CeedCallBackend(CeedBasisGetData(basis, &data)); 372 CeedCallCuda(ceed, cuModuleUnload(data->module)); 373 if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints)); 374 if (data->d_q_weight_1d) CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d)); 375 CeedCallBackend(CeedFree(&data->h_points_per_elem)); 376 if (data->d_points_per_elem) CeedCallCuda(ceed, cudaFree(data->d_points_per_elem)); 377 CeedCallCuda(ceed, cudaFree(data->d_interp_1d)); 378 CeedCallCuda(ceed, cudaFree(data->d_grad_1d)); 379 CeedCallCuda(ceed, cudaFree(data->d_collo_grad_1d)); 380 CeedCallCuda(ceed, cudaFree(data->d_chebyshev_interp_1d)); 381 CeedCallBackend(CeedFree(&data)); 382 return CEED_ERROR_SUCCESS; 383 } 384 385 //------------------------------------------------------------------------------ 386 // Create tensor basis 387 //------------------------------------------------------------------------------ 388 int CeedBasisCreateTensorH1_Cuda_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 389 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 390 Ceed ceed; 391 CeedInt num_comp; 392 const CeedInt q_bytes = Q_1d * sizeof(CeedScalar); 393 const CeedInt interp_bytes = q_bytes * P_1d; 394 CeedBasis_Cuda_shared *data; 395 396 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 397 CeedCallBackend(CeedCalloc(1, &data)); 398 399 // Copy basis data to GPU 400 if (q_weight_1d) { 401 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); 402 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice)); 403 } 404 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); 405 CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 406 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes)); 407 CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice)); 408 409 // Compute collocated gradient and copy to GPU 410 data->d_collo_grad_1d = NULL; 411 bool has_collocated_grad = dim == 3 && Q_1d >= P_1d; 412 413 if (has_collocated_grad) { 414 CeedScalar *collo_grad_1d; 415 416 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d)); 417 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d)); 418 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_collo_grad_1d, q_bytes * Q_1d)); 419 CeedCallCuda(ceed, cudaMemcpy(data->d_collo_grad_1d, collo_grad_1d, q_bytes * Q_1d, cudaMemcpyHostToDevice)); 420 CeedCallBackend(CeedFree(&collo_grad_1d)); 421 } 422 423 // Compile basis kernels 424 const char basis_kernel_source[] = "// Tensor basis source\n#include <ceed/jit-source/cuda/cuda-shared-basis-tensor.h>\n"; 425 426 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 427 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "T_1D", 428 CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), 429 "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_HAS_COLLOCATED_GRAD", has_collocated_grad)); 430 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 431 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose)); 432 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTransposeAdd", &data->InterpTransposeAdd)); 433 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 434 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTranspose", &data->GradTranspose)); 435 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTransposeAdd", &data->GradTransposeAdd)); 436 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 437 438 CeedCallBackend(CeedBasisSetData(basis, data)); 439 440 // Register backend functions 441 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Cuda_shared)); 442 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAddTensor_Cuda_shared)); 443 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAtPoints", CeedBasisApplyAtPoints_Cuda_shared)); 444 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAddAtPoints", CeedBasisApplyAddAtPoints_Cuda_shared)); 445 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda_shared)); 446 return CEED_ERROR_SUCCESS; 447 } 448 449 //------------------------------------------------------------------------------ 450