1 // Copyright (c) 2017-2026, 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/gen-tools.h> 13 #include <ceed/jit-tools.h> 14 #include <cuda_runtime.h> 15 16 #include <iostream> 17 #include <sstream> 18 #include <string> 19 20 #include "../cuda-ref/ceed-cuda-ref.h" 21 #include "../cuda-shared/ceed-cuda-shared.h" 22 #include "../cuda/ceed-cuda-common.h" 23 #include "../cuda/ceed-cuda-compile.h" 24 #include "ceed-cuda-gen.h" 25 26 struct FieldReuse_Cuda { 27 CeedInt index; 28 bool is_input; 29 CeedEvalMode eval_mode; 30 }; 31 32 //------------------------------------------------------------------------------ 33 // Determine type of operator 34 //------------------------------------------------------------------------------ 35 static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fields, CeedOperatorField *op_input_fields, 36 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, 37 CeedQFunctionField *qf_output_fields, CeedInt *max_P, CeedInt *max_P_1d, CeedInt *Q, CeedInt *Q_1d, 38 CeedInt *max_dim, bool *is_all_tensor, bool *use_3d_slices) { 39 // Check if all are tensor 40 *is_all_tensor = true; 41 for (CeedInt i = 0; i < num_input_fields; i++) { 42 CeedBasis basis; 43 44 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 45 if (basis != CEED_BASIS_NONE) { 46 bool is_field_tensor; 47 48 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 49 *is_all_tensor = *is_all_tensor && is_field_tensor; 50 } 51 CeedCallBackend(CeedBasisDestroy(&basis)); 52 } 53 for (CeedInt i = 0; i < num_output_fields; i++) { 54 CeedBasis basis; 55 56 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 57 if (basis != CEED_BASIS_NONE) { 58 bool is_field_tensor; 59 60 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 61 *is_all_tensor = *is_all_tensor && is_field_tensor; 62 } 63 CeedCallBackend(CeedBasisDestroy(&basis)); 64 } 65 66 // Find max_P, max_P_1d, Q, and Q_1d 67 bool is_all_3d = true; 68 69 *max_P = 0; 70 *max_P_1d = 0; 71 *Q = 0; 72 *Q_1d = 0; 73 for (CeedInt i = 0; i < num_input_fields; i++) { 74 CeedBasis basis; 75 76 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 77 if (basis != CEED_BASIS_NONE) { 78 bool is_field_tensor; 79 CeedInt field_dim = 0, field_P = 0, field_P_1d = 0, field_Q = 0, field_Q_1d = 0; 80 81 // Check if 3D 82 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 83 is_all_3d = is_all_3d && (field_dim == 3); 84 *max_dim = CeedIntMax(*max_dim, field_dim); 85 86 // Collect P, P_1d, Q, and Q_1d 87 CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P)); 88 *max_P = CeedIntMax(*max_P, field_P); 89 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 90 if (is_field_tensor) { 91 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 92 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 93 } 94 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q)); 95 CeedCheck(*Q == 0 || field_Q == *Q, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 96 *Q = field_Q; 97 if (is_field_tensor) { 98 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 99 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 100 *Q_1d = field_Q_1d; 101 } 102 } 103 CeedCallBackend(CeedBasisDestroy(&basis)); 104 } 105 for (CeedInt i = 0; i < num_output_fields; i++) { 106 CeedBasis basis; 107 108 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 109 if (basis != CEED_BASIS_NONE) { 110 bool is_field_tensor; 111 CeedInt field_dim = 0, field_P = 0, field_P_1d = 0, field_Q = 0, field_Q_1d = 0; 112 113 // Check if 3D 114 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 115 is_all_3d = is_all_3d && (field_dim == 3); 116 *max_dim = CeedIntMax(*max_dim, field_dim); 117 118 // Collect P, P_1d, Q, and Q_1d 119 CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P)); 120 *max_P = CeedIntMax(*max_P, field_P); 121 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 122 if (is_field_tensor) { 123 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 124 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 125 } 126 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q)); 127 CeedCheck(*Q == 0 || field_Q == *Q, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 128 *Q = field_Q; 129 if (is_field_tensor) { 130 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 131 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 132 *Q_1d = field_Q_1d; 133 } 134 } 135 CeedCallBackend(CeedBasisDestroy(&basis)); 136 } 137 138 // Only use 3D collocated gradient parallelization strategy when gradient is computed 139 *use_3d_slices = false; 140 if (is_all_3d && *is_all_tensor) { 141 bool was_grad_found = false; 142 143 for (CeedInt i = 0; i < num_input_fields; i++) { 144 CeedEvalMode eval_mode; 145 146 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 147 if (eval_mode == CEED_EVAL_GRAD) { 148 CeedBasis_Cuda_shared *basis_data; 149 CeedBasis basis; 150 151 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 152 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 153 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 154 was_grad_found = true; 155 CeedCallBackend(CeedBasisDestroy(&basis)); 156 } 157 } 158 for (CeedInt i = 0; i < num_output_fields; i++) { 159 CeedEvalMode eval_mode; 160 161 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 162 if (eval_mode == CEED_EVAL_GRAD) { 163 CeedBasis_Cuda_shared *basis_data; 164 CeedBasis basis; 165 166 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 167 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 168 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 169 was_grad_found = true; 170 CeedCallBackend(CeedBasisDestroy(&basis)); 171 } 172 } 173 } 174 return CEED_ERROR_SUCCESS; 175 } 176 177 //------------------------------------------------------------------------------ 178 // Setup fields 179 //------------------------------------------------------------------------------ 180 static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i, 181 CeedOperatorField op_field, CeedQFunctionField qf_field, FieldReuse_Cuda field_reuse, 182 CeedInt max_dim, CeedInt Q, CeedInt Q_1d, bool is_input, bool is_all_tensor, bool is_at_points, 183 bool use_3d_slices, bool skip_active_load) { 184 bool is_tensor = true, is_active = true; 185 CeedBasis basis; 186 187 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 188 if (basis != CEED_BASIS_NONE) CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 189 { 190 CeedVector vec; 191 192 CeedCallBackend(CeedOperatorFieldGetVector(op_field, &vec)); 193 is_active = vec == CEED_VECTOR_ACTIVE; 194 CeedCallBackend(CeedVectorDestroy(&vec)); 195 } 196 197 const char *field_name; 198 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 199 std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; 200 std::string option_name = (is_input ? "inputs" : "outputs"); 201 CeedEvalMode eval_mode = CEED_EVAL_NONE; 202 CeedInt elem_size = 0, num_comp = 0, dim = max_dim, P_1d = 0; 203 CeedElemRestriction elem_rstr; 204 CeedBasis_Cuda_shared *basis_data; 205 206 // Field reuse info 207 bool use_previous_field = field_reuse.index != -1; 208 209 CeedCallBackend(CeedOperatorFieldGetName(op_field, &field_name)); 210 code << tab << "// -- " << (is_input ? "Input" : "Output") << " field " << i << ": " << field_name << "\n"; 211 212 // Get field data 213 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 214 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 215 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 216 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 217 } 218 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 219 if (basis != CEED_BASIS_NONE) { 220 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 221 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 222 if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 223 else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); 224 } 225 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 226 227 // Set field constants 228 code << tab << "const CeedInt dim" << var_suffix << " = " << dim << ";\n"; 229 if (is_tensor && !is_all_tensor) { 230 CeedInt P = 0; 231 232 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 233 code << tab << "const CeedInt P" << var_suffix << " = " << (basis == CEED_BASIS_NONE ? Q : P) << ";\n"; 234 } 235 code << tab << "const CeedInt " << P_name << " = " << (basis == CEED_BASIS_NONE ? Q_1d : P_1d) << ";\n"; 236 if (eval_mode != CEED_EVAL_WEIGHT) { 237 code << tab << "const CeedInt num_comp" << var_suffix << " = " << num_comp << ";\n"; 238 } 239 240 // Load basis data 241 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 242 switch (eval_mode) { 243 case CEED_EVAL_NONE: 244 break; 245 case CEED_EVAL_INTERP: 246 if (is_at_points) { 247 // AtPoints 248 if (!basis_data->d_chebyshev_interp_1d) { 249 CeedSize interp_bytes; 250 CeedScalar *chebyshev_interp_1d; 251 252 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 253 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 254 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 255 CeedCallCuda(CeedBasisReturnCeed(basis), cudaMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); 256 CeedCallCuda(CeedBasisReturnCeed(basis), 257 cudaMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 258 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 259 } 260 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; 261 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; 262 } else { 263 // Standard quadrature 264 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 265 else data->B.outputs[i] = basis_data->d_interp_1d; 266 } 267 if (use_previous_field && !skip_active_load) { 268 std::string reuse_var = "s_B" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index)); 269 270 code << tab << "CeedScalar *s_B" << var_suffix << " = " << reuse_var << ";\n"; 271 } else { 272 bool is_collocated = false; 273 274 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated)); 275 if ((is_active && skip_active_load) || (is_collocated && !is_at_points)) { 276 code << tab << "CeedScalar *s_B" << var_suffix << " = NULL;\n"; 277 } else { 278 code << tab << "__shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; 279 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 280 } 281 } 282 break; 283 case CEED_EVAL_GRAD: 284 if (is_at_points) { 285 // AtPoints 286 if (!basis_data->d_chebyshev_interp_1d) { 287 CeedSize interp_bytes; 288 CeedScalar *chebyshev_interp_1d; 289 290 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 291 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 292 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 293 CeedCallCuda(CeedBasisReturnCeed(basis), cudaMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); 294 CeedCallCuda(CeedBasisReturnCeed(basis), 295 cudaMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice)); 296 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 297 } 298 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; 299 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; 300 } else { 301 // Standard quadrature 302 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 303 else data->B.outputs[i] = basis_data->d_interp_1d; 304 } 305 if (is_tensor) { 306 if (use_previous_field && !skip_active_load) { 307 std::string reuse_var = "s_B" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index)); 308 309 code << tab << "CeedScalar *s_B" << var_suffix << " = " << reuse_var << ";\n"; 310 } else { 311 bool is_collocated = false; 312 313 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated)); 314 if ((is_active && skip_active_load) || (is_collocated && !is_at_points)) { 315 code << tab << "CeedScalar *s_B" << var_suffix << " = NULL;\n"; 316 } else { 317 code << tab << "__shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; 318 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 319 } 320 } 321 } 322 if (is_at_points) break; // No G mat for AtPoints 323 if (use_3d_slices) { 324 if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; 325 else data->G.outputs[i] = basis_data->d_collo_grad_1d; 326 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) { 327 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index)); 328 329 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n"; 330 } else if (is_active && skip_active_load) { 331 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n"; 332 } else { 333 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; 334 code << tab << "LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 335 } 336 } else { 337 bool has_collo_grad = basis_data->d_collo_grad_1d; 338 339 if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 340 else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 341 if (has_collo_grad) { 342 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) { 343 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index)); 344 345 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n"; 346 } else if (is_active && skip_active_load) { 347 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n"; 348 } else { 349 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; 350 code << tab << "LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 351 } 352 } else { 353 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) { 354 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index)); 355 356 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n"; 357 } else if (is_active && skip_active_load) { 358 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n"; 359 } else { 360 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << P_name << "*" << Q_name << (is_tensor ? "" : "*dim") 361 << (is_tensor ? "" : var_suffix) << "];\n"; 362 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << (is_tensor ? "" : "*dim") << (is_tensor ? "" : var_suffix) << ">(data, G." 363 << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 364 } 365 } 366 } 367 break; 368 case CEED_EVAL_WEIGHT: 369 break; // No action 370 // LCOV_EXCL_START 371 case CEED_EVAL_DIV: 372 case CEED_EVAL_CURL: 373 break; // TODO: Not implemented 374 // LCOV_EXCL_STOP 375 } 376 CeedCallBackend(CeedBasisDestroy(&basis)); 377 return CEED_ERROR_SUCCESS; 378 } 379 380 //------------------------------------------------------------------------------ 381 // Restriction 382 //------------------------------------------------------------------------------ 383 static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i, 384 CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field, 385 CeedInt max_dim, CeedInt Q_1d, bool is_input, bool is_all_tensor, bool is_at_points, 386 bool use_3d_slices) { 387 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 388 std::string P_name = (is_all_tensor ? "P_1d" : "P") + var_suffix; 389 CeedEvalMode eval_mode = CEED_EVAL_NONE; 390 CeedInt elem_size = 0, num_comp = 0; 391 CeedSize l_size; 392 CeedRestrictionType rstr_type = CEED_RESTRICTION_STANDARD; 393 CeedElemRestriction_Cuda *rstr_data; 394 CeedElemRestriction elem_rstr; 395 396 // Get field data 397 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 398 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 399 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 400 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 401 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 402 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 403 } 404 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 405 406 // Restriction 407 if (is_input) { 408 // Input 409 if (field_input_buffer[i] != i) { 410 std::string buffer_name = "r_e_in_" + std::to_string(field_input_buffer[i]); 411 412 // Restriction was already done for previous input 413 code << tab << "CeedScalar *r_e" << var_suffix << " = " << buffer_name << ";\n"; 414 } else if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices && is_at_points)) { 415 if (eval_mode == CEED_EVAL_NONE && rstr_type != CEED_RESTRICTION_POINTS) { 416 // No basis action, so r_e_in_* in also r_q_in_* and needs to be allocated 417 code << tab << "CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; 418 } else if (rstr_type != CEED_RESTRICTION_POINTS) { 419 // Otherwise we're using the scratch space 420 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 421 } 422 switch (rstr_type) { 423 case CEED_RESTRICTION_STANDARD: { 424 CeedInt comp_stride; 425 426 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 427 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 428 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 429 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 430 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 431 code << tab << "ReadLVecStandard" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", " 432 << P_name << ">(data, l_size" << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix 433 << ");\n"; 434 break; 435 } 436 case CEED_RESTRICTION_STRIDED: { 437 bool has_backend_strides; 438 CeedInt num_elem; 439 440 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 441 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 442 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 443 444 if (!has_backend_strides) { 445 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 446 } 447 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1] 448 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n"; 449 code << tab << "ReadLVecStrided" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", " << P_name << ", strides" 450 << var_suffix << "_0, strides" << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, d" << var_suffix << ", r_e" 451 << var_suffix << ");\n"; 452 break; 453 } 454 case CEED_RESTRICTION_POINTS: { 455 CeedInt comp_stride; 456 457 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 458 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 459 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 460 break; 461 } 462 // LCOV_EXCL_START 463 case CEED_RESTRICTION_ORIENTED: 464 case CEED_RESTRICTION_CURL_ORIENTED: 465 break; // TODO: Not implemented 466 // LCOV_EXCL_STOP 467 } 468 } 469 } else { 470 // Output 471 switch (rstr_type) { 472 case CEED_RESTRICTION_STANDARD: { 473 CeedInt comp_stride; 474 475 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 476 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 477 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 478 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 479 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; 480 code << tab << "WriteLVecStandard" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", " 481 << P_name << ">(data, l_size" << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix 482 << ");\n"; 483 break; 484 } 485 case CEED_RESTRICTION_STRIDED: { 486 bool has_backend_strides; 487 CeedInt num_elem; 488 489 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 490 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 491 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 492 493 if (!has_backend_strides) { 494 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 495 } 496 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1] 497 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n"; 498 code << tab << "WriteLVecStrided" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", " << P_name << ", strides" 499 << var_suffix << "_0, strides" << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, r_e" << var_suffix << ", d" << var_suffix 500 << ");\n"; 501 break; 502 } 503 case CEED_RESTRICTION_POINTS: 504 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; 505 break; 506 // LCOV_EXCL_START 507 case CEED_RESTRICTION_ORIENTED: 508 case CEED_RESTRICTION_CURL_ORIENTED: 509 break; // TODO: Not implemented 510 // LCOV_EXCL_STOP 511 } 512 } 513 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 514 return CEED_ERROR_SUCCESS; 515 } 516 517 //------------------------------------------------------------------------------ 518 // Basis 519 //------------------------------------------------------------------------------ 520 static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i, 521 CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt max_dim, CeedInt Q_1d, 522 bool is_input, bool is_all_tensor, bool is_at_points, bool use_3d_slices) { 523 bool is_tensor = true, is_collocated = true; 524 CeedBasis basis; 525 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 526 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 527 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated)); 528 529 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 530 std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; 531 CeedEvalMode eval_mode = CEED_EVAL_NONE; 532 CeedInt dim = max_dim, elem_size = 0, num_comp = 0, P_1d = 0; 533 CeedElemRestriction elem_rstr; 534 535 // Get field data 536 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 537 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 538 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 539 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 540 } 541 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 542 if (basis != CEED_BASIS_NONE) { 543 CeedCallBackend(CeedBasisGetDimension(basis, &dim)); 544 if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 545 else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); 546 } 547 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 548 549 // Basis 550 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 551 if (is_input) { 552 switch (eval_mode) { 553 case CEED_EVAL_NONE: 554 if (!use_3d_slices && !is_at_points) { 555 code << tab << "CeedScalar *r_q" << var_suffix << " = r_e" << var_suffix << ";\n"; 556 } 557 break; 558 case CEED_EVAL_INTERP: 559 if (is_at_points) { 560 std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; 561 562 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; 563 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix 564 << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; 565 } else { 566 std::string function_name = is_tensor ? ((dim == 1 ? "Interp" : "InterpTensor") + std::string(is_collocated ? "CollocatedNodes" : "") + 567 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened")) 568 : "InterpNonTensor"; 569 std::string op_t_1d_name = (is_all_tensor || !is_tensor) ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name); 570 571 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; 572 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_e" 573 << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 574 } 575 break; 576 case CEED_EVAL_GRAD: 577 if (is_at_points) { 578 std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; 579 580 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; 581 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix 582 << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; 583 } else if (use_3d_slices) { 584 std::string function_name = 585 (dim > 1 ? "InterpTensor" : "Interp") + std::string(is_collocated ? "CollocatedNodes" : "") + std::to_string(dim) + "d"; 586 587 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 588 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix 589 << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 590 } else if (is_tensor) { 591 bool is_collocated_grad = dim == 3 && Q_1d >= P_1d; 592 std::string function_name = 593 (dim == 1 ? "Grad" : ("GradTensor" + std::string(is_collocated ? "CollocatedNodes" : (is_collocated_grad ? "Collocated" : "")))) + 594 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened"); 595 std::string op_t_1d_name = is_all_tensor ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name); 596 597 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "*" 598 << (is_all_tensor && dim >= 3 ? Q_name : "1") << "];\n"; 599 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_e" 600 << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; 601 } else { 602 std::string function_name = "GradNonTensor"; 603 604 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n"; 605 code << tab << function_name << "<num_comp" << var_suffix << ", dim" << var_suffix << ", " << P_name << ", " << Q_name 606 << ", OP_T_1D>(data, r_e" << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; 607 } 608 break; 609 case CEED_EVAL_WEIGHT: { 610 if (is_at_points) { 611 code << tab << "// Nothing to do AtPoints\n"; 612 } else { 613 CeedBasis_Cuda_shared *basis_data; 614 std::string function_name = is_tensor 615 ? ((dim == 1 ? "Weight" : "WeightTensor") + std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened")) 616 : "WeightNonTensor"; 617 618 code << tab << "CeedScalar r_q" << var_suffix << "[" << (is_all_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; 619 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 620 data->W = basis_data->d_q_weight_1d; 621 code << tab << function_name << "<" << P_name << ", " << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; 622 } 623 break; 624 } 625 // LCOV_EXCL_START 626 case CEED_EVAL_DIV: 627 case CEED_EVAL_CURL: 628 break; // TODO: Not implemented 629 // LCOV_EXCL_STOP 630 } 631 } else { 632 switch (eval_mode) { 633 case CEED_EVAL_NONE: 634 code << tab << "CeedScalar *r_e" << var_suffix << " = r_q" << var_suffix << ";\n"; 635 break; // No action 636 case CEED_EVAL_INTERP: 637 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 638 if (is_at_points) { 639 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; 640 641 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_c" << var_suffix 642 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 643 } else { 644 std::string function_name = 645 is_tensor ? ((dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::string(is_collocated ? "CollocatedNodes" : "") + 646 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened")) 647 : "InterpTransposeNonTensor"; 648 std::string op_t_1d_name = (is_all_tensor || !is_tensor) ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name); 649 650 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_q" 651 << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 652 } 653 break; 654 case CEED_EVAL_GRAD: 655 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 656 if (is_at_points) { 657 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; 658 659 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_c" << var_suffix 660 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 661 } else if (use_3d_slices) { 662 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::string(is_collocated ? "CollocatedNodes" : "") + 663 std::to_string(dim) + "d"; 664 665 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_q" << var_suffix 666 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 667 } else if (is_tensor) { 668 bool is_collocated_grad = dim == 3 && Q_1d >= P_1d; 669 std::string function_name = 670 (dim == 1 ? "GradTranspose" 671 : ("GradTransposeTensor" + std::string(is_collocated ? "CollocatedNodes" : (is_collocated_grad ? "Collocated" : "")))) + 672 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened"); 673 std::string op_t_1d_name = is_all_tensor ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name); 674 675 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_q" 676 << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; 677 } else { 678 std::string function_name = "GradTransposeNonTensor"; 679 680 code << tab << function_name << "<num_comp" << var_suffix << ", dim" << var_suffix << ", " << P_name << ", " << Q_name 681 << ", OP_T_1D>(data, r_q" << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; 682 } 683 break; 684 // LCOV_EXCL_START 685 case CEED_EVAL_WEIGHT: 686 break; // Should not occur 687 case CEED_EVAL_DIV: 688 case CEED_EVAL_CURL: 689 break; // TODO: Not implemented 690 // LCOV_EXCL_STOP 691 } 692 } 693 CeedCallBackend(CeedBasisDestroy(&basis)); 694 return CEED_ERROR_SUCCESS; 695 } 696 697 //------------------------------------------------------------------------------ 698 // QFunction 699 //------------------------------------------------------------------------------ 700 static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt max_dim, 701 CeedInt max_num_points, CeedInt num_input_fields, CeedOperatorField *op_input_fields, 702 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, 703 CeedOperatorField *op_output_fields, CeedQFunctionField *qf_output_fields, 704 std::string qfunction_name, CeedInt Q_1d, bool is_all_tensor, bool is_at_points, 705 bool use_3d_slices) { 706 std::string Q_name = is_all_tensor ? "Q_1d" : "Q"; 707 CeedEvalMode eval_mode = CEED_EVAL_NONE; 708 CeedElemRestriction elem_rstr; 709 710 // Setup output arrays 711 code << "\n"; 712 code << tab << "// -- Output field setup\n"; 713 for (CeedInt i = 0; i < num_output_fields; i++) { 714 const char *field_name; 715 std::string var_suffix = "_out_" + std::to_string(i); 716 717 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 718 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 719 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 720 switch (eval_mode) { 721 case CEED_EVAL_NONE: 722 if (is_at_points) { 723 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n"; 724 } else { 725 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (max_dim >= 3) ? Q_name : "1") 726 << "];\n"; 727 } 728 break; 729 case CEED_EVAL_INTERP: 730 if (is_at_points) { 731 // Accumulator for point data 732 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "];\n"; 733 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "; i++) r_c" << var_suffix 734 << "[i] = 0.0;\n"; 735 } else { 736 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (max_dim >= 3) ? Q_name : "1") 737 << "];\n"; 738 } 739 break; 740 case CEED_EVAL_GRAD: 741 if (is_at_points) { 742 // Accumulator for point data 743 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "];\n"; 744 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "; i++) r_c" << var_suffix 745 << "[i] = 0.0;\n"; 746 } else if (use_3d_slices) { 747 // Accumulator for gradient slices 748 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 749 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) r_q" << var_suffix << "[i] = 0.0;\n"; 750 } else { 751 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "*" 752 << (is_all_tensor && (max_dim >= 3) ? Q_name : "1") << "];\n"; 753 } 754 break; 755 case CEED_EVAL_WEIGHT: 756 break; 757 // LCOV_EXCL_START 758 case CEED_EVAL_DIV: 759 case CEED_EVAL_CURL: 760 break; // TODO: Not implemented 761 // LCOV_EXCL_STOP 762 } 763 } 764 765 if (is_at_points) { 766 // We need to handle batches of points 767 code << "\n"; 768 code << tab << "// Note: Using batches of points\n"; 769 code << tab << "const CeedInt point_loop_bound = (blockDim.x*blockDim.y) * ceil((1.0*max_num_points) / (blockDim.x*blockDim.y));\n\n"; 770 code << tab << "#pragma unroll\n"; 771 code << tab << "for (CeedInt i = threadIdx.x + threadIdx.y*blockDim.x; i < point_loop_bound; i += blockDim.x*blockDim.y) {\n"; 772 tab.push(); 773 code << tab << "const CeedInt p = i % max_num_points;\n\n"; 774 775 code << tab << "// -- Coordinates\n"; 776 code << tab << "CeedScalar r_x[max_dim];\n"; 777 code << tab << "ReadPoint<max_dim, coords_comp_stride, max_num_points>(data, elem, p, max_num_points, points.indices, points.coords, r_x);\n\n"; 778 779 code << tab << "// -- Input fields\n"; 780 for (CeedInt i = 0; i < num_input_fields; i++) { 781 const char *field_name; 782 std::string var_suffix = "_in_" + std::to_string(i); 783 std::string P_name = "P_1d" + var_suffix; 784 785 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name)); 786 code << tab << "// ---- Input field " << i << ": " << field_name << "\n"; 787 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 788 // Basis action 789 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 790 switch (eval_mode) { 791 case CEED_EVAL_NONE: 792 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 793 code << tab << "ReadPoint<num_comp" << var_suffix << ", comp_stride" << var_suffix 794 << ", max_num_points>(data, elem, p, max_num_points, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; 795 break; 796 case CEED_EVAL_INTERP: 797 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 798 code << tab << "InterpAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name 799 << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n"; 800 break; 801 case CEED_EVAL_GRAD: 802 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n"; 803 code << tab << "GradAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name 804 << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n"; 805 break; 806 case CEED_EVAL_WEIGHT: 807 code << tab << "CeedScalar r_s" << var_suffix << "[1];\n"; 808 code << tab << "r_s" << var_suffix << "[0] = 1.0;\n"; 809 break; 810 // LCOV_EXCL_START 811 case CEED_EVAL_DIV: 812 case CEED_EVAL_CURL: 813 break; // TODO: Not implemented 814 // LCOV_EXCL_STOP 815 } 816 } 817 code << "\n"; 818 code << tab << "// -- Output fields\n"; 819 for (CeedInt i = 0; i < num_output_fields; i++) { 820 const char *field_name; 821 std::string var_suffix = "_out_" + std::to_string(i); 822 823 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 824 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 825 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 826 // Basis action 827 switch (eval_mode) { 828 case CEED_EVAL_NONE: 829 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 830 break; 831 case CEED_EVAL_INTERP: 832 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 833 break; 834 case CEED_EVAL_GRAD: 835 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n"; 836 break; 837 // LCOV_EXCL_START 838 case CEED_EVAL_WEIGHT: 839 break; // Should not occur 840 case CEED_EVAL_DIV: 841 case CEED_EVAL_CURL: 842 break; // TODO: Not implemented 843 // LCOV_EXCL_STOP 844 } 845 } 846 847 } else if (use_3d_slices) { 848 // We treat quadrature points per slice in 3d to save registers 849 code << "\n"; 850 code << tab << "// Note: Using planes of 3D elements\n"; 851 code << tab << "#pragma unroll\n"; 852 code << tab << "for (CeedInt q = 0; q < " << Q_name << "; q++) {\n"; 853 tab.push(); 854 code << tab << "// -- Input fields\n"; 855 for (CeedInt i = 0; i < num_input_fields; i++) { 856 const char *field_name; 857 std::string var_suffix = "_in_" + std::to_string(i); 858 859 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name)); 860 code << tab << "// ---- Input field " << i << ": " << field_name << "\n"; 861 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 862 // Basis action 863 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 864 switch (eval_mode) { 865 case CEED_EVAL_NONE: 866 bool is_strided; 867 868 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 869 870 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 871 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 872 if (is_strided) { 873 bool has_backend_strides; 874 CeedInt num_elem, elem_size; 875 876 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 877 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 878 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 879 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 880 881 if (!has_backend_strides) { 882 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 883 } 884 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1] 885 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n"; 886 code << tab << "ReadEVecSliceStrided3d<num_comp" << var_suffix << ", " << Q_name << ", strides" << var_suffix << "_0, strides" 887 << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n"; 888 } else { 889 CeedSize l_size = 0; 890 CeedInt comp_stride; 891 CeedElemRestriction_Cuda *rstr_data; 892 893 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 894 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 895 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 896 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 897 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 898 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 899 code << tab << "ReadEVecSliceStandard3d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", " << Q_name << ">(data, l_size" 900 << var_suffix << ", elem, q, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; 901 } 902 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 903 break; 904 case CEED_EVAL_INTERP: 905 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 906 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) {\n"; 907 tab.push(); 908 code << tab << "r_s" << var_suffix << "[j] = r_q" << var_suffix << "[q + j*" << Q_name << "];\n"; 909 tab.pop(); 910 code << tab << "}\n"; 911 break; 912 case CEED_EVAL_GRAD: 913 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n"; 914 code << tab << "GradColloSlice3d<num_comp" << var_suffix << ", " << Q_name << ", OP_T_1D>(data, q, r_q" << var_suffix << ", s_G" 915 << var_suffix << ", r_s" << var_suffix << ");\n"; 916 break; 917 case CEED_EVAL_WEIGHT: 918 code << tab << "CeedScalar r_s" << var_suffix << "[1];\n"; 919 code << tab << "r_s" << var_suffix << "[0] = r_q" << var_suffix << "[q];\n"; 920 break; 921 // LCOV_EXCL_START 922 case CEED_EVAL_DIV: 923 case CEED_EVAL_CURL: 924 break; // TODO: Not implemented 925 // LCOV_EXCL_STOP 926 } 927 } 928 code << "\n"; 929 code << tab << "// -- Output fields\n"; 930 for (CeedInt i = 0; i < num_output_fields; i++) { 931 const char *field_name; 932 std::string var_suffix = "_out_" + std::to_string(i); 933 934 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 935 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 936 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 937 // Basis action 938 switch (eval_mode) { 939 case CEED_EVAL_NONE: 940 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 941 break; 942 case CEED_EVAL_INTERP: 943 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 944 break; 945 case CEED_EVAL_GRAD: 946 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n"; 947 break; 948 // LCOV_EXCL_START 949 case CEED_EVAL_WEIGHT: 950 break; // Should not occur 951 case CEED_EVAL_DIV: 952 case CEED_EVAL_CURL: 953 break; // TODO: Not implemented 954 // LCOV_EXCL_STOP 955 } 956 } 957 } else { 958 code << "\n"; 959 code << tab << "// Note: Using full elements\n"; 960 code << tab << "{\n"; 961 tab.push(); 962 code << tab << "// -- Input fields\n"; 963 for (CeedInt i = 0; i < num_input_fields; i++) { 964 const char *field_name; 965 966 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name)); 967 code << tab << "// ---- Input field " << i << ": " << field_name << "\n"; 968 code << tab << "CeedScalar *r_s_in_" << i << " = r_q_in_" << i << ";\n"; 969 } 970 code << tab << "// -- Output fields\n"; 971 for (CeedInt i = 0; i < num_output_fields; i++) { 972 const char *field_name; 973 974 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 975 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 976 code << tab << "CeedScalar *r_s_out_" << i << " = r_q_out_" << i << ";\n"; 977 } 978 } 979 980 // Input and output buffers 981 code << "\n"; 982 code << tab << "// -- QFunction inputs and outputs\n"; 983 code << tab << "// ---- Inputs\n"; 984 code << tab << "CeedScalar *inputs[" << CeedIntMax(num_input_fields, 1) << "];\n"; 985 for (CeedInt i = 0; i < num_input_fields; i++) { 986 const char *field_name; 987 988 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name)); 989 code << tab << "// ------ Input field " << i << ": " << field_name << "\n"; 990 code << tab << "inputs[" << i << "] = r_s_in_" << i << ";\n"; 991 } 992 code << tab << "// ---- Outputs\n"; 993 code << tab << "CeedScalar *outputs[" << CeedIntMax(num_output_fields, 1) << "];\n"; 994 for (CeedInt i = 0; i < num_output_fields; i++) { 995 const char *field_name; 996 997 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 998 code << tab << "// ------ Output field " << i << ": " << field_name << "\n"; 999 code << tab << "outputs[" << i << "] = r_s_out_" << i << ";\n"; 1000 } 1001 1002 // Apply QFunction 1003 code << "\n"; 1004 code << tab << "// -- Apply QFunction\n"; 1005 code << tab << "" << qfunction_name << "(ctx, "; 1006 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) { 1007 code << "1"; 1008 } else { 1009 code << Q_name; 1010 } 1011 code << ", inputs, outputs);\n"; 1012 1013 if (is_at_points) { 1014 // Map back to coefficients 1015 code << "\n"; 1016 code << tab << "// -- Output fields\n"; 1017 for (CeedInt i = 0; i < num_output_fields; i++) { 1018 const char *field_name; 1019 std::string var_suffix = "_out_" + std::to_string(i); 1020 std::string P_name = "P_1d" + var_suffix; 1021 1022 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 1023 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 1024 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1025 // Basis action 1026 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 1027 switch (eval_mode) { 1028 case CEED_EVAL_NONE: { 1029 CeedInt comp_stride; 1030 CeedElemRestriction elem_rstr; 1031 1032 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1033 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 1034 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1035 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 1036 code << tab << "WritePoint<num_comp" << var_suffix << ", comp_stride" << var_suffix 1037 << ", max_num_points>(data, elem, i, points.num_per_elem[elem], indices.outputs[" << i << "]" 1038 << ", r_s" << var_suffix << ", d" << var_suffix << ");\n"; 1039 break; 1040 } 1041 case CEED_EVAL_INTERP: 1042 code << tab << "if (i >= points.num_per_elem[elem]) {\n"; 1043 tab.push(); 1044 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n"; 1045 tab.pop(); 1046 code << tab << "}\n"; 1047 code << tab << "InterpTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name 1048 << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; 1049 break; 1050 case CEED_EVAL_GRAD: 1051 code << tab << "if (i >= points.num_per_elem[elem]) {\n"; 1052 tab.push(); 1053 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "*dim" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n"; 1054 tab.pop(); 1055 code << tab << "}\n"; 1056 code << tab << "GradTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name 1057 << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; 1058 break; 1059 // LCOV_EXCL_START 1060 case CEED_EVAL_WEIGHT: 1061 break; // Should not occur 1062 case CEED_EVAL_DIV: 1063 case CEED_EVAL_CURL: 1064 break; // TODO: Not implemented 1065 // LCOV_EXCL_STOP 1066 } 1067 } 1068 } else if (use_3d_slices) { 1069 // Copy or apply transpose grad, if needed 1070 code << "\n"; 1071 code << tab << "// -- Output fields\n"; 1072 for (CeedInt i = 0; i < num_output_fields; i++) { 1073 const char *field_name; 1074 std::string var_suffix = "_out_" + std::to_string(i); 1075 std::string P_name = "P_1d" + var_suffix; 1076 1077 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 1078 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 1079 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1080 // Basis action 1081 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 1082 switch (eval_mode) { 1083 case CEED_EVAL_NONE: 1084 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 1085 tab.push(); 1086 code << tab << "r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 1087 tab.pop(); 1088 code << tab << "}\n"; 1089 break; 1090 case CEED_EVAL_INTERP: 1091 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 1092 tab.push(); 1093 code << tab << "r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 1094 tab.pop(); 1095 code << tab << "}\n"; 1096 break; 1097 case CEED_EVAL_GRAD: 1098 code << tab << "GradColloSliceTranspose3d<num_comp" << var_suffix << ", " << Q_name << ", OP_T_1D>(data, q, r_s" << var_suffix << ", s_G" 1099 << var_suffix << ", r_q" << var_suffix << ");\n"; 1100 break; 1101 // LCOV_EXCL_START 1102 case CEED_EVAL_WEIGHT: 1103 break; // Should not occur 1104 case CEED_EVAL_DIV: 1105 case CEED_EVAL_CURL: 1106 break; // TODO: Not implemented 1107 // LCOV_EXCL_STOP 1108 } 1109 } 1110 } 1111 tab.pop(); 1112 code << tab << "}\n"; 1113 return CEED_ERROR_SUCCESS; 1114 } 1115 1116 //------------------------------------------------------------------------------ 1117 // Build single operator kernel 1118 //------------------------------------------------------------------------------ 1119 extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_build) { 1120 bool is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false; 1121 Ceed ceed; 1122 CeedInt Q = 0, Q_1d = 0, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0; 1123 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1124 CeedQFunction_Cuda_gen *qf_data; 1125 CeedQFunction qf; 1126 CeedOperatorField *op_input_fields, *op_output_fields; 1127 CeedOperator_Cuda_gen *data; 1128 std::ostringstream code; 1129 Tab tab; 1130 1131 CeedCallBackend(CeedOperatorGetData(op, &data)); 1132 { 1133 bool is_setup_done; 1134 1135 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 1136 if (is_setup_done) { 1137 *is_good_build = !data->use_fallback; 1138 return CEED_ERROR_SUCCESS; 1139 } 1140 } 1141 1142 // Check field compatibility 1143 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1144 { 1145 bool has_shared_bases = true; 1146 1147 for (CeedInt i = 0; i < num_input_fields; i++) { 1148 CeedBasis basis; 1149 1150 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1151 if (basis != CEED_BASIS_NONE) { 1152 bool is_tensor = true; 1153 const char *resource; 1154 char *resource_root; 1155 Ceed basis_ceed; 1156 1157 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 1158 is_all_tensor = is_all_tensor && is_tensor; 1159 is_all_nontensor = is_all_nontensor && !is_tensor; 1160 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); 1161 CeedCallBackend(CeedGetResource(basis_ceed, &resource)); 1162 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); 1163 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared"); 1164 CeedCallBackend(CeedFree(&resource_root)); 1165 CeedCallBackend(CeedDestroy(&basis_ceed)); 1166 } 1167 CeedCallBackend(CeedBasisDestroy(&basis)); 1168 } 1169 1170 for (CeedInt i = 0; i < num_output_fields; i++) { 1171 CeedBasis basis; 1172 1173 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1174 if (basis != CEED_BASIS_NONE) { 1175 bool is_tensor = true; 1176 const char *resource; 1177 char *resource_root; 1178 Ceed basis_ceed; 1179 1180 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 1181 is_all_tensor = is_all_tensor && is_tensor; 1182 is_all_nontensor = is_all_nontensor && !is_tensor; 1183 1184 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); 1185 CeedCallBackend(CeedGetResource(basis_ceed, &resource)); 1186 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); 1187 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared"); 1188 CeedCallBackend(CeedFree(&resource_root)); 1189 CeedCallBackend(CeedDestroy(&basis_ceed)); 1190 } 1191 CeedCallBackend(CeedBasisDestroy(&basis)); 1192 } 1193 // -- Fallback to ref if not all bases are shared 1194 if (!has_shared_bases) { 1195 *is_good_build = false; 1196 return CEED_ERROR_SUCCESS; 1197 } 1198 } 1199 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1200 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1201 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 1202 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1203 1204 // Get operator data 1205 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); 1206 { 1207 CeedInt max_P = 0, max_P_1d = 0; 1208 1209 CeedCallBackend(CeedOperatorBuildKernelData_Cuda_gen(ceed, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, 1210 op_output_fields, qf_output_fields, &max_P, &max_P_1d, &Q, &Q_1d, &max_dim, &is_all_tensor, 1211 &use_3d_slices)); 1212 data->max_P_1d = is_all_tensor ? max_P_1d : max_P; 1213 } 1214 if (is_at_points) { 1215 CeedInt coords_dim = 0; 1216 CeedElemRestriction_Cuda *rstr_data; 1217 CeedElemRestriction rstr_points = NULL; 1218 1219 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 1220 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 1221 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride)); 1222 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &coords_dim)); 1223 CeedCallBackend(CeedElemRestrictionGetData(rstr_points, &rstr_data)); 1224 data->points.indices = (CeedInt *)rstr_data->d_offsets; 1225 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1226 if (max_dim == 0) max_dim = coords_dim; 1227 if (Q_1d == 0) max_num_points = ceil(pow(max_num_points, 1.0 / max_dim)); 1228 } 1229 if (max_dim == 0) max_dim = 1; 1230 data->dim = max_dim; 1231 if (is_at_points) use_3d_slices = false; 1232 if (Q_1d == 0) { 1233 if (is_at_points) Q_1d = max_num_points; 1234 else CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q_1d)); 1235 } 1236 if (Q == 0) Q = Q_1d; 1237 data->Q = Q; 1238 data->Q_1d = Q_1d; 1239 1240 // Check for restriction only identity operator 1241 { 1242 bool is_identity_qf; 1243 1244 CeedCallBackend(CeedQFunctionIsIdentity(qf, &is_identity_qf)); 1245 if (is_identity_qf) { 1246 CeedEvalMode eval_mode_in, eval_mode_out; 1247 1248 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[0], &eval_mode_in)); 1249 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[0], &eval_mode_out)); 1250 CeedCheck(eval_mode_in != CEED_EVAL_NONE || eval_mode_out != CEED_EVAL_NONE, ceed, CEED_ERROR_BACKEND, 1251 "Backend does not implement restriction only identity operators"); 1252 } 1253 } 1254 1255 // Add atomicAdd function for old NVidia architectures 1256 { 1257 Ceed_Cuda *ceed_data; 1258 struct cudaDeviceProp prop; 1259 1260 CeedCallBackend(CeedGetData(ceed, &ceed_data)); 1261 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id)); 1262 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) { 1263 code << tab << "// AtomicAdd fallback source\n"; 1264 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n"; 1265 } 1266 } 1267 1268 // Load basis source files 1269 if (!is_all_nontensor) { 1270 code << tab << "// Tensor basis source\n"; 1271 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n"; 1272 } 1273 if (!is_all_tensor) { 1274 code << tab << "// Non-tensor basis source\n"; 1275 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-nontensor-templates.h>\n\n"; 1276 } 1277 if (!is_all_tensor && !is_all_nontensor) { 1278 code << "// Tensor basis source\n"; 1279 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-flattened-templates.h>\n\n"; 1280 } 1281 if (is_at_points) { 1282 code << "// AtPoints basis source\n"; 1283 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h>\n\n"; 1284 } 1285 code << "// CodeGen operator source\n"; 1286 code << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n"; 1287 1288 // Get QFunction name 1289 std::string qfunction_name(qf_data->qfunction_name); 1290 std::string operator_name; 1291 1292 operator_name = "CeedKernelCudaGenOperator_" + qfunction_name; 1293 1294 // Define CEED_Q_VLA 1295 code << "\n" << tab << "#undef CEED_Q_VLA\n"; 1296 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) { 1297 code << tab << "#define CEED_Q_VLA 1\n\n"; 1298 } else { 1299 code << tab << "#define CEED_Q_VLA " << Q_1d << "\n\n"; 1300 } 1301 1302 // Add user QFunction source 1303 { 1304 const char *source_path; 1305 1306 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path)); 1307 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file"); 1308 1309 code << tab << "// User QFunction source\n"; 1310 code << tab << "#include \"" << source_path << "\"\n\n"; 1311 } 1312 1313 // Setup 1314 code << "\n" << tab << "// -----------------------------------------------------------------------------\n"; 1315 code << tab << "// Operator Kernel\n"; 1316 code << tab << "// \n"; 1317 code << tab << "// d_[in,out]_i: CeedVector device array\n"; 1318 code << tab << "// r_[in,out]_e_i: Element vector register\n"; 1319 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n"; 1320 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n"; 1321 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; 1322 code << tab << "// \n"; 1323 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; 1324 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; 1325 code << tab << "// -----------------------------------------------------------------------------\n"; 1326 code << tab << "extern \"C\" __global__ void " << operator_name 1327 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda " 1328 "points) {\n"; 1329 tab.push(); 1330 1331 // Scratch buffers 1332 for (CeedInt i = 0; i < num_input_fields; i++) { 1333 CeedEvalMode eval_mode; 1334 1335 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1336 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 1337 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n"; 1338 } 1339 } 1340 for (CeedInt i = 0; i < num_output_fields; i++) { 1341 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n"; 1342 } 1343 1344 code << tab << "const CeedInt max_dim = " << max_dim << ";\n"; 1345 if (!is_all_tensor) { 1346 code << tab << "const CeedInt Q = " << Q << ";\n"; 1347 } 1348 if (!is_all_nontensor) { 1349 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n"; 1350 } 1351 if (is_at_points) { 1352 code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n"; 1353 code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; 1354 } 1355 1356 // Shared data 1357 code << tab << "extern __shared__ CeedScalar slice[];\n"; 1358 code << tab << "SharedData_Cuda data;\n"; 1359 code << tab << "data.t_id_x = threadIdx.x;\n"; 1360 code << tab << "data.t_id_y = threadIdx.y;\n"; 1361 code << tab << "data.t_id_z = threadIdx.z;\n"; 1362 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 1363 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n"; 1364 1365 // -- Determine input mat reuse 1366 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX]; 1367 1368 for (CeedInt i = 0; i < num_input_fields; i++) { 1369 input_matrix_reuse[i].index = -1; 1370 } 1371 for (CeedInt i = 0; i < num_input_fields; i++) { 1372 bool is_tensor = true; 1373 CeedEvalMode eval_mode_i; 1374 CeedBasis basis_i; 1375 1376 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i)); 1377 if (eval_mode_i == CEED_EVAL_WEIGHT) continue; 1378 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i)); 1379 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor)); 1380 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) { 1381 CeedEvalMode eval_mode_j; 1382 CeedBasis basis_j; 1383 1384 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 1385 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1386 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 1387 if (basis_i == basis_j) { 1388 if (is_tensor) { 1389 input_matrix_reuse[i].index = j; 1390 input_matrix_reuse[i].is_input = true; 1391 input_matrix_reuse[i].eval_mode = eval_mode_j; 1392 } else { 1393 // For non-tensor can only re-use with the same eval mode 1394 if (eval_mode_i == eval_mode_j) { 1395 input_matrix_reuse[i].index = j; 1396 input_matrix_reuse[i].is_input = true; 1397 input_matrix_reuse[i].eval_mode = eval_mode_j; 1398 } 1399 } 1400 } 1401 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1402 } 1403 CeedCallBackend(CeedBasisDestroy(&basis_i)); 1404 } 1405 1406 // -- Determine output mat reuse 1407 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX]; 1408 1409 for (CeedInt i = 0; i < num_output_fields; i++) { 1410 output_matrix_reuse[i].index = -1; 1411 } 1412 for (CeedInt i = 0; i < num_output_fields; i++) { 1413 bool is_tensor = true; 1414 CeedEvalMode eval_mode_i; 1415 CeedBasis basis_i; 1416 1417 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i)); 1418 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i)); 1419 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor)); 1420 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) { 1421 CeedEvalMode eval_mode_j; 1422 CeedBasis basis_j; 1423 1424 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 1425 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1426 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 1427 if (basis_i == basis_j) { 1428 if (is_tensor) { 1429 output_matrix_reuse[i].index = j; 1430 output_matrix_reuse[i].is_input = true; 1431 output_matrix_reuse[i].eval_mode = eval_mode_j; 1432 } else { 1433 // For non-tensor can only re-use with the same eval mode 1434 if (eval_mode_i == eval_mode_j) { 1435 output_matrix_reuse[i].index = j; 1436 output_matrix_reuse[i].is_input = true; 1437 output_matrix_reuse[i].eval_mode = eval_mode_j; 1438 } 1439 } 1440 } 1441 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1442 } 1443 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) { 1444 CeedEvalMode eval_mode_j; 1445 CeedBasis basis_j; 1446 1447 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j)); 1448 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1449 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j)); 1450 if (basis_i == basis_j) { 1451 if (is_tensor) { 1452 output_matrix_reuse[i].index = j; 1453 output_matrix_reuse[i].is_input = false; 1454 output_matrix_reuse[i].eval_mode = eval_mode_j; 1455 } else { 1456 // For non-tensor can only re-use with the same eval mode 1457 if (eval_mode_i == eval_mode_j) { 1458 output_matrix_reuse[i].index = j; 1459 output_matrix_reuse[i].is_input = false; 1460 output_matrix_reuse[i].eval_mode = eval_mode_j; 1461 } 1462 } 1463 } 1464 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1465 } 1466 CeedCallBackend(CeedBasisDestroy(&basis_i)); 1467 } 1468 1469 // Initialize constants, and matrices B and G 1470 code << "\n" << tab << "// Input field constants and basis data\n"; 1471 for (CeedInt i = 0; i < num_input_fields; i++) { 1472 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i], 1473 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, false)); 1474 } 1475 code << "\n" << tab << "// Output field constants and basis data\n"; 1476 for (CeedInt i = 0; i < num_output_fields; i++) { 1477 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i], 1478 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, false)); 1479 } 1480 1481 // Loop over all elements 1482 code << "\n" << tab << "// Element loop\n"; 1483 code << tab << "__syncthreads();\n"; 1484 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; 1485 tab.push(); 1486 1487 // -- Compute minimum buffer space needed 1488 CeedInt max_rstr_buffer_size = 1; 1489 1490 for (CeedInt i = 0; i < num_input_fields; i++) { 1491 CeedEvalMode eval_mode; 1492 1493 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1494 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) { 1495 CeedInt num_comp; 1496 CeedElemRestriction elem_rstr; 1497 1498 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1499 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1500 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 1501 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1502 } 1503 } 1504 for (CeedInt i = 0; i < num_output_fields; i++) { 1505 CeedEvalMode eval_mode; 1506 1507 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1508 if (eval_mode != CEED_EVAL_NONE) { 1509 CeedInt num_comp; 1510 CeedElemRestriction elem_rstr; 1511 1512 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1513 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1514 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 1515 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1516 } 1517 } 1518 code << tab << "// Scratch restriction buffer space\n"; 1519 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n"; 1520 1521 // -- Determine best input field processing order 1522 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX]; 1523 1524 for (CeedInt i = 0; i < num_input_fields; i++) { 1525 field_rstr_in_buffer[i] = -1; 1526 input_field_order[i] = -1; 1527 } 1528 { 1529 bool is_ordered[CEED_FIELD_MAX]; 1530 CeedInt curr_index = 0; 1531 1532 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false; 1533 for (CeedInt i = 0; i < num_input_fields; i++) { 1534 CeedVector vec_i; 1535 CeedElemRestriction rstr_i; 1536 1537 if (is_ordered[i]) continue; 1538 field_rstr_in_buffer[i] = i; 1539 is_ordered[i] = true; 1540 input_field_order[curr_index] = i; 1541 curr_index++; 1542 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i)); 1543 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT 1544 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i)); 1545 for (CeedInt j = i + 1; j < num_input_fields; j++) { 1546 CeedVector vec_j; 1547 CeedElemRestriction rstr_j; 1548 1549 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j)); 1550 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j)); 1551 if (rstr_i == rstr_j && vec_i == vec_j) { 1552 field_rstr_in_buffer[j] = i; 1553 is_ordered[j] = true; 1554 input_field_order[curr_index] = j; 1555 curr_index++; 1556 } 1557 CeedCallBackend(CeedVectorDestroy(&vec_j)); 1558 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j)); 1559 } 1560 CeedCallBackend(CeedVectorDestroy(&vec_i)); 1561 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i)); 1562 } 1563 } 1564 1565 // -- Input restriction and basis 1566 code << "\n" << tab << "// -- Input field restrictions and basis actions\n"; 1567 for (CeedInt i = 0; i < num_input_fields; i++) { 1568 const char *field_name; 1569 const CeedInt f = input_field_order[i]; 1570 1571 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name)); 1572 code << tab << "// ---- Input field " << f << ": " << field_name << "\n"; 1573 1574 // ---- Restriction 1575 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], 1576 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices)); 1577 1578 // ---- Basis action 1579 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true, 1580 is_all_tensor, is_at_points, use_3d_slices)); 1581 } 1582 1583 // -- Q function 1584 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields, 1585 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name, 1586 Q_1d, is_all_tensor, is_at_points, use_3d_slices)); 1587 1588 // -- Output basis and restriction 1589 code << "\n" << tab << "// -- Output field basis action and restrictions\n"; 1590 for (CeedInt i = 0; i < num_output_fields; i++) { 1591 const char *field_name; 1592 1593 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 1594 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 1595 1596 // ---- Basis action 1597 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d, false, 1598 is_all_tensor, is_at_points, use_3d_slices)); 1599 1600 // ---- Restriction 1601 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, i, NULL, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d, 1602 false, is_all_tensor, is_at_points, use_3d_slices)); 1603 } 1604 1605 // Close loop and function 1606 tab.pop(); 1607 code << tab << "}\n"; 1608 tab.pop(); 1609 code << tab << "}\n"; 1610 code << tab << "// -----------------------------------------------------------------------------\n\n"; 1611 1612 // Compile 1613 { 1614 bool is_compile_good = false; 1615 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d); 1616 1617 data->thread_1d = T_1d; 1618 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good, &data->module, 1, "OP_T_1D", T_1d)); 1619 if (is_compile_good) { 1620 *is_good_build = true; 1621 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, operator_name.c_str(), &data->op)); 1622 } else { 1623 *is_good_build = false; 1624 data->use_fallback = true; 1625 } 1626 } 1627 CeedCallBackend(CeedOperatorSetSetupDone(op)); 1628 CeedCallBackend(CeedDestroy(&ceed)); 1629 CeedCallBackend(CeedQFunctionDestroy(&qf)); 1630 return CEED_ERROR_SUCCESS; 1631 } 1632 1633 //------------------------------------------------------------------------------ 1634 // Build AtPoints assembly operator kernel 1635 //------------------------------------------------------------------------------ 1636 static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, bool is_full, bool *is_good_build) { 1637 bool is_all_tensor = true, is_at_points = false, use_3d_slices = false; 1638 Ceed ceed; 1639 CeedInt Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0; 1640 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1641 CeedQFunction_Cuda_gen *qf_data; 1642 CeedQFunction qf; 1643 CeedOperatorField *op_input_fields, *op_output_fields; 1644 CeedOperator_Cuda_gen *data; 1645 std::ostringstream code; 1646 Tab tab; 1647 1648 // Check compatibility 1649 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1650 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); 1651 CeedCheck(is_at_points, ceed, CEED_ERROR_BACKEND, "Only AtPoints operator assembly supported"); 1652 1653 // Retrieve operator data 1654 CeedCallBackend(CeedOperatorGetData(op, &data)); 1655 Q = data->Q; 1656 Q_1d = data->Q_1d; 1657 max_dim = data->dim; 1658 { 1659 CeedElemRestriction rstr_points = NULL; 1660 1661 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 1662 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 1663 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride)); 1664 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1665 } 1666 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1667 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 1668 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1669 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1670 1671 // Add atomicAdd function for old NVidia architectures 1672 { 1673 Ceed_Cuda *ceed_data; 1674 struct cudaDeviceProp prop; 1675 1676 CeedCallBackend(CeedGetData(ceed, &ceed_data)); 1677 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id)); 1678 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) { 1679 code << tab << "// AtomicAdd fallback source\n"; 1680 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n"; 1681 } 1682 } 1683 1684 // Load basis source files 1685 code << tab << "// Tensor basis source\n"; 1686 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n"; 1687 code << tab << "// AtPoints basis source\n"; 1688 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h>\n\n"; 1689 code << tab << "// CodeGen operator source\n"; 1690 code << tab << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n"; 1691 1692 // Get QFunction name 1693 std::string qfunction_name(qf_data->qfunction_name); 1694 std::string operator_name; 1695 1696 if (is_full) { 1697 operator_name = "CeedKernelCudaGenOperatorFullAssembly_" + qfunction_name; 1698 } else { 1699 operator_name = "CeedKernelCudaGenOperatorDiagonalAssembly_" + qfunction_name; 1700 } 1701 1702 // Define CEED_Q_VLA 1703 code << "\n" << tab << "#undef CEED_Q_VLA\n"; 1704 code << tab << "#define CEED_Q_VLA 1\n\n"; 1705 1706 // Add user QFunction source 1707 { 1708 const char *source_path; 1709 1710 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path)); 1711 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file"); 1712 1713 code << tab << "// User QFunction source\n"; 1714 code << tab << "#include \"" << source_path << "\"\n\n"; 1715 } 1716 1717 // Setup 1718 code << "\n" << tab << "// -----------------------------------------------------------------------------\n"; 1719 code << tab << "// Operator Assembly Kernel\n"; 1720 code << tab << "// \n"; 1721 code << tab << "// d_[in,out]_i: CeedVector device array\n"; 1722 code << tab << "// r_[in,out]_e_i: Element vector register\n"; 1723 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n"; 1724 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n"; 1725 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; 1726 code << tab << "// \n"; 1727 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; 1728 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; 1729 code << tab << "// -----------------------------------------------------------------------------\n"; 1730 code << tab << "extern \"C\" __global__ void " << operator_name 1731 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda " 1732 "points, CeedScalar *__restrict__ values_array) {\n"; 1733 tab.push(); 1734 1735 // Scratch buffers 1736 for (CeedInt i = 0; i < num_input_fields; i++) { 1737 CeedEvalMode eval_mode; 1738 1739 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1740 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 1741 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n"; 1742 } 1743 } 1744 for (CeedInt i = 0; i < num_output_fields; i++) { 1745 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n"; 1746 } 1747 1748 code << tab << "const CeedInt max_dim = " << max_dim << ";\n"; 1749 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n"; 1750 code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n"; 1751 code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; 1752 1753 // Shared data 1754 code << tab << "extern __shared__ CeedScalar slice[];\n"; 1755 code << tab << "SharedData_Cuda data;\n"; 1756 code << tab << "data.t_id_x = threadIdx.x;\n"; 1757 code << tab << "data.t_id_y = threadIdx.y;\n"; 1758 code << tab << "data.t_id_z = threadIdx.z;\n"; 1759 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 1760 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n"; 1761 1762 // -- Determine input mat reuse 1763 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX]; 1764 1765 for (CeedInt i = 0; i < num_input_fields; i++) { 1766 input_matrix_reuse[i].index = -1; 1767 } 1768 for (CeedInt i = 0; i < num_input_fields; i++) { 1769 CeedEvalMode eval_mode_i; 1770 CeedBasis basis_i; 1771 1772 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i)); 1773 if (eval_mode_i == CEED_EVAL_WEIGHT) continue; 1774 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i)); 1775 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) { 1776 CeedEvalMode eval_mode_j; 1777 CeedBasis basis_j; 1778 1779 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 1780 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1781 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 1782 if (basis_i == basis_j) { 1783 input_matrix_reuse[i].index = j; 1784 input_matrix_reuse[i].is_input = true; 1785 input_matrix_reuse[i].eval_mode = eval_mode_j; 1786 } 1787 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1788 } 1789 CeedCallBackend(CeedBasisDestroy(&basis_i)); 1790 } 1791 1792 // -- Determine output mat reuse 1793 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX]; 1794 1795 for (CeedInt i = 0; i < num_output_fields; i++) { 1796 output_matrix_reuse[i].index = -1; 1797 } 1798 for (CeedInt i = 0; i < num_output_fields; i++) { 1799 CeedEvalMode eval_mode_i; 1800 CeedBasis basis_i; 1801 1802 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i)); 1803 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i)); 1804 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) { 1805 CeedEvalMode eval_mode_j; 1806 CeedBasis basis_j; 1807 1808 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 1809 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1810 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 1811 if (basis_i == basis_j) { 1812 output_matrix_reuse[i].index = j; 1813 output_matrix_reuse[i].is_input = true; 1814 output_matrix_reuse[i].eval_mode = eval_mode_j; 1815 } 1816 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1817 } 1818 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) { 1819 CeedEvalMode eval_mode_j; 1820 CeedBasis basis_j; 1821 1822 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j)); 1823 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 1824 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j)); 1825 if (basis_i == basis_j) { 1826 output_matrix_reuse[i].index = j; 1827 output_matrix_reuse[i].is_input = false; 1828 output_matrix_reuse[i].eval_mode = eval_mode_j; 1829 } 1830 CeedCallBackend(CeedBasisDestroy(&basis_j)); 1831 } 1832 CeedCallBackend(CeedBasisDestroy(&basis_i)); 1833 } 1834 1835 // Initialize constants, and matrices B and G 1836 code << "\n" << tab << "// Input field constants and basis data\n"; 1837 for (CeedInt i = 0; i < num_input_fields; i++) { 1838 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i], 1839 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, false)); 1840 } 1841 code << "\n" << tab << "// Output field constants and basis data\n"; 1842 for (CeedInt i = 0; i < num_output_fields; i++) { 1843 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i], 1844 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, false)); 1845 } 1846 1847 // Loop over all elements 1848 code << "\n" << tab << "// Element loop\n"; 1849 code << tab << "__syncthreads();\n"; 1850 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; 1851 tab.push(); 1852 1853 // -- Compute minimum buffer space needed 1854 CeedInt max_rstr_buffer_size = 1; 1855 1856 for (CeedInt i = 0; i < num_input_fields; i++) { 1857 CeedEvalMode eval_mode; 1858 1859 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1860 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) { 1861 CeedInt num_comp; 1862 CeedElemRestriction elem_rstr; 1863 1864 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1865 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1866 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 1867 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1868 } 1869 } 1870 for (CeedInt i = 0; i < num_output_fields; i++) { 1871 CeedEvalMode eval_mode; 1872 1873 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1874 if (eval_mode != CEED_EVAL_NONE) { 1875 CeedInt num_comp; 1876 CeedElemRestriction elem_rstr; 1877 1878 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1879 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1880 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 1881 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1882 } 1883 } 1884 code << tab << "// Scratch restriction buffer space\n"; 1885 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n"; 1886 1887 // -- Determine best input field processing order 1888 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX]; 1889 1890 for (CeedInt i = 0; i < num_input_fields; i++) { 1891 field_rstr_in_buffer[i] = -1; 1892 input_field_order[i] = -1; 1893 } 1894 { 1895 bool is_ordered[CEED_FIELD_MAX]; 1896 CeedInt curr_index = 0; 1897 1898 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false; 1899 for (CeedInt i = 0; i < num_input_fields; i++) { 1900 CeedVector vec_i; 1901 CeedElemRestriction rstr_i; 1902 1903 if (is_ordered[i]) continue; 1904 field_rstr_in_buffer[i] = i; 1905 is_ordered[i] = true; 1906 input_field_order[curr_index] = i; 1907 curr_index++; 1908 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i)); 1909 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT 1910 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i)); 1911 for (CeedInt j = i + 1; j < num_input_fields; j++) { 1912 CeedVector vec_j; 1913 CeedElemRestriction rstr_j; 1914 1915 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j)); 1916 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j)); 1917 if (rstr_i == rstr_j && vec_i == vec_j) { 1918 field_rstr_in_buffer[j] = i; 1919 is_ordered[j] = true; 1920 input_field_order[curr_index] = j; 1921 curr_index++; 1922 } 1923 CeedCallBackend(CeedVectorDestroy(&vec_j)); 1924 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j)); 1925 } 1926 CeedCallBackend(CeedVectorDestroy(&vec_i)); 1927 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i)); 1928 } 1929 } 1930 1931 // -- Input restriction and basis 1932 code << "\n" << tab << "// -- Input field restrictions and basis actions\n"; 1933 CeedInt active_field_index = -1; 1934 1935 for (CeedInt i = 0; i < num_input_fields; i++) { 1936 bool is_active = false; 1937 const char *field_name; 1938 const CeedInt f = input_field_order[i]; 1939 1940 { 1941 CeedVector vec; 1942 1943 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec)); 1944 is_active = vec == CEED_VECTOR_ACTIVE; 1945 CeedCallBackend(CeedVectorDestroy(&vec)); 1946 } 1947 1948 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name)); 1949 code << tab << "// ---- Input field " << f << ": " << field_name << "\n"; 1950 1951 if (is_active) { 1952 std::string var_suffix = "_in_" + std::to_string(f); 1953 1954 code << tab << "// Active field - no restriction or basis action here\n"; 1955 if (active_field_index == -1) { 1956 active_field_index = f; 1957 code << tab << "CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? "P_1d" + var_suffix : "1") 1958 << "] = {0.0};\n"; 1959 } else { 1960 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_in_" << active_field_index << ";\n"; 1961 } 1962 } else { 1963 // ---- Restriction 1964 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], 1965 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices)); 1966 1967 // ---- Basis action 1968 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true, 1969 is_all_tensor, is_at_points, use_3d_slices)); 1970 } 1971 } 1972 1973 // -- Loop over active field 1974 std::string active_var_suffix = "_in_" + std::to_string(active_field_index); 1975 1976 code << "\n" << tab << "// Loop over nodes in active field\n"; 1977 code << tab << "for (CeedInt n = 0; n < num_comp" << active_var_suffix << "*P_1d" << active_var_suffix 1978 << (max_dim > 1 ? "*P_1d" + active_var_suffix : "") << (max_dim > 2 ? "*P_1d" + active_var_suffix : "") << "; n++) {\n"; 1979 tab.push(); 1980 1981 // -- Set current active node and component to 1 1982 code << tab << "// Set current active node and component to 1.0\n"; 1983 code << tab << "SetEVecStandard" << max_dim << "d_Single<num_comp" << active_var_suffix << ", P_1d" << active_var_suffix << ">(data, n, 1.0, r_e" 1984 << active_var_suffix << ");\n\n"; 1985 1986 for (CeedInt i = 0; i < num_input_fields; i++) { 1987 bool is_active = false; 1988 const char *field_name; 1989 const CeedInt f = input_field_order[i]; 1990 1991 { 1992 CeedVector vec; 1993 1994 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec)); 1995 is_active = vec == CEED_VECTOR_ACTIVE; 1996 CeedCallBackend(CeedVectorDestroy(&vec)); 1997 } 1998 if (!is_active) continue; 1999 2000 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name)); 2001 code << tab << "// ---- Input field " << f << ": " << field_name << "\n"; 2002 2003 // ---- Basis action 2004 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true, 2005 is_all_tensor, is_at_points, use_3d_slices)); 2006 } 2007 2008 // -- Q function 2009 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields, 2010 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name, 2011 Q_1d, is_all_tensor, is_at_points, use_3d_slices)); 2012 2013 // -- Output basis and restriction 2014 code << "\n" << tab << "// -- Output field basis action and restrictions\n"; 2015 for (CeedInt i = 0; i < num_output_fields; i++) { 2016 bool is_active = false; 2017 const char *field_name; 2018 2019 { 2020 CeedVector vec; 2021 2022 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 2023 is_active = vec == CEED_VECTOR_ACTIVE; 2024 CeedCallBackend(CeedVectorDestroy(&vec)); 2025 } 2026 if (!is_active) continue; 2027 2028 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 2029 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 2030 2031 // ---- Basis action 2032 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d, false, 2033 is_all_tensor, is_at_points, use_3d_slices)); 2034 2035 // ---- Restriction 2036 if (is_full) { 2037 std::string var_suffix = "_out_" + std::to_string(i); 2038 CeedInt comp_stride; 2039 CeedSize l_size; 2040 CeedElemRestriction elem_rstr; 2041 2042 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 2043 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 2044 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 2045 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 2046 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 2047 code << tab << "WriteLVecStandard" << max_dim << "d_Assembly<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", P_1d" + var_suffix 2048 << ">(data, l_size" << var_suffix << ", elem, n, r_e" << var_suffix << ", values_array);\n"; 2049 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 2050 } else { 2051 std::string var_suffix = "_out_" + std::to_string(i); 2052 CeedInt comp_stride; 2053 CeedSize l_size; 2054 CeedElemRestriction elem_rstr; 2055 2056 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 2057 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 2058 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 2059 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 2060 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 2061 code << tab << "WriteLVecStandard" << max_dim << "d_Single<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", P_1d" + var_suffix 2062 << ">(data, l_size" << var_suffix << ", elem, n, indices.outputs[" << i << "], r_e" << var_suffix << ", values_array);\n"; 2063 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 2064 } 2065 } 2066 2067 // -- Reset current active node and component 2068 code << "\n" << tab << "// Reset current active node and component to 0.0\n"; 2069 code << tab << "SetEVecStandard" << max_dim << "d_Single<num_comp" << active_var_suffix << ", P_1d" << active_var_suffix << ">(data, n, 0.0, r_e" 2070 << active_var_suffix << ");\n"; 2071 2072 // -- End of loop over active field 2073 tab.pop(); 2074 code << tab << "}\n"; 2075 2076 // Close loop and function 2077 tab.pop(); 2078 code << tab << "}\n"; 2079 tab.pop(); 2080 code << tab << "}\n"; 2081 code << tab << "// -----------------------------------------------------------------------------\n\n"; 2082 2083 // Compile 2084 { 2085 bool is_compile_good = false; 2086 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d); 2087 2088 data->thread_1d = T_1d; 2089 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good, 2090 is_full ? &data->module_assemble_full : &data->module_assemble_diagonal, 1, "OP_T_1D", T_1d)); 2091 if (is_compile_good) { 2092 *is_good_build = true; 2093 CeedCallBackend(CeedGetKernel_Cuda(ceed, is_full ? data->module_assemble_full : data->module_assemble_diagonal, operator_name.c_str(), 2094 is_full ? &data->assemble_full : &data->assemble_diagonal)); 2095 } else { 2096 *is_good_build = false; 2097 data->use_assembly_fallback = true; 2098 } 2099 } 2100 CeedCallBackend(CeedDestroy(&ceed)); 2101 CeedCallBackend(CeedQFunctionDestroy(&qf)); 2102 return CEED_ERROR_SUCCESS; 2103 } 2104 2105 extern "C" int CeedOperatorBuildKernelDiagonalAssemblyAtPoints_Cuda_gen(CeedOperator op, bool *is_good_build) { 2106 return CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(op, false, is_good_build); 2107 } 2108 2109 extern "C" int CeedOperatorBuildKernelFullAssemblyAtPoints_Cuda_gen(CeedOperator op, bool *is_good_build) { 2110 return CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(op, true, is_good_build); 2111 } 2112 2113 //------------------------------------------------------------------------------ 2114 // Build QFunction assembly operator kernel 2115 //------------------------------------------------------------------------------ 2116 extern "C" int CeedOperatorBuildKernelLinearAssembleQFunction_Cuda_gen(CeedOperator op, bool *is_good_build) { 2117 bool is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false; 2118 Ceed ceed; 2119 CeedInt Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0; 2120 CeedQFunctionField *qf_input_fields, *qf_output_fields; 2121 CeedQFunction_Cuda_gen *qf_data; 2122 CeedQFunction qf; 2123 CeedOperatorField *op_input_fields, *op_output_fields; 2124 CeedOperator_Cuda_gen *data; 2125 std::ostringstream code; 2126 Tab tab; 2127 2128 // Check compatibility 2129 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 2130 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); 2131 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "AtPoints QFunction assembly is not supported"); 2132 2133 // Check field compatibility 2134 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 2135 { 2136 bool has_shared_bases = true; 2137 2138 for (CeedInt i = 0; i < num_input_fields; i++) { 2139 CeedBasis basis; 2140 2141 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 2142 if (basis != CEED_BASIS_NONE) { 2143 bool is_tensor = true; 2144 const char *resource; 2145 char *resource_root; 2146 Ceed basis_ceed; 2147 2148 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 2149 is_all_tensor = is_all_tensor && is_tensor; 2150 is_all_nontensor = is_all_nontensor && !is_tensor; 2151 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); 2152 CeedCallBackend(CeedGetResource(basis_ceed, &resource)); 2153 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); 2154 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared"); 2155 CeedCallBackend(CeedFree(&resource_root)); 2156 CeedCallBackend(CeedDestroy(&basis_ceed)); 2157 } 2158 CeedCallBackend(CeedBasisDestroy(&basis)); 2159 } 2160 2161 for (CeedInt i = 0; i < num_output_fields; i++) { 2162 CeedBasis basis; 2163 2164 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 2165 if (basis != CEED_BASIS_NONE) { 2166 bool is_tensor = true; 2167 const char *resource; 2168 char *resource_root; 2169 Ceed basis_ceed; 2170 2171 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); 2172 is_all_tensor = is_all_tensor && is_tensor; 2173 is_all_nontensor = is_all_nontensor && !is_tensor; 2174 2175 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); 2176 CeedCallBackend(CeedGetResource(basis_ceed, &resource)); 2177 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); 2178 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared"); 2179 CeedCallBackend(CeedFree(&resource_root)); 2180 CeedCallBackend(CeedDestroy(&basis_ceed)); 2181 } 2182 CeedCallBackend(CeedBasisDestroy(&basis)); 2183 } 2184 } 2185 2186 // Retrieve operator data 2187 CeedCallBackend(CeedOperatorGetData(op, &data)); 2188 Q = data->Q; 2189 Q_1d = data->Q_1d; 2190 max_dim = data->dim; 2191 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 2192 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 2193 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 2194 2195 // Add atomicAdd function for old NVidia architectures 2196 { 2197 Ceed_Cuda *ceed_data; 2198 struct cudaDeviceProp prop; 2199 2200 CeedCallBackend(CeedGetData(ceed, &ceed_data)); 2201 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id)); 2202 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) { 2203 code << tab << "// AtomicAdd fallback source\n"; 2204 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n"; 2205 } 2206 } 2207 2208 // Load basis source files 2209 if (!is_all_nontensor) { 2210 code << tab << "// Tensor basis source\n"; 2211 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n"; 2212 } 2213 if (!is_all_tensor) { 2214 code << tab << "// Non-tensor basis source\n"; 2215 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-nontensor-templates.h>\n\n"; 2216 } 2217 if (!is_all_tensor && !is_all_nontensor) { 2218 code << "// Tensor basis source\n"; 2219 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-flattened-templates.h>\n\n"; 2220 } 2221 code << "// CodeGen operator source\n"; 2222 code << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n"; 2223 2224 // Get QFunction name 2225 std::string qfunction_name(qf_data->qfunction_name); 2226 std::string operator_name; 2227 2228 operator_name = "CeedKernelCudaGenQFunctionAssembly_" + qfunction_name; 2229 2230 // Define CEED_Q_VLA 2231 code << "\n" << tab << "#undef CEED_Q_VLA\n"; 2232 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) { 2233 code << tab << "#define CEED_Q_VLA 1\n\n"; 2234 } else { 2235 code << tab << "#define CEED_Q_VLA " << Q_1d << "\n\n"; 2236 } 2237 2238 // Add user QFunction source 2239 { 2240 const char *source_path; 2241 2242 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path)); 2243 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file"); 2244 2245 code << tab << "// User QFunction source\n"; 2246 code << tab << "#include \"" << source_path << "\"\n\n"; 2247 } 2248 2249 // Setup 2250 code << "\n" << tab << "// -----------------------------------------------------------------------------\n"; 2251 code << tab << "// Operator Assembly Kernel\n"; 2252 code << tab << "// \n"; 2253 code << tab << "// d_[in,out]_i: CeedVector device array\n"; 2254 code << tab << "// r_[in,out]_e_i: Element vector register\n"; 2255 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n"; 2256 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n"; 2257 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; 2258 code << tab << "// \n"; 2259 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; 2260 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; 2261 code << tab << "// -----------------------------------------------------------------------------\n"; 2262 code << tab << "extern \"C\" __global__ void " << operator_name 2263 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda " 2264 "points, CeedScalar *__restrict__ values_array) {\n"; 2265 tab.push(); 2266 2267 // Scratch buffers 2268 for (CeedInt i = 0; i < num_input_fields; i++) { 2269 CeedEvalMode eval_mode; 2270 2271 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 2272 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 2273 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n"; 2274 } 2275 } 2276 for (CeedInt i = 0; i < num_output_fields; i++) { 2277 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n"; 2278 } 2279 2280 code << tab << "const CeedInt max_dim = " << max_dim << ";\n"; 2281 if (!is_all_tensor) { 2282 code << tab << "const CeedInt Q = " << Q << ";\n"; 2283 } 2284 if (!is_all_nontensor) { 2285 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n"; 2286 } 2287 2288 // Shared data 2289 code << tab << "extern __shared__ CeedScalar slice[];\n"; 2290 code << tab << "SharedData_Cuda data;\n"; 2291 code << tab << "data.t_id_x = threadIdx.x;\n"; 2292 code << tab << "data.t_id_y = threadIdx.y;\n"; 2293 code << tab << "data.t_id_z = threadIdx.z;\n"; 2294 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 2295 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n"; 2296 2297 // -- Determine input mat reuse 2298 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX]; 2299 2300 for (CeedInt i = 0; i < num_input_fields; i++) { 2301 input_matrix_reuse[i].index = -1; 2302 } 2303 for (CeedInt i = 0; i < num_input_fields; i++) { 2304 bool is_tensor = true; 2305 CeedEvalMode eval_mode_i; 2306 CeedBasis basis_i; 2307 2308 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i)); 2309 if (eval_mode_i == CEED_EVAL_WEIGHT) continue; 2310 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i)); 2311 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor)); 2312 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) { 2313 CeedEvalMode eval_mode_j; 2314 CeedBasis basis_j; 2315 2316 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 2317 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 2318 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 2319 if (basis_i == basis_j) { 2320 if (is_tensor) { 2321 input_matrix_reuse[i].index = j; 2322 input_matrix_reuse[i].is_input = true; 2323 input_matrix_reuse[i].eval_mode = eval_mode_j; 2324 } else { 2325 // For non-tensor can only re-use with the same eval mode 2326 if (eval_mode_i == eval_mode_j) { 2327 input_matrix_reuse[i].index = j; 2328 input_matrix_reuse[i].is_input = true; 2329 input_matrix_reuse[i].eval_mode = eval_mode_j; 2330 } 2331 } 2332 } 2333 CeedCallBackend(CeedBasisDestroy(&basis_j)); 2334 } 2335 CeedCallBackend(CeedBasisDestroy(&basis_i)); 2336 } 2337 2338 // -- Determine output mat reuse 2339 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX]; 2340 2341 for (CeedInt i = 0; i < num_output_fields; i++) { 2342 output_matrix_reuse[i].index = -1; 2343 } 2344 for (CeedInt i = 0; i < num_output_fields; i++) { 2345 bool is_tensor = true; 2346 CeedEvalMode eval_mode_i; 2347 CeedBasis basis_i; 2348 2349 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i)); 2350 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i)); 2351 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor)); 2352 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) { 2353 CeedEvalMode eval_mode_j; 2354 CeedBasis basis_j; 2355 2356 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j)); 2357 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 2358 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j)); 2359 if (basis_i == basis_j) { 2360 if (is_tensor) { 2361 output_matrix_reuse[i].index = j; 2362 output_matrix_reuse[i].is_input = true; 2363 output_matrix_reuse[i].eval_mode = eval_mode_j; 2364 } else { 2365 // For non-tensor can only re-use with the same eval mode 2366 if (eval_mode_i == eval_mode_j) { 2367 output_matrix_reuse[i].index = j; 2368 output_matrix_reuse[i].is_input = true; 2369 output_matrix_reuse[i].eval_mode = eval_mode_j; 2370 } 2371 } 2372 } 2373 CeedCallBackend(CeedBasisDestroy(&basis_j)); 2374 } 2375 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) { 2376 CeedEvalMode eval_mode_j; 2377 CeedBasis basis_j; 2378 2379 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j)); 2380 if (eval_mode_j == CEED_EVAL_WEIGHT) continue; 2381 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j)); 2382 if (basis_i == basis_j) { 2383 if (is_tensor) { 2384 output_matrix_reuse[i].index = j; 2385 output_matrix_reuse[i].is_input = false; 2386 output_matrix_reuse[i].eval_mode = eval_mode_j; 2387 } else { 2388 // For non-tensor can only re-use with the same eval mode 2389 if (eval_mode_i == eval_mode_j) { 2390 output_matrix_reuse[i].index = j; 2391 output_matrix_reuse[i].is_input = false; 2392 output_matrix_reuse[i].eval_mode = eval_mode_j; 2393 } 2394 } 2395 } 2396 CeedCallBackend(CeedBasisDestroy(&basis_j)); 2397 } 2398 CeedCallBackend(CeedBasisDestroy(&basis_i)); 2399 } 2400 2401 // Initialize constants, and matrices B and G 2402 code << "\n" << tab << "// Input field constants and basis data\n"; 2403 for (CeedInt i = 0; i < num_input_fields; i++) { 2404 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i], 2405 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, true)); 2406 } 2407 code << "\n" << tab << "// Output field constants and basis data\n"; 2408 for (CeedInt i = 0; i < num_output_fields; i++) { 2409 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i], 2410 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, true)); 2411 } 2412 2413 // Loop over all elements 2414 code << "\n" << tab << "// Element loop\n"; 2415 code << tab << "__syncthreads();\n"; 2416 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; 2417 tab.push(); 2418 2419 // -- Compute minimum buffer space needed 2420 CeedInt max_rstr_buffer_size = 1; 2421 2422 for (CeedInt i = 0; i < num_input_fields; i++) { 2423 CeedEvalMode eval_mode; 2424 2425 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 2426 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) { 2427 CeedInt num_comp; 2428 CeedElemRestriction elem_rstr; 2429 2430 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 2431 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 2432 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 2433 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 2434 } 2435 } 2436 for (CeedInt i = 0; i < num_output_fields; i++) { 2437 CeedEvalMode eval_mode; 2438 2439 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 2440 if (eval_mode != CEED_EVAL_NONE) { 2441 CeedInt num_comp; 2442 CeedElemRestriction elem_rstr; 2443 2444 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 2445 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 2446 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1)); 2447 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 2448 } 2449 } 2450 code << tab << "// Scratch restriction buffer space\n"; 2451 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n"; 2452 2453 // -- Determine best input field processing order 2454 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX]; 2455 2456 for (CeedInt i = 0; i < num_input_fields; i++) { 2457 field_rstr_in_buffer[i] = -1; 2458 input_field_order[i] = -1; 2459 } 2460 { 2461 bool is_ordered[CEED_FIELD_MAX]; 2462 CeedInt curr_index = 0; 2463 2464 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false; 2465 for (CeedInt i = 0; i < num_input_fields; i++) { 2466 CeedVector vec_i; 2467 CeedElemRestriction rstr_i; 2468 2469 if (is_ordered[i]) continue; 2470 field_rstr_in_buffer[i] = i; 2471 is_ordered[i] = true; 2472 input_field_order[curr_index] = i; 2473 curr_index++; 2474 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i)); 2475 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT 2476 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i)); 2477 for (CeedInt j = i + 1; j < num_input_fields; j++) { 2478 CeedVector vec_j; 2479 CeedElemRestriction rstr_j; 2480 2481 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j)); 2482 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j)); 2483 if (rstr_i == rstr_j && vec_i == vec_j) { 2484 field_rstr_in_buffer[j] = i; 2485 is_ordered[j] = true; 2486 input_field_order[curr_index] = j; 2487 curr_index++; 2488 } 2489 CeedCallBackend(CeedVectorDestroy(&vec_j)); 2490 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j)); 2491 } 2492 CeedCallBackend(CeedVectorDestroy(&vec_i)); 2493 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i)); 2494 } 2495 } 2496 2497 // -- Input restriction and basis 2498 code << "\n" << tab << "// -- Input field restrictions and basis actions\n"; 2499 CeedInt num_active_in = 0, num_active_out = 0, qf_assembly_size_out = 0; 2500 CeedInt active_fields_in[CEED_FIELD_MAX], active_fields_out[CEED_FIELD_MAX]; 2501 2502 for (CeedInt i = 0; i < num_input_fields; i++) { 2503 bool is_active = false; 2504 const char *field_name; 2505 const CeedInt f = input_field_order[i]; 2506 2507 { 2508 CeedVector vec; 2509 2510 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec)); 2511 is_active = vec == CEED_VECTOR_ACTIVE; 2512 CeedCallBackend(CeedVectorDestroy(&vec)); 2513 } 2514 2515 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name)); 2516 code << tab << "// ---- Input field " << f << ": " << field_name << "\n"; 2517 2518 if (is_active) { 2519 CeedEvalMode eval_mode; 2520 CeedInt field_size; 2521 2522 active_fields_in[num_active_in] = f; 2523 num_active_in++; 2524 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[f], &field_size)); 2525 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[f], &eval_mode)); 2526 if (eval_mode == CEED_EVAL_GRAD) { 2527 code << tab << "CeedScalar r_q_in_" << f << "[num_comp_in_" << f << "*" << "dim_in_" << f << "*" 2528 << (is_all_tensor && (max_dim >= 3) ? "Q_1d" : "1") << "] = {0.};\n"; 2529 } else { 2530 code << tab << "CeedScalar r_q_in_" << f << "[num_comp_in_" << f << "*" << (is_all_tensor && (max_dim >= 3) ? "Q_1d" : "1") << "] = {0.};\n"; 2531 } 2532 code << tab << "const CeedInt field_size_in_" << f << " = " << field_size << ";\n"; 2533 } else { 2534 // ---- Restriction 2535 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], 2536 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices)); 2537 2538 // ---- Basis action 2539 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true, 2540 is_all_tensor, is_at_points, use_3d_slices)); 2541 } 2542 } 2543 code << tab << "const CeedInt field_sizes_in[" << num_active_in << "] = {"; 2544 for (CeedInt i = 0; i < num_active_in; i++) { 2545 code << "field_size_in_" << active_fields_in[i] << (i < num_active_in - 1 ? ", " : ""); 2546 } 2547 code << "};\n"; 2548 code << tab << "CeedScalar * r_q_in[" << num_active_in << "] = {"; 2549 for (CeedInt i = 0; i < num_active_in; i++) { 2550 code << "r_q_in_" << active_fields_in[i] << (i < num_active_in - 1 ? ", " : ""); 2551 } 2552 code << "};\n"; 2553 2554 for (CeedInt i = 0; i < num_output_fields; i++) { 2555 bool is_active = false; 2556 2557 { 2558 CeedVector vec; 2559 2560 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 2561 is_active = vec == CEED_VECTOR_ACTIVE; 2562 CeedCallBackend(CeedVectorDestroy(&vec)); 2563 } 2564 if (is_active) { 2565 const char *field_name; 2566 CeedInt field_size; 2567 2568 active_fields_out[num_active_out] = i; 2569 num_active_out++; 2570 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 2571 qf_assembly_size_out += field_size; 2572 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 2573 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 2574 code << tab << "const CeedInt field_size_out_" << i << " = " << field_size << ";\n"; 2575 } 2576 } 2577 code << tab << "const CeedInt field_sizes_out[" << num_active_out << "] = {"; 2578 for (CeedInt i = 0; i < num_active_out; i++) { 2579 code << "field_size_out_" << active_fields_out[i] << (i < num_active_out - 1 ? ", " : ""); 2580 } 2581 code << "};\n"; 2582 code << tab << "const CeedInt total_size_out = " << qf_assembly_size_out << ";\n"; 2583 2584 // -- Loop over active field 2585 code << "\n" << tab << "CeedInt input_offset = 0;\n"; 2586 code << tab << "// Loop over active QFunction input fields\n"; 2587 code << tab << "const CeedInt num_active_in = " << num_active_in << ";\n"; 2588 code << tab << "for (CeedInt a = 0; a < num_active_in; a++) {\n"; 2589 tab.push(); 2590 2591 // -- Loop over size of active field 2592 code << "\n" << tab << "// Loop over current active input field size\n"; 2593 code << tab << "const CeedInt field_size_in = field_sizes_in[a];\n"; 2594 code << tab << "for (CeedInt s = 0; s < field_size_in; s++) {\n"; 2595 tab.push(); 2596 2597 // -- Set current active point and component to 1 2598 code << tab << "// Set current active point and component to 1.0\n"; 2599 if (is_all_tensor && (max_dim >= 3)) { 2600 code << tab << "for (CeedInt i = 0; i < Q_1d; i++) r_q_in[a][i + s * Q_1d] = 1.0;\n"; 2601 } else { 2602 code << tab << "r_q_in[a][s] = 1.0;\n"; 2603 } 2604 2605 // -- Q function 2606 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields, 2607 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name, 2608 Q_1d, is_all_tensor, is_at_points, use_3d_slices)); 2609 2610 // -- Output basis and restriction 2611 code << "\n" << tab << "// -- Output field basis action and restrictions\n"; 2612 CeedScalar offset = 0; 2613 2614 for (CeedInt i = 0; i < num_output_fields; i++) { 2615 bool is_active = false; 2616 const char *field_name; 2617 2618 { 2619 CeedVector vec; 2620 2621 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 2622 is_active = vec == CEED_VECTOR_ACTIVE; 2623 CeedCallBackend(CeedVectorDestroy(&vec)); 2624 } 2625 if (!is_active) continue; 2626 2627 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name)); 2628 code << tab << "// ---- Output field " << i << ": " << field_name << "\n"; 2629 2630 // ---- Restriction 2631 CeedInt field_size; 2632 2633 code << tab << "WriteLVecStandard" << (is_all_tensor ? max_dim : 1) << "d_QFAssembly<total_size_out, field_size_out_" << i << ", " 2634 << (is_all_tensor ? "Q_1d" : "Q") << ">(data, num_elem, elem, input_offset + s, " << offset << ", r_q_out_" << i << ", values_array);\n"; 2635 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 2636 offset += field_size; 2637 } 2638 2639 // -- Reset current active node and component 2640 code << "\n" << tab << "// Reset current active node and component to 0.0\n"; 2641 if (is_all_tensor && (max_dim >= 3)) { 2642 code << tab << "for (CeedInt i = 0; i < Q_1d; i++) r_q_in[a][i + s * Q_1d] = 0.0;\n"; 2643 } else { 2644 code << tab << "r_q_in[a][s] = 0.0;\n"; 2645 } 2646 2647 // -- End of loop over size of active field 2648 tab.pop(); 2649 code << tab << "}\n"; 2650 code << tab << "input_offset += field_size_in;\n"; 2651 2652 // -- End of loop over active field 2653 tab.pop(); 2654 code << tab << "}\n"; 2655 2656 // Close loop and function 2657 tab.pop(); 2658 code << tab << "}\n"; 2659 tab.pop(); 2660 code << tab << "}\n"; 2661 code << tab << "// -----------------------------------------------------------------------------\n\n"; 2662 2663 // Compile 2664 { 2665 bool is_compile_good = false; 2666 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d); 2667 2668 data->thread_1d = T_1d; 2669 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good, &data->module_assemble_qfunction, 1, "OP_T_1D", T_1d)); 2670 if (is_compile_good) { 2671 *is_good_build = true; 2672 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module_assemble_qfunction, operator_name.c_str(), &data->assemble_qfunction)); 2673 } else { 2674 *is_good_build = false; 2675 data->use_assembly_fallback = true; 2676 } 2677 } 2678 CeedCallBackend(CeedDestroy(&ceed)); 2679 CeedCallBackend(CeedQFunctionDestroy(&qf)); 2680 return CEED_ERROR_SUCCESS; 2681 } 2682 2683 //------------------------------------------------------------------------------ 2684