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