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