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 #define CEED_DEBUG_COLOR 12 9 10 #include <ceed.h> 11 #include <ceed/backend.h> 12 #include <ceed/jit-tools.h> 13 #include <cuda_runtime.h> 14 15 #include <iostream> 16 #include <sstream> 17 #include <string> 18 19 #include "../cuda-ref/ceed-cuda-ref.h" 20 #include "../cuda-shared/ceed-cuda-shared.h" 21 #include "../cuda/ceed-cuda-common.h" 22 #include "../cuda/ceed-cuda-compile.h" 23 #include "ceed-cuda-gen.h" 24 25 //------------------------------------------------------------------------------ 26 // Determine type of operator 27 //------------------------------------------------------------------------------ 28 static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fields, CeedOperatorField *op_input_fields, 29 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, 30 CeedQFunctionField *qf_output_fields, CeedInt *max_P_1d, CeedInt *Q_1d, CeedInt *dim, bool *is_tensor, 31 bool *use_3d_slices) { 32 // Find dim, P_1d, Q_1d 33 *max_P_1d = 0; 34 *Q_1d = 0; 35 *dim = 0; 36 *is_tensor = true; 37 for (CeedInt i = 0; i < num_input_fields; i++) { 38 CeedBasis basis; 39 40 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 41 if (basis != CEED_BASIS_NONE) { 42 bool is_field_tensor; 43 CeedInt field_P_1d = 0, field_Q_1d = 0, field_dim = 0; 44 45 // Collect dim, P_1d, and Q_1d 46 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 47 CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 48 *is_tensor = *is_tensor && is_field_tensor; 49 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 50 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 51 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 52 CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 53 *dim = field_dim; 54 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 55 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 56 *Q_1d = field_Q_1d; 57 } 58 } 59 for (CeedInt i = 0; i < num_output_fields; i++) { 60 CeedBasis basis; 61 62 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 63 if (basis != CEED_BASIS_NONE) { 64 bool is_field_tensor; 65 CeedInt field_P_1d = 0, field_Q_1d = 0, field_dim = 0; 66 67 // Collect dim, P_1d, and Q_1d 68 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 69 CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 70 *is_tensor = *is_tensor && is_field_tensor; 71 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 72 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 73 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 74 CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 75 *dim = field_dim; 76 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 77 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 78 *Q_1d = field_Q_1d; 79 } 80 } 81 82 // Only use 3D collocated gradient parallelization strategy when gradient is computed 83 *use_3d_slices = false; 84 if (*dim == 3) { 85 bool was_grad_found = false; 86 87 for (CeedInt i = 0; i < num_input_fields; i++) { 88 CeedEvalMode eval_mode; 89 90 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 91 if (eval_mode == CEED_EVAL_GRAD) { 92 CeedBasis_Cuda_shared *basis_data; 93 CeedBasis basis; 94 95 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 96 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 97 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 98 was_grad_found = true; 99 } 100 } 101 for (CeedInt i = 0; i < num_output_fields; i++) { 102 CeedEvalMode eval_mode; 103 104 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 105 if (eval_mode == CEED_EVAL_GRAD) { 106 CeedBasis_Cuda_shared *basis_data; 107 CeedBasis basis; 108 109 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 110 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 111 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 112 was_grad_found = true; 113 } 114 } 115 } 116 return CEED_ERROR_SUCCESS; 117 } 118 119 //------------------------------------------------------------------------------ 120 // Setup fields 121 //------------------------------------------------------------------------------ 122 static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedOperatorField op_field, 123 CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool use_3d_slices) { 124 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 125 std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; 126 std::string option_name = (is_input ? "inputs" : "outputs"); 127 CeedEvalMode eval_mode = CEED_EVAL_NONE; 128 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 129 CeedElemRestriction elem_rstr; 130 CeedBasis_Cuda_shared *basis_data; 131 CeedBasis basis; 132 133 code << " // -- " << (is_input ? "Input" : "Output") << " field " << i << "\n"; 134 135 // Get field data 136 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 137 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 138 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 139 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 140 } 141 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 142 if (basis != CEED_BASIS_NONE) { 143 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 144 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 145 } 146 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 147 148 // Set field constants 149 if (eval_mode != CEED_EVAL_WEIGHT) { 150 code << " const CeedInt " << P_name << " = " << (basis == CEED_BASIS_NONE ? Q_1d : P_1d) << ";\n"; 151 code << " const CeedInt num_comp" << var_suffix << " = " << num_comp << ";\n"; 152 } 153 154 // Load basis data 155 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 156 switch (eval_mode) { 157 case CEED_EVAL_NONE: 158 break; 159 case CEED_EVAL_INTERP: 160 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 161 else data->B.outputs[i] = basis_data->d_interp_1d; 162 code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; 163 code << " loadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 164 break; 165 case CEED_EVAL_GRAD: 166 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 167 else data->B.outputs[i] = basis_data->d_interp_1d; 168 code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; 169 code << " loadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 170 if (use_3d_slices) { 171 if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; 172 else data->G.outputs[i] = basis_data->d_collo_grad_1d; 173 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; 174 code << " loadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 175 } else { 176 bool has_collo_grad = basis_data->d_collo_grad_1d; 177 178 if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 179 else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 180 if (has_collo_grad) { 181 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; 182 code << " loadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 183 } else { 184 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * P_1d << "];\n"; 185 code << " loadMatrix<" << P_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 186 } 187 } 188 break; 189 case CEED_EVAL_WEIGHT: 190 break; // No action 191 // LCOV_EXCL_START 192 case CEED_EVAL_DIV: 193 break; // TODO: Not implemented 194 case CEED_EVAL_CURL: 195 break; // TODO: Not implemented 196 // LCOV_EXCL_STOP 197 } 198 return CEED_ERROR_SUCCESS; 199 } 200 201 //------------------------------------------------------------------------------ 202 // Restriction 203 //------------------------------------------------------------------------------ 204 static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedInt dim, 205 CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, 206 bool use_3d_slices) { 207 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 208 std::string P_name = "P_1d" + var_suffix; 209 CeedEvalMode eval_mode = CEED_EVAL_NONE; 210 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 211 CeedSize l_size; 212 CeedElemRestriction_Cuda *rstr_data; 213 CeedElemRestriction elem_rstr; 214 CeedBasis basis; 215 216 // Get field data 217 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 218 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 219 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 220 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 221 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 222 } 223 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 224 if (basis != CEED_BASIS_NONE) { 225 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 226 } 227 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 228 229 // Restriction 230 if (is_input) { 231 // Input 232 if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices)) { 233 bool is_strided; 234 235 code << " CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; 236 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 237 if (!is_strided) { 238 CeedInt comp_stride; 239 240 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 241 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 242 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 243 code << " // CompStride: " << comp_stride << "\n"; 244 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 245 code << " readDofsOffset" << dim << "d<num_comp" << var_suffix << ", " << comp_stride << ", " << P_name << ">(data, l_size" << var_suffix 246 << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; 247 } else { 248 bool has_backend_strides; 249 CeedInt num_elem; 250 251 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 252 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 253 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 254 255 if (!has_backend_strides) { 256 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 257 } 258 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 259 code << " readDofsStrided" << dim << "d<num_comp" << var_suffix << ", " << P_name << "," << strides[0] << "," << strides[1] << "," 260 << strides[2] << ">(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; 261 } 262 } 263 } else { 264 // Output 265 bool is_strided; 266 267 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 268 if (!is_strided) { 269 CeedInt comp_stride; 270 271 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 272 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 273 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 274 code << " // CompStride: " << comp_stride << "\n"; 275 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; 276 code << " writeDofsOffset" << dim << "d<num_comp" << var_suffix << ", " << comp_stride << ", " << P_name << ">(data, l_size" << var_suffix 277 << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; 278 } else { 279 bool has_backend_strides; 280 CeedInt num_elem; 281 282 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 283 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 284 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 285 286 if (!has_backend_strides) { 287 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 288 } 289 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 290 code << " writeDofsStrided" << dim << "d<num_comp" << var_suffix << ", " << P_name << "," << strides[0] << "," << strides[1] << "," 291 << strides[2] << ">(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; 292 } 293 } 294 return CEED_ERROR_SUCCESS; 295 } 296 297 //------------------------------------------------------------------------------ 298 // Basis 299 //------------------------------------------------------------------------------ 300 static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedInt dim, 301 CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, 302 bool use_3d_slices) { 303 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 304 std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; 305 CeedEvalMode eval_mode = CEED_EVAL_NONE; 306 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 307 CeedElemRestriction elem_rstr; 308 CeedBasis basis; 309 310 // Get field data 311 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 312 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 313 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 314 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 315 } 316 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 317 if (basis != CEED_BASIS_NONE) { 318 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 319 } 320 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 321 322 // Basis 323 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 324 if (is_input) { 325 switch (eval_mode) { 326 case CEED_EVAL_NONE: 327 if (!use_3d_slices) { 328 code << " CeedScalar *r_q" << var_suffix << " = r_e" << var_suffix << ";\n"; 329 } 330 break; 331 case CEED_EVAL_INTERP: 332 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 333 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 334 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 335 break; 336 case CEED_EVAL_GRAD: 337 if (use_3d_slices) { 338 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 339 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 340 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 341 } else { 342 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; 343 code << " Grad" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d<num_comp" << var_suffix 344 << ", P_1d" << var_suffix << ", " << Q_name << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" 345 << var_suffix << ");\n"; 346 } 347 break; 348 case CEED_EVAL_WEIGHT: { 349 CeedBasis_Cuda_shared *basis_data; 350 351 code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; 352 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 353 data->W = basis_data->d_q_weight_1d; 354 code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; 355 break; 356 } 357 // LCOV_EXCL_START 358 case CEED_EVAL_DIV: 359 break; // TODO: Not implemented 360 case CEED_EVAL_CURL: 361 break; // TODO: Not implemented 362 // LCOV_EXCL_STOP 363 } 364 } else { 365 switch (eval_mode) { 366 case CEED_EVAL_NONE: 367 code << " CeedScalar *r_e" << var_suffix << " = r_q" << var_suffix << ";\n"; 368 break; // No action 369 case CEED_EVAL_INTERP: 370 code << " CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; 371 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 372 << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 373 break; 374 case CEED_EVAL_GRAD: 375 code << " CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; 376 if (use_3d_slices) { 377 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 378 << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 379 } else { 380 code << " GradTranspose" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d<num_comp" 381 << var_suffix << ", " << P_name << "," << Q_name << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix 382 << ", r_e" << var_suffix << ");\n"; 383 } 384 break; 385 // LCOV_EXCL_START 386 case CEED_EVAL_WEIGHT: 387 break; // Should not occur 388 case CEED_EVAL_DIV: 389 break; // TODO: Not implemented 390 case CEED_EVAL_CURL: 391 break; // TODO: Not implemented 392 // LCOV_EXCL_STOP 393 } 394 } 395 return CEED_ERROR_SUCCESS; 396 } 397 398 //------------------------------------------------------------------------------ 399 // QFunction 400 //------------------------------------------------------------------------------ 401 static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt dim, CeedInt num_input_fields, 402 CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, 403 CeedInt num_output_fields, CeedOperatorField *op_output_fields, 404 CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, 405 bool use_3d_slices) { 406 std::string Q_name = "Q_1d"; 407 CeedEvalMode eval_mode = CEED_EVAL_NONE; 408 CeedElemRestriction elem_rstr; 409 410 // Setup output arays 411 code << "\n // -- Output field setup\n"; 412 for (CeedInt i = 0; i < num_output_fields; i++) { 413 std::string var_suffix = "_out_" + std::to_string(i); 414 415 code << " // ---- Output field " << i << "\n"; 416 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 417 if (eval_mode == CEED_EVAL_NONE || eval_mode == CEED_EVAL_INTERP) { 418 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 419 } 420 if (eval_mode == CEED_EVAL_GRAD) { 421 if (use_3d_slices) { 422 // Accumulator for gradient slices 423 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 424 code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) {\n"; 425 code << " r_q" << var_suffix << "[i] = 0.0;\n"; 426 code << " }\n"; 427 } else { 428 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; 429 } 430 } 431 } 432 433 // We treat quadrature points per slice in 3d to save registers 434 if (use_3d_slices) { 435 code << "\n // Note: Using planes of 3D elements\n"; 436 code << "#pragma unroll\n"; 437 code << " for (CeedInt q = 0; q < " << Q_name << "; q++) {\n"; 438 code << " // -- Input fields\n"; 439 for (CeedInt i = 0; i < num_input_fields; i++) { 440 std::string var_suffix = "_in_" + std::to_string(i); 441 442 code << " // ---- Input field " << i << "\n"; 443 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 444 // Basis action 445 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 446 switch (eval_mode) { 447 case CEED_EVAL_NONE: 448 bool is_strided; 449 450 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 451 452 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 453 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 454 if (is_strided) { 455 bool has_backend_strides; 456 CeedInt num_elem, elem_size; 457 458 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 459 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 460 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 461 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 462 463 if (!has_backend_strides) { 464 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 465 } 466 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 467 code << " readSliceQuadsStrided3d<num_comp" << var_suffix << ", " << Q_name << "," << strides[0] << "," << strides[1] << "," 468 << strides[2] << ">(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n"; 469 } else { 470 CeedSize l_size = 0; 471 CeedInt comp_stride; 472 CeedElemRestriction_Cuda *rstr_data; 473 474 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 475 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 476 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 477 code << " // CompStride: " << comp_stride << "\n"; 478 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 479 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 480 code << " readSliceQuadsOffset3d<num_comp" << var_suffix << ", " << comp_stride << ", " << Q_name << ">(data, l_size" << var_suffix 481 << ", elem, q, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; 482 } 483 break; 484 case CEED_EVAL_INTERP: 485 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 486 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) {\n"; 487 code << " r_s" << var_suffix << "[j] = r_q" << var_suffix << "[q + j*" << Q_name << "];\n"; 488 code << " }\n"; 489 break; 490 case CEED_EVAL_GRAD: 491 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 492 code << " gradCollo3d<num_comp" << var_suffix << ", " << Q_name << ">(data, q, r_q" << var_suffix << ", s_G" << var_suffix << ", r_s" 493 << var_suffix << ");\n"; 494 break; 495 case CEED_EVAL_WEIGHT: 496 code << " CeedScalar r_s" << var_suffix << "[1];\n"; 497 code << " r_s" << var_suffix << "[0] = r_q" << var_suffix << "[q];\n"; 498 break; // No action 499 // LCOV_EXCL_START 500 case CEED_EVAL_DIV: 501 break; // TODO: Not implemented 502 case CEED_EVAL_CURL: 503 break; // TODO: Not implemented 504 // LCOV_EXCL_STOP 505 } 506 } 507 code << "\n // -- Output fields\n"; 508 for (CeedInt i = 0; i < num_output_fields; i++) { 509 std::string var_suffix = "_out_" + std::to_string(i); 510 511 code << " // ---- Output field " << i << "\n"; 512 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 513 // Basis action 514 switch (eval_mode) { 515 case CEED_EVAL_NONE: 516 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 517 break; // No action 518 case CEED_EVAL_INTERP: 519 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 520 break; 521 case CEED_EVAL_GRAD: 522 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 523 break; 524 // LCOV_EXCL_START 525 case CEED_EVAL_WEIGHT: 526 break; // Should not occur 527 case CEED_EVAL_DIV: 528 break; // TODO: Not implemented 529 case CEED_EVAL_CURL: 530 break; // TODO: Not implemented 531 // LCOV_EXCL_STOP 532 } 533 } 534 } else { 535 code << "\n // Note: Using full elements\n"; 536 code << " {\n"; 537 code << " // -- Input fields\n"; 538 for (CeedInt i = 0; i < num_input_fields; i++) { 539 code << " // ---- Input field " << i << "\n"; 540 code << " CeedScalar *r_s_in_" << i << " = r_q_in_" << i << ";\n"; 541 } 542 code << " // -- Output fields\n"; 543 for (CeedInt i = 0; i < num_output_fields; i++) { 544 code << " // ---- Output field " << i << "\n"; 545 code << " CeedScalar *r_s_out_" << i << " = r_q_out_" << i << ";\n"; 546 } 547 } 548 549 // Input and output buffers 550 code << "\n // -- QFunction inputs and outputs\n"; 551 code << " // ---- Inputs\n"; 552 code << " CeedScalar *inputs[" << CeedIntMax(num_input_fields, 1) << "];\n"; 553 for (CeedInt i = 0; i < num_input_fields; i++) { 554 code << " // ------ Input field " << i << "\n"; 555 code << " inputs[" << i << "] = r_s_in_" << i << ";\n"; 556 } 557 code << " // ---- Outputs\n"; 558 code << " CeedScalar *outputs[" << CeedIntMax(num_output_fields, 1) << "];\n"; 559 for (CeedInt i = 0; i < num_output_fields; i++) { 560 code << " // ------ Output field " << i << "\n"; 561 code << " outputs[" << i << "] = r_s_out_" << i << ";\n"; 562 } 563 564 // Apply QFunction 565 code << "\n // -- Apply QFunction\n"; 566 code << " " << qfunction_name << "(ctx, "; 567 if (dim != 3 || use_3d_slices) { 568 code << "1"; 569 } else { 570 code << "Q_1d"; 571 } 572 code << ", inputs, outputs);\n"; 573 574 // Copy or apply transpose grad, if needed 575 if (use_3d_slices) { 576 code << " // -- Output fields\n"; 577 for (CeedInt i = 0; i < num_output_fields; i++) { 578 std::string var_suffix = "_out_" + std::to_string(i); 579 std::string P_name = "P_1d" + var_suffix; 580 581 code << " // ---- Output field " << i << "\n"; 582 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 583 // Basis action 584 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 585 switch (eval_mode) { 586 case CEED_EVAL_NONE: 587 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 588 code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 589 code << " }\n"; 590 break; // No action 591 case CEED_EVAL_INTERP: 592 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 593 code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 594 code << " }\n"; 595 break; 596 case CEED_EVAL_GRAD: 597 code << " gradColloTranspose3d<num_comp" << var_suffix << ", " << Q_name << ">(data, q, r_s" << var_suffix << ", s_G" << var_suffix 598 << ", r_q" << var_suffix << ");\n"; 599 break; 600 // LCOV_EXCL_START 601 case CEED_EVAL_WEIGHT: 602 break; // Should not occur 603 case CEED_EVAL_DIV: 604 break; // TODO: Not implemented 605 case CEED_EVAL_CURL: 606 break; // TODO: Not implemented 607 // LCOV_EXCL_STOP 608 } 609 } 610 } 611 code << " }\n"; 612 return CEED_ERROR_SUCCESS; 613 } 614 615 //------------------------------------------------------------------------------ 616 // Build single operator kernel 617 //------------------------------------------------------------------------------ 618 extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { 619 bool is_tensor = true, use_3d_slices = false; 620 Ceed ceed; 621 CeedInt Q_1d, num_input_fields, num_output_fields, dim = 1; 622 CeedQFunctionField *qf_input_fields, *qf_output_fields; 623 CeedQFunction_Cuda_gen *qf_data; 624 CeedQFunction qf; 625 CeedOperatorField *op_input_fields, *op_output_fields; 626 CeedOperator_Cuda_gen *data; 627 std::ostringstream code; 628 629 { 630 bool is_setup_done; 631 632 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 633 if (is_setup_done) return CEED_ERROR_SUCCESS; 634 } 635 636 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 637 CeedCallBackend(CeedOperatorGetData(op, &data)); 638 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 639 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 640 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 641 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 642 643 // Get operator data 644 CeedCallBackend(CeedOperatorBuildKernelData_Cuda_gen(ceed, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, op_output_fields, 645 qf_output_fields, &data->max_P_1d, &Q_1d, &dim, &is_tensor, &use_3d_slices)); 646 if (dim == 0) dim = 1; 647 data->dim = dim; 648 if (Q_1d == 0) { 649 CeedInt Q; 650 651 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 652 Q_1d = Q; 653 } 654 data->Q_1d = Q_1d; 655 656 // Check for restriction only identity operator 657 { 658 bool is_identity_qf; 659 660 CeedCallBackend(CeedQFunctionIsIdentity(qf, &is_identity_qf)); 661 if (is_identity_qf) { 662 CeedEvalMode eval_mode_in, eval_mode_out; 663 664 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[0], &eval_mode_in)); 665 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[0], &eval_mode_out)); 666 CeedCheck(eval_mode_in != CEED_EVAL_NONE || eval_mode_out != CEED_EVAL_NONE, ceed, CEED_ERROR_BACKEND, 667 "Backend does not implement restriction only identity operators"); 668 } 669 } 670 671 // Add atomicAdd function for old NVidia architectures 672 { 673 Ceed_Cuda *ceed_data; 674 struct cudaDeviceProp prop; 675 676 CeedCallBackend(CeedGetData(ceed, &ceed_data)); 677 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id)); 678 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) { 679 char *atomic_add_source; 680 const char *atomic_add_path; 681 682 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-atomic-add-fallback.h", &atomic_add_path)); 683 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Atomic Add Source -----\n"); 684 CeedCallBackend(CeedLoadSourceToBuffer(ceed, atomic_add_path, &atomic_add_source)); 685 code << atomic_add_source; 686 CeedCallBackend(CeedFree(&atomic_add_path)); 687 CeedCallBackend(CeedFree(&atomic_add_source)); 688 } 689 } 690 691 // Load basis source files 692 // TODO: Add non-tensor, AtPoints 693 { 694 char *tensor_basis_kernel_source; 695 const char *tensor_basis_kernel_path; 696 697 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h", &tensor_basis_kernel_path)); 698 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Tensor Basis Kernel Source -----\n"); 699 CeedCallBackend(CeedLoadSourceToBuffer(ceed, tensor_basis_kernel_path, &tensor_basis_kernel_source)); 700 code << tensor_basis_kernel_source; 701 CeedCallBackend(CeedFree(&tensor_basis_kernel_path)); 702 CeedCallBackend(CeedFree(&tensor_basis_kernel_source)); 703 } 704 { 705 char *cuda_gen_template_source; 706 const char *cuda_gen_template_path; 707 708 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-gen-templates.h", &cuda_gen_template_path)); 709 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Cuda-Gen Template Source -----\n"); 710 CeedCallBackend(CeedLoadSourceToBuffer(ceed, cuda_gen_template_path, &cuda_gen_template_source)); 711 code << cuda_gen_template_source; 712 CeedCallBackend(CeedFree(&cuda_gen_template_path)); 713 CeedCallBackend(CeedFree(&cuda_gen_template_source)); 714 } 715 716 // Get QFunction name 717 std::string qfunction_name(qf_data->qfunction_name); 718 std::string operator_name; 719 720 operator_name = "CeedKernelCudaGenOperator_" + qfunction_name; 721 722 // Define CEED_Q_VLA 723 code << "\n#undef CEED_Q_VLA\n"; 724 if (dim != 3 || use_3d_slices) { 725 code << "#define CEED_Q_VLA 1\n\n"; 726 } else { 727 code << "#define CEED_Q_VLA " << Q_1d << "\n\n"; 728 } 729 730 // Add user QFunction source 731 { 732 std::string qfunction_source(qf_data->qfunction_source); 733 734 code << qfunction_source; 735 } 736 737 // Setup 738 code << "\n// -----------------------------------------------------------------------------\n"; 739 code << "// Operator Kernel\n"; 740 code << "// \n"; 741 code << "// d_[in,out]_i: CeedVector device array\n"; 742 code << "// r_[in,out]_e_i: Element vector register\n"; 743 code << "// r_[in,out]_q_i: Quadrature space vector register\n"; 744 code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; 745 code << "// \n"; 746 code << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; 747 code << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; 748 code << "// -----------------------------------------------------------------------------\n"; 749 code << "extern \"C\" __global__ void " << operator_name 750 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W) {\n"; 751 752 // Scratch buffers 753 for (CeedInt i = 0; i < num_input_fields; i++) { 754 CeedEvalMode eval_mode; 755 756 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 757 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 758 code << " const CeedScalar *d_in_" << i << " = fields.inputs[" << i << "];\n"; 759 } 760 } 761 for (CeedInt i = 0; i < num_output_fields; i++) { 762 code << " CeedScalar *d_out_" << i << " = fields.outputs[" << i << "];\n"; 763 } 764 765 code << " const CeedInt dim = " << dim << ";\n"; 766 code << " const CeedInt Q_1d = " << Q_1d << ";\n"; 767 768 // Shared data 769 code << " extern __shared__ CeedScalar slice[];\n"; 770 code << " SharedData_Cuda data;\n"; 771 code << " data.t_id_x = threadIdx.x;\n"; 772 code << " data.t_id_y = threadIdx.y;\n"; 773 code << " data.t_id_z = threadIdx.z;\n"; 774 code << " data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 775 code << " data.slice = slice + data.t_id_z*T_1D" << (dim > 1 ? "*T_1D" : "") << ";\n"; 776 777 // Initialize constants, and matrices B and G 778 code << "\n // Input field constants and basis data\n"; 779 for (CeedInt i = 0; i < num_input_fields; i++) { 780 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, use_3d_slices)); 781 } 782 code << "\n // Output field constants and basis data\n"; 783 for (CeedInt i = 0; i < num_output_fields; i++) { 784 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); 785 } 786 787 // Loop over all elements 788 code << "\n // Element loop\n"; 789 code << " __syncthreads();\n"; 790 code << " for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; 791 792 // -- Input restriction and basis 793 code << " // -- Input field restrictions and basis actions\n"; 794 for (CeedInt i = 0; i < num_input_fields; i++) { 795 code << " // ---- Input field " << i << "\n"; 796 797 // ---- Restriction 798 CeedCallBackend( 799 CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, i, dim, op_input_fields[i], qf_input_fields[i], Q_1d, true, use_3d_slices)); 800 801 // ---- Basis action 802 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, i, dim, op_input_fields[i], qf_input_fields[i], Q_1d, true, use_3d_slices)); 803 } 804 805 // -- Q function 806 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, dim, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, 807 op_output_fields, qf_output_fields, qfunction_name, Q_1d, use_3d_slices)); 808 809 // -- Output basis and restriction 810 code << "\n // -- Output field basis action and restrictions\n"; 811 for (CeedInt i = 0; i < num_output_fields; i++) { 812 code << " // ---- Output field " << i << "\n"; 813 814 // ---- Basis action 815 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); 816 817 // ---- Restriction 818 CeedCallBackend( 819 CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); 820 } 821 822 // Close loop and function 823 code << " }\n"; 824 code << "}\n"; 825 code << "// -----------------------------------------------------------------------------\n\n"; 826 827 // View kernel for debugging 828 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "Generated Operator Kernels:\n"); 829 CeedDebug(ceed, code.str().c_str()); 830 831 CeedCallBackend(CeedCompile_Cuda(ceed, code.str().c_str(), &data->module, 1, "T_1D", CeedIntMax(Q_1d, data->max_P_1d))); 832 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, operator_name.c_str(), &data->op)); 833 834 CeedCallBackend(CeedOperatorSetSetupDone(op)); 835 return CEED_ERROR_SUCCESS; 836 } 837 838 //------------------------------------------------------------------------------ 839