1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #define CEED_DEBUG_COLOR 12 9 10 #include <ceed.h> 11 #include <ceed/backend.h> 12 #include <ceed/jit-tools.h> 13 14 #include <iostream> 15 #include <sstream> 16 #include <string> 17 18 #include "../hip-ref/ceed-hip-ref.h" 19 #include "../hip-shared/ceed-hip-shared.h" 20 #include "../hip/ceed-hip-common.h" 21 #include "../hip/ceed-hip-compile.h" 22 #include "ceed-hip-gen.h" 23 24 //------------------------------------------------------------------------------ 25 // Calculate the block size used for launching the operator kernel 26 //------------------------------------------------------------------------------ 27 extern "C" int BlockGridCalculate_Hip_gen(const CeedInt dim, const CeedInt num_elem, const CeedInt P_1d, const CeedInt Q_1d, CeedInt *block_sizes) { 28 const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 29 if (dim == 1) { 30 CeedInt elems_per_block = 64 * thread_1d > 256 ? 256 / thread_1d : 64; 31 32 elems_per_block = elems_per_block > 0 ? elems_per_block : 1; 33 block_sizes[0] = thread_1d; 34 block_sizes[1] = 1; 35 block_sizes[2] = elems_per_block; 36 } else if (dim == 2) { 37 const CeedInt elems_per_block = thread_1d < 4 ? 16 : 2; 38 39 block_sizes[0] = thread_1d; 40 block_sizes[1] = thread_1d; 41 block_sizes[2] = elems_per_block; 42 } else if (dim == 3) { 43 const CeedInt elems_per_block = thread_1d < 6 ? 4 : (thread_1d < 8 ? 2 : 1); 44 45 block_sizes[0] = thread_1d; 46 block_sizes[1] = thread_1d; 47 block_sizes[2] = elems_per_block; 48 } 49 return CEED_ERROR_SUCCESS; 50 } 51 52 //------------------------------------------------------------------------------ 53 // Determine type of operator 54 //------------------------------------------------------------------------------ 55 static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fields, CeedOperatorField *op_input_fields, 56 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, 57 CeedQFunctionField *qf_output_fields, CeedInt *max_P_1d, CeedInt *Q_1d, CeedInt *dim, bool *is_tensor, 58 bool *use_3d_slices) { 59 // Find dim, P_1d, Q_1d 60 *max_P_1d = 0; 61 *Q_1d = 0; 62 *dim = 0; 63 *is_tensor = true; 64 for (CeedInt i = 0; i < num_input_fields; i++) { 65 CeedBasis basis; 66 67 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 68 if (basis != CEED_BASIS_NONE) { 69 bool is_field_tensor; 70 CeedInt field_P_1d = 0, field_Q_1d = 0, field_dim = 0; 71 72 // Collect dim, P_1d, and Q_1d 73 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 74 CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 75 *is_tensor = *is_tensor && is_field_tensor; 76 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 77 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 78 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 79 CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 80 *dim = field_dim; 81 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 82 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 83 *Q_1d = field_Q_1d; 84 } 85 CeedCallBackend(CeedBasisDestroy(&basis)); 86 } 87 for (CeedInt i = 0; i < num_output_fields; i++) { 88 CeedBasis basis; 89 90 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 91 if (basis != CEED_BASIS_NONE) { 92 bool is_field_tensor; 93 CeedInt field_P_1d = 0, field_Q_1d = 0, field_dim = 0; 94 95 // Collect dim, P_1d, and Q_1d 96 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); 97 CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 98 *is_tensor = *is_tensor && is_field_tensor; 99 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); 100 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); 101 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); 102 CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 103 *dim = field_dim; 104 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); 105 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); 106 *Q_1d = field_Q_1d; 107 } 108 CeedCallBackend(CeedBasisDestroy(&basis)); 109 } 110 111 // Only use 3D collocated gradient parallelization strategy when gradient is computed 112 *use_3d_slices = false; 113 if (*dim == 3) { 114 bool was_grad_found = false; 115 116 for (CeedInt i = 0; i < num_input_fields; i++) { 117 CeedEvalMode eval_mode; 118 119 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 120 if (eval_mode == CEED_EVAL_GRAD) { 121 CeedBasis_Hip_shared *basis_data; 122 CeedBasis basis; 123 124 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 125 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 126 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 127 was_grad_found = true; 128 CeedCallBackend(CeedBasisDestroy(&basis)); 129 } 130 } 131 for (CeedInt i = 0; i < num_output_fields; i++) { 132 CeedEvalMode eval_mode; 133 134 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 135 if (eval_mode == CEED_EVAL_GRAD) { 136 CeedBasis_Hip_shared *basis_data; 137 CeedBasis basis; 138 139 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 140 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 141 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); 142 was_grad_found = true; 143 CeedCallBackend(CeedBasisDestroy(&basis)); 144 } 145 } 146 } 147 return CEED_ERROR_SUCCESS; 148 } 149 150 //------------------------------------------------------------------------------ 151 // Setup fields 152 //------------------------------------------------------------------------------ 153 static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedOperatorField op_field, 154 CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { 155 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 156 std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; 157 std::string option_name = (is_input ? "inputs" : "outputs"); 158 CeedEvalMode eval_mode = CEED_EVAL_NONE; 159 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 160 CeedElemRestriction elem_rstr; 161 CeedBasis_Hip_shared *basis_data; 162 CeedBasis basis; 163 164 code << " // -- " << (is_input ? "Input" : "Output") << " field " << i << "\n"; 165 166 // Get field data 167 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 168 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 169 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 170 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 171 } 172 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 173 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 174 if (basis != CEED_BASIS_NONE) { 175 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 176 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 177 } 178 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 179 180 // Set field constants 181 if (eval_mode != CEED_EVAL_WEIGHT) { 182 code << " const CeedInt " << P_name << " = " << (basis == CEED_BASIS_NONE ? Q_1d : P_1d) << ";\n"; 183 code << " const CeedInt num_comp" << var_suffix << " = " << num_comp << ";\n"; 184 } 185 186 // Load basis data 187 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 188 switch (eval_mode) { 189 case CEED_EVAL_NONE: 190 break; 191 case CEED_EVAL_INTERP: 192 if (is_at_points) { 193 // AtPoints 194 if (!basis_data->d_chebyshev_interp_1d) { 195 CeedSize interp_bytes; 196 CeedScalar *chebyshev_interp_1d; 197 198 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 199 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 200 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 201 CeedCallHip(CeedBasisReturnCeed(basis), hipMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); 202 CeedCallHip(CeedBasisReturnCeed(basis), 203 hipMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, hipMemcpyHostToDevice)); 204 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 205 } 206 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; 207 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; 208 } else { 209 // Standard quadrature 210 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 211 else data->B.outputs[i] = basis_data->d_interp_1d; 212 } 213 code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; 214 code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 215 break; 216 case CEED_EVAL_GRAD: 217 if (is_at_points) { 218 // AtPoints 219 if (!basis_data->d_chebyshev_interp_1d) { 220 CeedSize interp_bytes; 221 CeedScalar *chebyshev_interp_1d; 222 223 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); 224 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 225 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 226 CeedCallHip(CeedBasisReturnCeed(basis), hipMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); 227 CeedCallHip(CeedBasisReturnCeed(basis), 228 hipMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, hipMemcpyHostToDevice)); 229 CeedCallBackend(CeedFree(&chebyshev_interp_1d)); 230 } 231 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; 232 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; 233 } else { 234 // Standard quadrature 235 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; 236 else data->B.outputs[i] = basis_data->d_interp_1d; 237 } 238 code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; 239 code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; 240 if (is_at_points) break; // No G mat for AtPoints 241 if (use_3d_slices) { 242 if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; 243 else data->G.outputs[i] = basis_data->d_collo_grad_1d; 244 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; 245 code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 246 } else { 247 bool has_collo_grad = basis_data->d_collo_grad_1d; 248 249 if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 250 else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; 251 if (has_collo_grad) { 252 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; 253 code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 254 } else { 255 code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * P_1d << "];\n"; 256 code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; 257 } 258 } 259 break; 260 case CEED_EVAL_WEIGHT: 261 break; // No action 262 // LCOV_EXCL_START 263 case CEED_EVAL_DIV: 264 case CEED_EVAL_CURL: 265 break; // TODO: Not implemented 266 // LCOV_EXCL_STOP 267 } 268 CeedCallBackend(CeedBasisDestroy(&basis)); 269 return CEED_ERROR_SUCCESS; 270 } 271 272 //------------------------------------------------------------------------------ 273 // Restriction 274 //------------------------------------------------------------------------------ 275 static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, 276 CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field, 277 CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { 278 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 279 std::string P_name = "P_1d" + var_suffix; 280 CeedEvalMode eval_mode = CEED_EVAL_NONE; 281 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 282 CeedSize l_size; 283 CeedRestrictionType rstr_type = CEED_RESTRICTION_STANDARD; 284 CeedElemRestriction_Hip *rstr_data; 285 CeedElemRestriction elem_rstr; 286 CeedBasis basis; 287 288 // Get field data 289 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 290 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 291 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 292 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 293 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 294 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 295 } 296 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 297 if (basis != CEED_BASIS_NONE) { 298 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 299 } 300 CeedCallBackend(CeedBasisDestroy(&basis)); 301 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 302 303 // Restriction 304 if (is_input) { 305 // Input 306 if (field_input_buffer[i] != i) { 307 std::string buffer_name = "r_e_in_" + std::to_string(field_input_buffer[i]); 308 309 // Restriction was already done for previous input 310 code << " CeedScalar *r_e" << var_suffix << " = " << buffer_name << ";\n"; 311 } else if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices && is_at_points)) { 312 if (eval_mode == CEED_EVAL_NONE && rstr_type != CEED_RESTRICTION_POINTS) { 313 // No basis action, so r_e_in_* in also r_q_in_* and needs to be allocated 314 code << " CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; 315 } else if (rstr_type != CEED_RESTRICTION_POINTS) { 316 // Otherwise we're using the scratch space 317 code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 318 } 319 switch (rstr_type) { 320 case CEED_RESTRICTION_STANDARD: { 321 CeedInt comp_stride; 322 323 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 324 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 325 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 326 code << " // CompStride: " << comp_stride << "\n"; 327 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 328 code << " ReadLVecStandard" << dim << "d<num_comp" << var_suffix << ", " << comp_stride << ", " << P_name << ">(data, l_size" 329 << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; 330 break; 331 } 332 case CEED_RESTRICTION_STRIDED: { 333 bool has_backend_strides; 334 CeedInt num_elem; 335 336 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 337 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 338 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 339 340 if (!has_backend_strides) { 341 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 342 } 343 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 344 code << " ReadLVecStrided" << dim << "d<num_comp" << var_suffix << ", " << P_name << "," << strides[0] << "," << strides[1] << "," 345 << strides[2] << ">(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; 346 break; 347 } 348 case CEED_RESTRICTION_POINTS: { 349 CeedInt comp_stride; 350 351 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 352 code << " const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 353 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 354 break; 355 } 356 // LCOV_EXCL_START 357 case CEED_RESTRICTION_ORIENTED: 358 case CEED_RESTRICTION_CURL_ORIENTED: 359 break; // TODO: Not implemented 360 // LCOV_EXCL_STOP 361 } 362 } 363 } else { 364 // Output 365 switch (rstr_type) { 366 case CEED_RESTRICTION_STANDARD: { 367 CeedInt comp_stride; 368 369 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 370 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 371 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 372 code << " // CompStride: " << comp_stride << "\n"; 373 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; 374 code << " WriteLVecStandard" << dim << "d<num_comp" << var_suffix << ", " << comp_stride << ", " << P_name << ">(data, l_size" 375 << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; 376 break; 377 } 378 case CEED_RESTRICTION_STRIDED: { 379 bool has_backend_strides; 380 CeedInt num_elem; 381 382 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 383 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 384 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 385 386 if (!has_backend_strides) { 387 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 388 } 389 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 390 code << " WriteLVecStrided" << dim << "d<num_comp" << var_suffix << ", " << P_name << "," << strides[0] << "," << strides[1] << "," 391 << strides[2] << ">(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; 392 break; 393 } 394 case CEED_RESTRICTION_POINTS: 395 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; 396 break; 397 // LCOV_EXCL_START 398 case CEED_RESTRICTION_ORIENTED: 399 case CEED_RESTRICTION_CURL_ORIENTED: 400 break; // TODO: Not implemented 401 // LCOV_EXCL_STOP 402 } 403 } 404 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 405 return CEED_ERROR_SUCCESS; 406 } 407 408 //------------------------------------------------------------------------------ 409 // Basis 410 //------------------------------------------------------------------------------ 411 static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, 412 CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, 413 bool is_at_points, bool use_3d_slices) { 414 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); 415 std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; 416 CeedEvalMode eval_mode = CEED_EVAL_NONE; 417 CeedInt elem_size = 0, num_comp = 0, P_1d = 0; 418 CeedElemRestriction elem_rstr; 419 CeedBasis basis; 420 421 // Get field data 422 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr)); 423 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) { 424 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 425 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 426 } 427 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 428 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); 429 if (basis != CEED_BASIS_NONE) { 430 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); 431 } 432 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); 433 434 // Basis 435 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 436 if (is_input) { 437 switch (eval_mode) { 438 case CEED_EVAL_NONE: 439 if (!use_3d_slices && !is_at_points) { 440 code << " CeedScalar *r_q" << var_suffix << " = r_e" << var_suffix << ";\n"; 441 } 442 break; 443 case CEED_EVAL_INTERP: 444 if (is_at_points) { 445 code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; 446 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 447 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; 448 } else { 449 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 450 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 451 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 452 } 453 break; 454 case CEED_EVAL_GRAD: 455 if (is_at_points) { 456 code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; 457 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 458 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; 459 } else if (use_3d_slices) { 460 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 461 code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", P_1d" << var_suffix << ", " << Q_name 462 << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; 463 } else { 464 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; 465 code << " Grad" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d<num_comp" << var_suffix 466 << ", P_1d" << var_suffix << ", " << Q_name << ">(data, r_e" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" 467 << var_suffix << ");\n"; 468 } 469 break; 470 case CEED_EVAL_WEIGHT: { 471 if (is_at_points) { 472 code << " // Nothing to do AtPoints\n"; 473 } else { 474 CeedBasis_Hip_shared *basis_data; 475 476 code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; 477 CeedCallBackend(CeedBasisGetData(basis, &basis_data)); 478 data->W = basis_data->d_q_weight_1d; 479 code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; 480 } 481 break; 482 } 483 // LCOV_EXCL_START 484 case CEED_EVAL_DIV: 485 case CEED_EVAL_CURL: 486 break; // TODO: Not implemented 487 // LCOV_EXCL_STOP 488 } 489 } else { 490 switch (eval_mode) { 491 case CEED_EVAL_NONE: 492 code << " CeedScalar *r_e" << var_suffix << " = r_q" << var_suffix << ";\n"; 493 break; // No action 494 case CEED_EVAL_INTERP: 495 code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 496 if (is_at_points) { 497 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 498 << ">(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 499 } else { 500 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 501 << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 502 } 503 break; 504 case CEED_EVAL_GRAD: 505 code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; 506 if (is_at_points) { 507 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 508 << ">(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 509 } else if (use_3d_slices) { 510 code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name 511 << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; 512 } else { 513 code << " GradTranspose" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d<num_comp" 514 << var_suffix << ", " << P_name << "," << Q_name << ">(data, r_q" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix 515 << ", r_e" << var_suffix << ");\n"; 516 } 517 break; 518 // LCOV_EXCL_START 519 case CEED_EVAL_WEIGHT: 520 break; // Should not occur 521 case CEED_EVAL_DIV: 522 case CEED_EVAL_CURL: 523 break; // TODO: Not implemented 524 // LCOV_EXCL_STOP 525 } 526 } 527 CeedCallBackend(CeedBasisDestroy(&basis)); 528 return CEED_ERROR_SUCCESS; 529 } 530 531 //------------------------------------------------------------------------------ 532 // QFunction 533 //------------------------------------------------------------------------------ 534 static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt dim, CeedInt max_num_points, 535 CeedInt num_input_fields, CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, 536 CeedInt num_output_fields, CeedOperatorField *op_output_fields, 537 CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, bool is_at_points, 538 bool use_3d_slices) { 539 std::string Q_name = "Q_1d"; 540 CeedEvalMode eval_mode = CEED_EVAL_NONE; 541 CeedElemRestriction elem_rstr; 542 543 // Setup output arrays 544 code << "\n // -- Output field setup\n"; 545 for (CeedInt i = 0; i < num_output_fields; i++) { 546 std::string var_suffix = "_out_" + std::to_string(i); 547 548 code << " // ---- Output field " << i << "\n"; 549 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 550 switch (eval_mode) { 551 case CEED_EVAL_NONE: 552 if (is_at_points) { 553 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n"; 554 } else { 555 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 556 } 557 break; 558 case CEED_EVAL_INTERP: 559 if (is_at_points) { 560 // Accumulator for point data 561 code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; 562 code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; 563 code << " r_c" << var_suffix << "[i] = 0.0;\n"; 564 code << " }\n"; 565 } else { 566 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 567 } 568 break; 569 case CEED_EVAL_GRAD: 570 if (is_at_points) { 571 // Accumulator for point data 572 code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "*dim];\n"; 573 code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; 574 code << " r_c" << var_suffix << "[i] = 0.0;\n"; 575 code << " }\n"; 576 } else if (use_3d_slices) { 577 // Accumulator for gradient slices 578 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; 579 code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) {\n"; 580 code << " r_q" << var_suffix << "[i] = 0.0;\n"; 581 code << " }\n"; 582 } else { 583 code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; 584 } 585 break; 586 case CEED_EVAL_WEIGHT: 587 break; 588 // LCOV_EXCL_START 589 case CEED_EVAL_DIV: 590 case CEED_EVAL_CURL: 591 break; // TODO: Not implemented 592 // LCOV_EXCL_STOP 593 } 594 } 595 596 if (is_at_points) { 597 // We need to handle batches of points 598 code << "\n // Note: Using batches of points\n"; 599 code << " const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));\n\n"; 600 code << " #pragma unroll\n"; 601 code << " for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {\n"; 602 code << " const CeedInt p = i % max_num_points;\n\n"; 603 604 code << " // -- Coordinates\n"; 605 code << " CeedScalar r_x[dim];\n"; 606 code << " ReadPoint<dim, coords_comp_stride, max_num_points>(data, elem, p, max_num_points, points.indices, points.coords, r_x);\n\n"; 607 608 code << " // -- Input fields\n"; 609 for (CeedInt i = 0; i < num_input_fields; i++) { 610 std::string var_suffix = "_in_" + std::to_string(i); 611 612 code << " // ---- Input field " << i << "\n"; 613 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 614 // Basis action 615 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 616 switch (eval_mode) { 617 case CEED_EVAL_NONE: 618 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 619 code << " ReadPoint<num_comp" << var_suffix << ", comp_stride" << var_suffix 620 << ", max_num_points>(data, elem, p, max_num_points, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; 621 break; 622 case CEED_EVAL_INTERP: 623 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 624 code << " InterpAtPoints" << dim << "d<num_comp" << var_suffix << ", max_num_points, " << Q_name << ">(data, i, r_c" << var_suffix 625 << ", r_x, r_s" << var_suffix << ");\n"; 626 break; 627 case CEED_EVAL_GRAD: 628 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 629 code << " GradAtPoints" << dim << "d<num_comp" << var_suffix << ", max_num_points, " << Q_name << ">(data, i, r_c" << var_suffix 630 << ", r_x, r_s" << var_suffix << ");\n"; 631 break; 632 case CEED_EVAL_WEIGHT: 633 code << " CeedScalar r_s" << var_suffix << "[1];\n"; 634 code << " r_s" << var_suffix << "[0] = 1.0;\n"; 635 break; 636 // LCOV_EXCL_START 637 case CEED_EVAL_DIV: 638 case CEED_EVAL_CURL: 639 break; // TODO: Not implemented 640 // LCOV_EXCL_STOP 641 } 642 } 643 code << "\n // -- Output fields\n"; 644 for (CeedInt i = 0; i < num_output_fields; i++) { 645 std::string var_suffix = "_out_" + std::to_string(i); 646 647 code << " // ---- Output field " << i << "\n"; 648 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 649 // Basis action 650 switch (eval_mode) { 651 case CEED_EVAL_NONE: 652 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 653 break; 654 case CEED_EVAL_INTERP: 655 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 656 break; 657 case CEED_EVAL_GRAD: 658 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 659 break; 660 // LCOV_EXCL_START 661 case CEED_EVAL_WEIGHT: 662 break; // Should not occur 663 case CEED_EVAL_DIV: 664 case CEED_EVAL_CURL: 665 break; // TODO: Not implemented 666 // LCOV_EXCL_STOP 667 } 668 } 669 670 } else if (use_3d_slices) { 671 // We treat quadrature points per slice in 3d to save registers 672 code << "\n // Note: Using planes of 3D elements\n"; 673 code << " #pragma unroll\n"; 674 code << " for (CeedInt q = 0; q < " << Q_name << "; q++) {\n"; 675 code << " // -- Input fields\n"; 676 for (CeedInt i = 0; i < num_input_fields; i++) { 677 std::string var_suffix = "_in_" + std::to_string(i); 678 679 code << " // ---- Input field " << i << "\n"; 680 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 681 // Basis action 682 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 683 switch (eval_mode) { 684 case CEED_EVAL_NONE: 685 bool is_strided; 686 687 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 688 689 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 690 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 691 if (is_strided) { 692 bool has_backend_strides; 693 CeedInt num_elem, elem_size; 694 695 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 696 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides)); 697 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); 698 CeedInt strides[3] = {1, elem_size * num_elem, elem_size}; 699 700 if (!has_backend_strides) { 701 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); 702 } 703 code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; 704 code << " ReadEVecSliceStrided3d<num_comp" << var_suffix << ", " << Q_name << "," << strides[0] << "," << strides[1] << "," 705 << strides[2] << ">(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n"; 706 } else { 707 CeedSize l_size = 0; 708 CeedInt comp_stride; 709 CeedElemRestriction_Hip *rstr_data; 710 711 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size)); 712 code << " const CeedInt l_size" << var_suffix << " = " << l_size << ";\n"; 713 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 714 code << " // CompStride: " << comp_stride << "\n"; 715 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data)); 716 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; 717 code << " ReadEVecSliceStandard3d<num_comp" << var_suffix << ", " << comp_stride << ", " << Q_name << ">(data, l_size" << var_suffix 718 << ", elem, q, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; 719 } 720 break; 721 case CEED_EVAL_INTERP: 722 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 723 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) {\n"; 724 code << " r_s" << var_suffix << "[j] = r_q" << var_suffix << "[q + j*" << Q_name << "];\n"; 725 code << " }\n"; 726 break; 727 case CEED_EVAL_GRAD: 728 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 729 code << " GradColloSlice3d<num_comp" << var_suffix << ", " << Q_name << ">(data, q, r_q" << var_suffix << ", s_G" << var_suffix 730 << ", r_s" << var_suffix << ");\n"; 731 break; 732 case CEED_EVAL_WEIGHT: 733 code << " CeedScalar r_s" << var_suffix << "[1];\n"; 734 code << " r_s" << var_suffix << "[0] = r_q" << var_suffix << "[q];\n"; 735 break; 736 // LCOV_EXCL_START 737 case CEED_EVAL_DIV: 738 case CEED_EVAL_CURL: 739 break; // TODO: Not implemented 740 // LCOV_EXCL_STOP 741 } 742 } 743 code << "\n // -- Output fields\n"; 744 for (CeedInt i = 0; i < num_output_fields; i++) { 745 std::string var_suffix = "_out_" + std::to_string(i); 746 747 code << " // ---- Output field " << i << "\n"; 748 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 749 // Basis action 750 switch (eval_mode) { 751 case CEED_EVAL_NONE: 752 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 753 break; 754 case CEED_EVAL_INTERP: 755 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; 756 break; 757 case CEED_EVAL_GRAD: 758 code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; 759 break; 760 // LCOV_EXCL_START 761 case CEED_EVAL_WEIGHT: 762 break; // Should not occur 763 case CEED_EVAL_DIV: 764 case CEED_EVAL_CURL: 765 break; // TODO: Not implemented 766 // LCOV_EXCL_STOP 767 } 768 } 769 } else { 770 code << "\n // Note: Using full elements\n"; 771 code << " {\n"; 772 code << " // -- Input fields\n"; 773 for (CeedInt i = 0; i < num_input_fields; i++) { 774 code << " // ---- Input field " << i << "\n"; 775 code << " CeedScalar *r_s_in_" << i << " = r_q_in_" << i << ";\n"; 776 } 777 code << " // -- Output fields\n"; 778 for (CeedInt i = 0; i < num_output_fields; i++) { 779 code << " // ---- Output field " << i << "\n"; 780 code << " CeedScalar *r_s_out_" << i << " = r_q_out_" << i << ";\n"; 781 } 782 } 783 784 // Input and output buffers 785 code << "\n // -- QFunction inputs and outputs\n"; 786 code << " // ---- Inputs\n"; 787 code << " CeedScalar *inputs[" << CeedIntMax(num_input_fields, 1) << "];\n"; 788 for (CeedInt i = 0; i < num_input_fields; i++) { 789 code << " // ------ Input field " << i << "\n"; 790 code << " inputs[" << i << "] = r_s_in_" << i << ";\n"; 791 } 792 code << " // ---- Outputs\n"; 793 code << " CeedScalar *outputs[" << CeedIntMax(num_output_fields, 1) << "];\n"; 794 for (CeedInt i = 0; i < num_output_fields; i++) { 795 code << " // ------ Output field " << i << "\n"; 796 code << " outputs[" << i << "] = r_s_out_" << i << ";\n"; 797 } 798 799 // Apply QFunction 800 code << "\n // -- Apply QFunction\n"; 801 code << " " << qfunction_name << "(ctx, "; 802 if (dim != 3 || is_at_points || use_3d_slices) { 803 code << "1"; 804 } else { 805 code << "Q_1d"; 806 } 807 code << ", inputs, outputs);\n"; 808 809 if (is_at_points) { 810 // Map back to coefficients 811 code << "\n // -- Output fields\n"; 812 for (CeedInt i = 0; i < num_output_fields; i++) { 813 std::string var_suffix = "_out_" + std::to_string(i); 814 std::string P_name = "P_1d" + var_suffix; 815 816 code << " // ---- Output field " << i << "\n"; 817 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 818 // Basis action 819 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 820 switch (eval_mode) { 821 case CEED_EVAL_NONE: { 822 CeedInt comp_stride; 823 CeedElemRestriction elem_rstr; 824 825 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 826 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); 827 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 828 code << " const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; 829 code << " WritePoint<num_comp" << var_suffix << ", comp_stride" << var_suffix 830 << ", max_num_points>(data, elem, i, points.num_per_elem[elem], indices.outputs[" << i << "]" 831 << ", r_s" << var_suffix << ", d" << var_suffix << ");\n"; 832 break; 833 } 834 case CEED_EVAL_INTERP: 835 code << " if (i >= points.num_per_elem[elem]) {\n"; 836 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n"; 837 code << " }\n"; 838 code << " InterpTransposeAtPoints" << dim << "d<num_comp" << var_suffix << ", max_num_points, " << Q_name << ">(data, i, r_s" 839 << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; 840 break; 841 case CEED_EVAL_GRAD: 842 code << " if (i >= points.num_per_elem[elem]) {\n"; 843 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "*dim; j++) r_s" << var_suffix << "[j] = 0.0;\n"; 844 code << " }\n"; 845 code << " GradTransposeAtPoints" << dim << "d<num_comp" << var_suffix << ", max_num_points, " << Q_name << ">(data, i, r_s" 846 << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; 847 break; 848 // LCOV_EXCL_START 849 case CEED_EVAL_WEIGHT: 850 break; // Should not occur 851 case CEED_EVAL_DIV: 852 case CEED_EVAL_CURL: 853 break; // TODO: Not implemented 854 // LCOV_EXCL_STOP 855 } 856 } 857 } else if (use_3d_slices) { 858 // Copy or apply transpose grad, if needed 859 code << "\n // -- Output fields\n"; 860 for (CeedInt i = 0; i < num_output_fields; i++) { 861 std::string var_suffix = "_out_" + std::to_string(i); 862 std::string P_name = "P_1d" + var_suffix; 863 864 code << " // ---- Output field " << i << "\n"; 865 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 866 // Basis action 867 code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; 868 switch (eval_mode) { 869 case CEED_EVAL_NONE: 870 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 871 code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 872 code << " }\n"; 873 break; 874 case CEED_EVAL_INTERP: 875 code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; 876 code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; 877 code << " }\n"; 878 break; 879 case CEED_EVAL_GRAD: 880 code << " GradColloSliceTranspose3d<num_comp" << var_suffix << ", " << Q_name << ">(data, q, r_s" << var_suffix << ", s_G" 881 << var_suffix << ", r_q" << var_suffix << ");\n"; 882 break; 883 // LCOV_EXCL_START 884 case CEED_EVAL_WEIGHT: 885 break; // Should not occur 886 case CEED_EVAL_DIV: 887 case CEED_EVAL_CURL: 888 break; // TODO: Not implemented 889 // LCOV_EXCL_STOP 890 } 891 } 892 } 893 code << " }\n"; 894 return CEED_ERROR_SUCCESS; 895 } 896 897 //------------------------------------------------------------------------------ 898 // Build single operator kernel 899 //------------------------------------------------------------------------------ 900 extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { 901 bool is_tensor = true, is_at_points = false, use_3d_slices = false; 902 Ceed ceed; 903 CeedInt Q_1d, num_input_fields, num_output_fields, dim = 1, max_num_points = 0, coords_comp_stride = 0; 904 CeedQFunctionField *qf_input_fields, *qf_output_fields; 905 CeedQFunction_Hip_gen *qf_data; 906 CeedQFunction qf; 907 CeedOperatorField *op_input_fields, *op_output_fields; 908 CeedOperator_Hip_gen *data; 909 std::ostringstream code; 910 911 { 912 bool is_setup_done; 913 914 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 915 if (is_setup_done) return CEED_ERROR_SUCCESS; 916 } 917 918 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 919 CeedCallBackend(CeedOperatorGetData(op, &data)); 920 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 921 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 922 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 923 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 924 925 // Get operator data 926 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); 927 CeedCallBackend(CeedOperatorBuildKernelData_Hip_gen(ceed, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, op_output_fields, 928 qf_output_fields, &data->max_P_1d, &Q_1d, &dim, &is_tensor, &use_3d_slices)); 929 if (dim == 0) dim = 1; 930 data->dim = dim; 931 if (is_at_points) { 932 CeedElemRestriction_Hip *rstr_data; 933 CeedElemRestriction rstr_points = NULL; 934 935 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 936 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 937 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride)); 938 CeedCallBackend(CeedElemRestrictionGetData(rstr_points, &rstr_data)); 939 data->points.indices = (CeedInt *)rstr_data->d_offsets; 940 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 941 } 942 if (is_at_points) use_3d_slices = false; 943 if (Q_1d == 0) { 944 if (is_at_points) Q_1d = max_num_points; 945 else CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q_1d)); 946 } 947 data->Q_1d = Q_1d; 948 949 // Check for restriction only identity operator 950 { 951 bool is_identity_qf; 952 953 CeedCallBackend(CeedQFunctionIsIdentity(qf, &is_identity_qf)); 954 if (is_identity_qf) { 955 CeedEvalMode eval_mode_in, eval_mode_out; 956 957 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[0], &eval_mode_in)); 958 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[0], &eval_mode_out)); 959 CeedCheck(eval_mode_in != CEED_EVAL_NONE || eval_mode_out != CEED_EVAL_NONE, ceed, CEED_ERROR_BACKEND, 960 "Backend does not implement restriction only identity operators"); 961 } 962 } 963 964 // Load basis source files 965 // TODO: Add non-tensor, AtPoints 966 code << "// Tensor basis source\n"; 967 code << "#include <ceed/jit-source/hip/hip-shared-basis-tensor-templates.h>\n\n"; 968 code << "// AtPoints basis source\n"; 969 code << "#include <ceed/jit-source/hip/hip-shared-basis-tensor-at-points-templates.h>\n\n"; 970 code << "// CodeGen operator source\n"; 971 code << "#include <ceed/jit-source/hip/hip-gen-templates.h>\n\n"; 972 973 // Get QFunction name 974 std::string qfunction_name(qf_data->qfunction_name); 975 std::string operator_name; 976 977 operator_name = "CeedKernelHipGenOperator_" + qfunction_name; 978 979 // Define CEED_Q_VLA 980 code << "\n#undef CEED_Q_VLA\n"; 981 if (dim != 3 || is_at_points || use_3d_slices) { 982 code << "#define CEED_Q_VLA 1\n\n"; 983 } else { 984 code << "#define CEED_Q_VLA " << Q_1d << "\n\n"; 985 } 986 987 // Add user QFunction source 988 { 989 const char *source_path; 990 991 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path)); 992 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/hip/gen backend requires QFunction source code file"); 993 994 code << "// User QFunction source\n"; 995 code << "#include \"" << source_path << "\"\n\n"; 996 } 997 998 // Setup 999 code << "\n// -----------------------------------------------------------------------------\n"; 1000 code << "// Operator Kernel\n"; 1001 code << "// \n"; 1002 code << "// d_[in,out]_i: CeedVector device array\n"; 1003 code << "// r_[in,out]_e_i: Element vector register\n"; 1004 code << "// r_[in,out]_q_i: Quadrature space vector register\n"; 1005 code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficents register\n"; 1006 code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; 1007 code << "// \n"; 1008 code << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; 1009 code << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; 1010 code << "// -----------------------------------------------------------------------------\n"; 1011 code << "\nextern \"C\" __launch_bounds__(BLOCK_SIZE)\n"; 1012 code << "__global__ void " << operator_name 1013 << "(CeedInt num_elem, void* ctx, FieldsInt_Hip indices, Fields_Hip fields, Fields_Hip B, Fields_Hip G, CeedScalar* W, Points_Hip points) {\n"; 1014 1015 // Scratch buffers 1016 for (CeedInt i = 0; i < num_input_fields; i++) { 1017 CeedEvalMode eval_mode; 1018 1019 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1020 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 1021 code << " const CeedScalar *d_in_" << i << " = fields.inputs[" << i << "];\n"; 1022 } 1023 } 1024 for (CeedInt i = 0; i < num_output_fields; i++) { 1025 code << " CeedScalar *d_out_" << i << " = fields.outputs[" << i << "];\n"; 1026 } 1027 1028 code << " const CeedInt dim = " << dim << ";\n"; 1029 code << " const CeedInt Q_1d = " << Q_1d << ";\n"; 1030 if (is_at_points) { 1031 code << " const CeedInt max_num_points = " << max_num_points << ";\n"; 1032 code << " const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; 1033 } 1034 1035 // Shared data 1036 code << " extern __shared__ CeedScalar slice[];\n"; 1037 code << " SharedData_Hip data;\n"; 1038 code << " data.t_id_x = threadIdx.x;\n"; 1039 code << " data.t_id_y = threadIdx.y;\n"; 1040 code << " data.t_id_z = threadIdx.z;\n"; 1041 code << " data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 1042 code << " data.slice = slice + data.t_id_z*T_1D" << (dim > 1 ? "*T_1D" : "") << ";\n"; 1043 1044 // Initialize constants, and matrices B and G 1045 code << "\n // Input field constants and basis data\n"; 1046 for (CeedInt i = 0; i < num_input_fields; i++) { 1047 CeedCallBackend( 1048 CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_at_points, use_3d_slices)); 1049 } 1050 code << "\n // Output field constants and basis data\n"; 1051 for (CeedInt i = 0; i < num_output_fields; i++) { 1052 CeedCallBackend( 1053 CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); 1054 } 1055 1056 // Loop over all elements 1057 code << "\n // Element loop\n"; 1058 code << " __syncthreads();\n"; 1059 code << " for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; 1060 1061 // -- Compute minimum buffer space needed 1062 CeedInt max_rstr_buffer_size = 1; 1063 1064 for (CeedInt i = 0; i < num_input_fields; i++) { 1065 CeedInt num_comp, elem_size; 1066 CeedElemRestriction elem_rstr; 1067 1068 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1069 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1070 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1071 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); 1072 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1073 } 1074 for (CeedInt i = 0; i < num_output_fields; i++) { 1075 CeedInt num_comp, elem_size; 1076 CeedElemRestriction elem_rstr; 1077 1078 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1079 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1080 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1081 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); 1082 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); 1083 } 1084 code << " // Scratch restriction buffer space\n"; 1085 code << " CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n"; 1086 1087 // -- Determine best input field processing order 1088 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX]; 1089 1090 for (CeedInt i = 0; i < num_input_fields; i++) { 1091 field_rstr_in_buffer[i] = -1; 1092 input_field_order[i] = -1; 1093 } 1094 { 1095 bool is_ordered[CEED_FIELD_MAX]; 1096 CeedInt curr_index = 0; 1097 1098 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false; 1099 for (CeedInt i = 0; i < num_input_fields; i++) { 1100 CeedVector vec_i; 1101 CeedElemRestriction rstr_i; 1102 1103 if (is_ordered[i]) continue; 1104 field_rstr_in_buffer[i] = i; 1105 is_ordered[i] = true; 1106 input_field_order[curr_index] = i; 1107 curr_index++; 1108 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i)); 1109 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT 1110 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i)); 1111 for (CeedInt j = i + 1; j < num_input_fields; j++) { 1112 CeedVector vec_j; 1113 CeedElemRestriction rstr_j; 1114 1115 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j)); 1116 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j)); 1117 if (rstr_i == rstr_j && vec_i == vec_j) { 1118 field_rstr_in_buffer[j] = i; 1119 is_ordered[j] = true; 1120 input_field_order[curr_index] = j; 1121 curr_index++; 1122 } 1123 CeedCallBackend(CeedVectorDestroy(&vec_j)); 1124 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j)); 1125 } 1126 CeedCallBackend(CeedVectorDestroy(&vec_i)); 1127 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i)); 1128 } 1129 } 1130 1131 // -- Input restriction and basis 1132 code << "\n // -- Input field restrictions and basis actions\n"; 1133 for (CeedInt i = 0; i < num_input_fields; i++) { 1134 CeedInt f = input_field_order[i]; 1135 1136 code << " // ---- Input field " << f << "\n"; 1137 1138 // ---- Restriction 1139 CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, f, dim, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], Q_1d, 1140 true, is_at_points, use_3d_slices)); 1141 1142 // ---- Basis action 1143 CeedCallBackend( 1144 CeedOperatorBuildKernelBasis_Hip_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_at_points, use_3d_slices)); 1145 } 1146 1147 // -- Q function 1148 CeedCallBackend(CeedOperatorBuildKernelQFunction_Hip_gen(code, data, dim, max_num_points, num_input_fields, op_input_fields, qf_input_fields, 1149 num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_at_points, 1150 use_3d_slices)); 1151 1152 // -- Output basis and restriction 1153 code << "\n // -- Output field basis action and restrictions\n"; 1154 for (CeedInt i = 0; i < num_output_fields; i++) { 1155 code << " // ---- Output field " << i << "\n"; 1156 1157 // ---- Basis action 1158 CeedCallBackend( 1159 CeedOperatorBuildKernelBasis_Hip_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); 1160 1161 // ---- Restriction 1162 CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, i, dim, NULL, op_output_fields[i], qf_output_fields[i], Q_1d, false, 1163 is_at_points, use_3d_slices)); 1164 } 1165 1166 // Close loop and function 1167 code << " }\n"; 1168 code << "}\n"; 1169 code << "// -----------------------------------------------------------------------------\n\n"; 1170 1171 CeedInt block_sizes[3] = {0, 0, 0}; 1172 CeedInt num_elem; 1173 1174 // Compile 1175 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1176 CeedCallBackend(BlockGridCalculate_Hip_gen(dim, num_elem, data->max_P_1d, Q_1d, block_sizes)); 1177 CeedCallBackend(CeedCompile_Hip(ceed, code.str().c_str(), &data->module, 2, "T_1D", block_sizes[0], "BLOCK_SIZE", 1178 block_sizes[0] * block_sizes[1] * block_sizes[2])); 1179 CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, operator_name.c_str(), &data->op)); 1180 CeedCallBackend(CeedOperatorSetSetupDone(op)); 1181 CeedCallBackend(CeedDestroy(&ceed)); 1182 CeedCallBackend(CeedQFunctionDestroy(&qf)); 1183 return CEED_ERROR_SUCCESS; 1184 } 1185 1186 //------------------------------------------------------------------------------ 1187