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