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