1 // Copyright (c) 2017-2022, 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 16 #include "../cuda/ceed-cuda-common.h" 17 #include "../cuda/ceed-cuda-compile.h" 18 #include "ceed-cuda-shared.h" 19 20 //------------------------------------------------------------------------------ 21 // Device initalization 22 //------------------------------------------------------------------------------ 23 int CeedInit_CudaInterp(CeedScalar *d_B, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B); 24 int CeedInit_CudaGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); 25 int CeedInit_CudaCollocatedGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); 26 27 //------------------------------------------------------------------------------ 28 // Apply basis 29 //------------------------------------------------------------------------------ 30 int CeedBasisApplyTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, 31 CeedVector v) { 32 Ceed ceed; 33 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 34 Ceed_Cuda *ceed_Cuda; 35 CeedCallBackend(CeedGetData(ceed, &ceed_Cuda)); 36 CeedBasis_Cuda_shared *data; 37 CeedCallBackend(CeedBasisGetData(basis, &data)); 38 CeedInt dim, num_comp; 39 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 40 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 41 42 // Read vectors 43 const CeedScalar *d_u; 44 CeedScalar *d_v; 45 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 46 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 47 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 48 49 // Apply basis operation 50 switch (eval_mode) { 51 case CEED_EVAL_INTERP: { 52 CeedInt P_1d, Q_1d; 53 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 54 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 55 CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 56 CeedCallBackend(CeedInit_CudaInterp(data->d_interp_1d, P_1d, Q_1d, &data->c_B)); 57 void *interp_args[] = {(void *)&num_elem, &data->c_B, &d_u, &d_v}; 58 if (dim == 1) { 59 CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 60 1)); // avoid >512 total threads 61 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 62 CeedInt shared_mem = elems_per_block * thread_1d * sizeof(CeedScalar); 63 64 if (t_mode == CEED_TRANSPOSE) { 65 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->InterpTranspose, grid, thread_1d, 1, elems_per_block, shared_mem, interp_args)); 66 } else { 67 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, 1, elems_per_block, shared_mem, interp_args)); 68 } 69 } else if (dim == 2) { 70 const CeedInt opt_elems[7] = {0, 32, 8, 6, 4, 2, 8}; 71 // elems_per_block must be at least 1 72 CeedInt elems_per_block = CeedIntMax(thread_1d < 7 ? opt_elems[thread_1d] / num_comp : 1, 1); 73 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 74 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 75 76 if (t_mode == CEED_TRANSPOSE) { 77 CeedCallBackend( 78 CeedRunKernelDimShared_Cuda(ceed, data->InterpTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 79 } else { 80 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 81 } 82 } else if (dim == 3) { 83 CeedInt elems_per_block = 1; 84 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 85 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 86 87 if (t_mode == CEED_TRANSPOSE) { 88 CeedCallBackend( 89 CeedRunKernelDimShared_Cuda(ceed, data->InterpTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 90 } else { 91 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args)); 92 } 93 } 94 } break; 95 case CEED_EVAL_GRAD: { 96 CeedInt P_1d, Q_1d; 97 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 98 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 99 CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 100 if (data->d_collo_grad_1d) { 101 CeedCallBackend(CeedInit_CudaCollocatedGrad(data->d_interp_1d, data->d_collo_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); 102 } else { 103 CeedCallBackend(CeedInit_CudaGrad(data->d_interp_1d, data->d_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); 104 } 105 void *grad_args[] = {(void *)&num_elem, &data->c_B, &data->c_G, &d_u, &d_v}; 106 if (dim == 1) { 107 CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 108 1)); // avoid >512 total threads 109 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 110 CeedInt shared_mem = elems_per_block * thread_1d * sizeof(CeedScalar); 111 112 if (t_mode == CEED_TRANSPOSE) { 113 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->GradTranspose, grid, thread_1d, 1, elems_per_block, shared_mem, grad_args)); 114 } else { 115 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, 1, elems_per_block, shared_mem, grad_args)); 116 } 117 } else if (dim == 2) { 118 const CeedInt opt_elems[7] = {0, 32, 8, 6, 4, 2, 8}; 119 // elems_per_block must be at least 1 120 CeedInt elems_per_block = CeedIntMax(thread_1d < 7 ? opt_elems[thread_1d] / num_comp : 1, 1); 121 CeedInt grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 122 CeedInt shared_mem = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar); 123 124 if (t_mode == CEED_TRANSPOSE) { 125 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->GradTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 126 } else { 127 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 128 } 129 } else if (dim == 3) { 130 CeedInt elems_per_block = 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, data->GradTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 136 } else { 137 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args)); 138 } 139 } 140 } break; 141 case CEED_EVAL_WEIGHT: { 142 CeedInt Q_1d; 143 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 144 void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v}; 145 if (dim == 1) { 146 const CeedInt elems_per_block = 32 / Q_1d; 147 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 148 149 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, elems_per_block, 1, weight_args)); 150 } else if (dim == 2) { 151 const CeedInt opt_elems = 32 / (Q_1d * Q_1d); 152 const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1; 153 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 154 155 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args)); 156 } else if (dim == 3) { 157 const CeedInt opt_elems = 32 / (Q_1d * Q_1d); 158 const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1; 159 const CeedInt grid_size = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 160 161 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args)); 162 } 163 } break; 164 // LCOV_EXCL_START 165 // Evaluate the divergence to/from the quadrature points 166 case CEED_EVAL_DIV: 167 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_DIV not supported"); 168 // Evaluate the curl to/from the quadrature points 169 case CEED_EVAL_CURL: 170 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported"); 171 // Take no action, BasisApply should not have been called 172 case CEED_EVAL_NONE: 173 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 174 // LCOV_EXCL_STOP 175 } 176 177 // Restore vectors 178 if (eval_mode != CEED_EVAL_WEIGHT) { 179 CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 180 } 181 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 182 return CEED_ERROR_SUCCESS; 183 } 184 185 //------------------------------------------------------------------------------ 186 // Destroy basis 187 //------------------------------------------------------------------------------ 188 static int CeedBasisDestroy_Cuda_shared(CeedBasis basis) { 189 Ceed ceed; 190 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 191 192 CeedBasis_Cuda_shared *data; 193 CeedCallBackend(CeedBasisGetData(basis, &data)); 194 195 CeedCallCuda(ceed, cuModuleUnload(data->module)); 196 197 CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d)); 198 CeedCallCuda(ceed, cudaFree(data->d_interp_1d)); 199 CeedCallCuda(ceed, cudaFree(data->d_grad_1d)); 200 CeedCallCuda(ceed, cudaFree(data->d_collo_grad_1d)); 201 202 CeedCallBackend(CeedFree(&data)); 203 204 return CEED_ERROR_SUCCESS; 205 } 206 207 //------------------------------------------------------------------------------ 208 // Create tensor basis 209 //------------------------------------------------------------------------------ 210 int CeedBasisCreateTensorH1_Cuda_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 211 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 212 Ceed ceed; 213 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 214 CeedBasis_Cuda_shared *data; 215 CeedCallBackend(CeedCalloc(1, &data)); 216 217 // Copy basis data to GPU 218 const CeedInt q_bytes = Q_1d * sizeof(CeedScalar); 219 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); 220 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice)); 221 222 const CeedInt interp_bytes = q_bytes * P_1d; 223 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); 224 CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 225 226 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes)); 227 CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice)); 228 229 // Compute collocated gradient and copy to GPU 230 data->d_collo_grad_1d = NULL; 231 bool has_collocated_grad = dim == 3 && Q_1d >= P_1d; 232 if (has_collocated_grad) { 233 CeedScalar *collo_grad_1d; 234 CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d)); 235 CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d)); 236 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_collo_grad_1d, q_bytes * Q_1d)); 237 CeedCallCuda(ceed, cudaMemcpy(data->d_collo_grad_1d, collo_grad_1d, q_bytes * Q_1d, cudaMemcpyHostToDevice)); 238 CeedCallBackend(CeedFree(&collo_grad_1d)); 239 } 240 241 // Compile basis kernels 242 CeedInt num_comp; 243 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 244 char *basis_kernel_path, *basis_kernel_source; 245 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-shared-basis-tensor.h", &basis_kernel_path)); 246 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n"); 247 CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source)); 248 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete -----\n"); 249 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "T_1D", 250 CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), 251 "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_HAS_COLLOCATED_GRAD", has_collocated_grad)); 252 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 253 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose)); 254 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 255 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTranspose", &data->GradTranspose)); 256 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 257 CeedCallBackend(CeedFree(&basis_kernel_path)); 258 CeedCallBackend(CeedFree(&basis_kernel_source)); 259 260 CeedCallBackend(CeedBasisSetData(basis, data)); 261 262 // Register backend functions 263 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Cuda_shared)); 264 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda_shared)); 265 return CEED_ERROR_SUCCESS; 266 } 267 268 //------------------------------------------------------------------------------ 269