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 14 #include "../cuda/ceed-cuda-common.h" 15 #include "../cuda/ceed-cuda-compile.h" 16 #include "ceed-cuda-ref.h" 17 18 //------------------------------------------------------------------------------ 19 // Basis apply - tensor 20 //------------------------------------------------------------------------------ 21 int CeedBasisApply_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 22 Ceed ceed; 23 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 24 Ceed_Cuda *ceed_Cuda; 25 CeedCallBackend(CeedGetData(ceed, &ceed_Cuda)); 26 CeedBasis_Cuda *data; 27 CeedCallBackend(CeedBasisGetData(basis, &data)); 28 const CeedInt transpose = t_mode == CEED_TRANSPOSE; 29 const int max_block_size = 32; 30 31 // Read vectors 32 const CeedScalar *d_u; 33 CeedScalar *d_v; 34 if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 35 else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); 36 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 37 38 // Clear v for transpose operation 39 if (t_mode == CEED_TRANSPOSE) { 40 CeedSize length; 41 CeedCallBackend(CeedVectorGetLength(v, &length)); 42 CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar))); 43 } 44 CeedInt Q_1d, dim; 45 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 46 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 47 48 // Basis action 49 switch (eval_mode) { 50 case CEED_EVAL_INTERP: { 51 void *interp_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_interp_1d, &d_u, &d_v}; 52 const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size); 53 54 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Interp, num_elem, block_size, interp_args)); 55 } break; 56 case CEED_EVAL_GRAD: { 57 void *grad_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_interp_1d, &data->d_grad_1d, &d_u, &d_v}; 58 const CeedInt block_size = max_block_size; 59 60 CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Grad, num_elem, block_size, grad_args)); 61 } break; 62 case CEED_EVAL_WEIGHT: { 63 void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v}; 64 const int block_size_x = Q_1d; 65 const int block_size_y = dim >= 2 ? Q_1d : 1; 66 67 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, num_elem, block_size_x, block_size_y, 1, weight_args)); 68 } break; 69 // LCOV_EXCL_START 70 // Evaluate the divergence to/from the quadrature points 71 case CEED_EVAL_DIV: 72 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_DIV not supported"); 73 // Evaluate the curl to/from the quadrature points 74 case CEED_EVAL_CURL: 75 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported"); 76 // Take no action, BasisApply should not have been called 77 case CEED_EVAL_NONE: 78 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 79 // LCOV_EXCL_STOP 80 } 81 82 // Restore vectors 83 if (eval_mode != CEED_EVAL_WEIGHT) { 84 CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 85 } 86 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 87 return CEED_ERROR_SUCCESS; 88 } 89 90 //------------------------------------------------------------------------------ 91 // Basis apply - non-tensor 92 //------------------------------------------------------------------------------ 93 int CeedBasisApplyNonTensor_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, 94 CeedVector v) { 95 Ceed ceed; 96 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 97 Ceed_Cuda *ceed_Cuda; 98 CeedCallBackend(CeedGetData(ceed, &ceed_Cuda)); 99 CeedBasisNonTensor_Cuda *data; 100 CeedCallBackend(CeedBasisGetData(basis, &data)); 101 CeedInt num_nodes, num_qpts; 102 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 103 CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); 104 const CeedInt transpose = t_mode == CEED_TRANSPOSE; 105 int elems_per_block = 1; 106 int grid = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0); 107 108 // Read vectors 109 const CeedScalar *d_u; 110 CeedScalar *d_v; 111 if (eval_mode != CEED_EVAL_WEIGHT) { 112 CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 113 } 114 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 115 116 // Clear v for transpose operation 117 if (t_mode == CEED_TRANSPOSE) { 118 CeedSize length; 119 CeedCallBackend(CeedVectorGetLength(v, &length)); 120 CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar))); 121 } 122 123 // Apply basis operation 124 switch (eval_mode) { 125 case CEED_EVAL_INTERP: { 126 void *interp_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_interp, &d_u, &d_v}; 127 const int block_size_x = transpose ? num_nodes : num_qpts; 128 129 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Interp, grid, block_size_x, 1, elems_per_block, interp_args)); 130 } break; 131 case CEED_EVAL_GRAD: { 132 void *grad_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_grad, &d_u, &d_v}; 133 const int block_size_x = transpose ? num_nodes : num_qpts; 134 135 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Grad, grid, block_size_x, 1, elems_per_block, grad_args)); 136 } break; 137 case CEED_EVAL_WEIGHT: { 138 void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight, &d_v}; 139 140 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid, num_qpts, 1, elems_per_block, weight_args)); 141 } break; 142 // LCOV_EXCL_START 143 // Evaluate the divergence to/from the quadrature points 144 case CEED_EVAL_DIV: 145 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_DIV not supported"); 146 // Evaluate the curl to/from the quadrature points 147 case CEED_EVAL_CURL: 148 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported"); 149 // Take no action, BasisApply should not have been called 150 case CEED_EVAL_NONE: 151 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context"); 152 // LCOV_EXCL_STOP 153 } 154 155 // Restore vectors 156 if (eval_mode != CEED_EVAL_WEIGHT) { 157 CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 158 } 159 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 160 return CEED_ERROR_SUCCESS; 161 } 162 163 //------------------------------------------------------------------------------ 164 // Destroy tensor basis 165 //------------------------------------------------------------------------------ 166 static int CeedBasisDestroy_Cuda(CeedBasis basis) { 167 Ceed ceed; 168 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 169 170 CeedBasis_Cuda *data; 171 CeedCallBackend(CeedBasisGetData(basis, &data)); 172 173 CeedCallCuda(ceed, cuModuleUnload(data->module)); 174 175 CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d)); 176 CeedCallCuda(ceed, cudaFree(data->d_interp_1d)); 177 CeedCallCuda(ceed, cudaFree(data->d_grad_1d)); 178 CeedCallBackend(CeedFree(&data)); 179 180 return CEED_ERROR_SUCCESS; 181 } 182 183 //------------------------------------------------------------------------------ 184 // Destroy non-tensor basis 185 //------------------------------------------------------------------------------ 186 static int CeedBasisDestroyNonTensor_Cuda(CeedBasis basis) { 187 Ceed ceed; 188 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 189 190 CeedBasisNonTensor_Cuda *data; 191 CeedCallBackend(CeedBasisGetData(basis, &data)); 192 193 CeedCallCuda(ceed, cuModuleUnload(data->module)); 194 195 CeedCallCuda(ceed, cudaFree(data->d_q_weight)); 196 CeedCallCuda(ceed, cudaFree(data->d_interp)); 197 CeedCallCuda(ceed, cudaFree(data->d_grad)); 198 CeedCallBackend(CeedFree(&data)); 199 200 return CEED_ERROR_SUCCESS; 201 } 202 203 //------------------------------------------------------------------------------ 204 // Create tensor 205 //------------------------------------------------------------------------------ 206 int CeedBasisCreateTensorH1_Cuda(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d, 207 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) { 208 Ceed ceed; 209 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 210 CeedBasis_Cuda *data; 211 CeedCallBackend(CeedCalloc(1, &data)); 212 213 // Copy data to GPU 214 const CeedInt q_bytes = Q_1d * sizeof(CeedScalar); 215 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); 216 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice)); 217 218 const CeedInt interp_bytes = q_bytes * P_1d; 219 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); 220 CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 221 222 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes)); 223 CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice)); 224 225 // Compile basis kernels 226 CeedInt num_comp; 227 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 228 char *basis_kernel_path, *basis_kernel_source; 229 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-tensor.h", &basis_kernel_path)); 230 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n"); 231 CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source)); 232 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n"); 233 CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 7, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN", 234 num_comp * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, 235 "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim))); 236 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 237 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 238 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 239 CeedCallBackend(CeedFree(&basis_kernel_path)); 240 CeedCallBackend(CeedFree(&basis_kernel_source)); 241 242 CeedCallBackend(CeedBasisSetData(basis, data)); 243 244 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Cuda)); 245 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda)); 246 return CEED_ERROR_SUCCESS; 247 } 248 249 //------------------------------------------------------------------------------ 250 // Create non-tensor 251 //------------------------------------------------------------------------------ 252 int CeedBasisCreateH1_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad, 253 const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 254 Ceed ceed; 255 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 256 CeedBasisNonTensor_Cuda *data; 257 CeedCallBackend(CeedCalloc(1, &data)); 258 259 // Copy basis data to GPU 260 const CeedInt q_bytes = num_qpts * sizeof(CeedScalar); 261 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes)); 262 CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice)); 263 264 const CeedInt interp_bytes = q_bytes * num_nodes; 265 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes)); 266 CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice)); 267 268 const CeedInt grad_bytes = q_bytes * num_nodes * dim; 269 CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad, grad_bytes)); 270 CeedCallCuda(ceed, cudaMemcpy(data->d_grad, grad, grad_bytes, cudaMemcpyHostToDevice)); 271 272 // Compile basis kernels 273 CeedInt num_comp; 274 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 275 char *basis_kernel_path, *basis_kernel_source; 276 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path)); 277 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n"); 278 CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source)); 279 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n"); 280 CeedCallCuda(ceed, CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 4, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_DIM", dim, 281 "BASIS_NUM_COMP", num_comp)); 282 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); 283 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); 284 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); 285 CeedCallBackend(CeedFree(&basis_kernel_path)); 286 CeedCallBackend(CeedFree(&basis_kernel_source)); 287 288 CeedCallBackend(CeedBasisSetData(basis, data)); 289 290 // Register backend functions 291 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda)); 292 CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda)); 293 return CEED_ERROR_SUCCESS; 294 } 295 296 //------------------------------------------------------------------------------ 297