From 9f6cf23f05437e07130e84cd8d6d707fab1db3c2 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Tue, 17 Dec 2024 12:36:46 -0700 Subject: [PATCH] hip - AtPoints for hip/gen --- .../cuda-gen/ceed-cuda-gen-operator-build.cpp | 2 +- .../hip-gen/ceed-hip-gen-operator-build.cpp | 422 ++++++++++++++---- backends/hip-gen/ceed-hip-gen-operator.c | 49 ++ backends/hip-gen/ceed-hip-gen.c | 1 + backends/hip-gen/ceed-hip-gen.h | 1 + .../ceed/jit-source/hip/hip-gen-templates.h | 32 ++ include/ceed/jit-source/hip/hip-types.h | 7 + 7 files changed, 430 insertions(+), 84 deletions(-) diff --git a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp index cafcbd6d94..99da342f66 100644 --- a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp +++ b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp @@ -1156,7 +1156,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { code << "}\n"; code << "// -----------------------------------------------------------------------------\n\n"; - // View kernel for debugging + // Compile CeedCallBackend(CeedCompile_Cuda(ceed, code.str().c_str(), &data->module, 1, "T_1D", CeedIntMax(Q_1d, data->max_P_1d))); CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, operator_name.c_str(), &data->op)); CeedCallBackend(CeedOperatorSetSetupDone(op)); diff --git a/backends/hip-gen/ceed-hip-gen-operator-build.cpp b/backends/hip-gen/ceed-hip-gen-operator-build.cpp index 9922cb3a69..bd56a77350 100644 --- a/backends/hip-gen/ceed-hip-gen-operator-build.cpp +++ b/backends/hip-gen/ceed-hip-gen-operator-build.cpp @@ -25,25 +25,25 @@ // Calculate the block size used for launching the operator kernel //------------------------------------------------------------------------------ extern "C" int BlockGridCalculate_Hip_gen(const CeedInt dim, const CeedInt num_elem, const CeedInt P_1d, const CeedInt Q_1d, CeedInt *block_sizes) { - const CeedInt thread1d = CeedIntMax(Q_1d, P_1d); + const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); if (dim == 1) { - CeedInt elems_per_block = 64 * thread1d > 256 ? 256 / thread1d : 64; + CeedInt elems_per_block = 64 * thread_1d > 256 ? 256 / thread_1d : 64; elems_per_block = elems_per_block > 0 ? elems_per_block : 1; - block_sizes[0] = thread1d; + block_sizes[0] = thread_1d; block_sizes[1] = 1; block_sizes[2] = elems_per_block; } else if (dim == 2) { - const CeedInt elems_per_block = thread1d < 4 ? 16 : 2; + const CeedInt elems_per_block = thread_1d < 4 ? 16 : 2; - block_sizes[0] = thread1d; - block_sizes[1] = thread1d; + block_sizes[0] = thread_1d; + block_sizes[1] = thread_1d; block_sizes[2] = elems_per_block; } else if (dim == 3) { - const CeedInt elems_per_block = thread1d < 6 ? 4 : (thread1d < 8 ? 2 : 1); + const CeedInt elems_per_block = thread_1d < 6 ? 4 : (thread_1d < 8 ? 2 : 1); - block_sizes[0] = thread1d; - block_sizes[1] = thread1d; + block_sizes[0] = thread_1d; + block_sizes[1] = thread_1d; block_sizes[2] = elems_per_block; } return CEED_ERROR_SUCCESS; @@ -82,6 +82,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } + CeedCallBackend(CeedBasisDestroy(&basis)); } for (CeedInt i = 0; i < num_output_fields; i++) { CeedBasis basis; @@ -104,6 +105,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } + CeedCallBackend(CeedBasisDestroy(&basis)); } // Only use 3D collocated gradient parallelization strategy when gradient is computed @@ -123,6 +125,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel CeedCallBackend(CeedBasisGetData(basis, &basis_data)); *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); was_grad_found = true; + CeedCallBackend(CeedBasisDestroy(&basis)); } } for (CeedInt i = 0; i < num_output_fields; i++) { @@ -137,6 +140,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel CeedCallBackend(CeedBasisGetData(basis, &basis_data)); *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true); was_grad_found = true; + CeedCallBackend(CeedBasisDestroy(&basis)); } } } @@ -147,7 +151,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel // Setup fields //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedOperatorField op_field, - CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool use_3d_slices) { + CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; std::string option_name = (is_input ? "inputs" : "outputs"); @@ -165,6 +169,7 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); } + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { CeedCallBackend(CeedBasisGetData(basis, &basis_data)); @@ -184,16 +189,55 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce case CEED_EVAL_NONE: break; case CEED_EVAL_INTERP: - if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; - else data->B.outputs[i] = basis_data->d_interp_1d; + if (is_at_points) { + // AtPoints + if (!basis_data->d_chebyshev_interp_1d) { + CeedSize interp_bytes; + CeedScalar *chebyshev_interp_1d; + + interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); + CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); + CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); + CeedCallHip(CeedBasisReturnCeed(basis), hipMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); + CeedCallHip(CeedBasisReturnCeed(basis), + hipMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, hipMemcpyHostToDevice)); + CeedCallBackend(CeedFree(&chebyshev_interp_1d)); + } + if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; + else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; + } else { + // Standard quadrature + if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; + else data->B.outputs[i] = basis_data->d_interp_1d; + } code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; break; case CEED_EVAL_GRAD: - if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; - else data->B.outputs[i] = basis_data->d_interp_1d; + if (is_at_points) { + // AtPoints + if (!basis_data->d_chebyshev_interp_1d) { + CeedSize interp_bytes; + CeedScalar *chebyshev_interp_1d; + + interp_bytes = P_1d * Q_1d * sizeof(CeedScalar); + CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); + CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); + CeedCallHip(CeedBasisReturnCeed(basis), hipMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes)); + CeedCallHip(CeedBasisReturnCeed(basis), + hipMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, hipMemcpyHostToDevice)); + CeedCallBackend(CeedFree(&chebyshev_interp_1d)); + } + if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d; + else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d; + } else { + // Standard quadrature + if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; + else data->B.outputs[i] = basis_data->d_interp_1d; + } code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; + if (is_at_points) break; // No G mat for AtPoints if (use_3d_slices) { if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; else data->G.outputs[i] = basis_data->d_collo_grad_1d; @@ -217,11 +261,11 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce break; // No action // LCOV_EXCL_START case CEED_EVAL_DIV: - break; // TODO: Not implemented case CEED_EVAL_CURL: break; // TODO: Not implemented // LCOV_EXCL_STOP } + CeedCallBackend(CeedBasisDestroy(&basis)); return CEED_ERROR_SUCCESS; } @@ -230,7 +274,7 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field, - CeedInt Q_1d, bool is_input, bool use_3d_slices) { + CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); std::string P_name = "P_1d" + var_suffix; CeedEvalMode eval_mode = CEED_EVAL_NONE; @@ -253,6 +297,7 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, if (basis != CEED_BASIS_NONE) { CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); } + CeedCallBackend(CeedBasisDestroy(&basis)); CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); // Restriction @@ -263,11 +308,11 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, // Restriction was already done for previous input code << " CeedScalar *r_e" << var_suffix << " = " << buffer_name << ";\n"; - } else if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices)) { - if (eval_mode == CEED_EVAL_NONE) { + } else if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices && is_at_points)) { + if (eval_mode == CEED_EVAL_NONE && rstr_type != CEED_RESTRICTION_POINTS) { // No basis action, so r_e_in_* in also r_q_in_* and needs to be allocated code << " CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n"; - } else { + } else if (rstr_type != CEED_RESTRICTION_POINTS) { // Otherwise we're using the scratch space code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; } @@ -300,10 +345,17 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, << strides[2] << ">(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; break; } + case CEED_RESTRICTION_POINTS: { + CeedInt comp_stride; + + CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); + code << " const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; + data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; + break; + } // LCOV_EXCL_START case CEED_RESTRICTION_ORIENTED: case CEED_RESTRICTION_CURL_ORIENTED: - case CEED_RESTRICTION_POINTS: break; // TODO: Not implemented // LCOV_EXCL_STOP } @@ -339,14 +391,17 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, << strides[2] << ">(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; break; } + case CEED_RESTRICTION_POINTS: + data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; + break; // LCOV_EXCL_START case CEED_RESTRICTION_ORIENTED: case CEED_RESTRICTION_CURL_ORIENTED: - case CEED_RESTRICTION_POINTS: break; // TODO: Not implemented // LCOV_EXCL_STOP } } + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); return CEED_ERROR_SUCCESS; } @@ -355,7 +410,7 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, - bool use_3d_slices) { + bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; CeedEvalMode eval_mode = CEED_EVAL_NONE; @@ -369,6 +424,7 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); } + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); @@ -380,17 +436,27 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp if (is_input) { switch (eval_mode) { case CEED_EVAL_NONE: - if (!use_3d_slices) { + if (!use_3d_slices && !is_at_points) { code << " CeedScalar *r_q" << var_suffix << " = r_e" << var_suffix << ";\n"; } break; case CEED_EVAL_INTERP: - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + if (is_at_points) { + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; + code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + } else { + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + } break; case CEED_EVAL_GRAD: - if (use_3d_slices) { + if (is_at_points) { + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; + code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + } else if (use_3d_slices) { code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; @@ -402,12 +468,16 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp } break; case CEED_EVAL_WEIGHT: { - CeedBasis_Hip_shared *basis_data; + if (is_at_points) { + code << " // Nothing to do AtPoints\n"; + } else { + CeedBasis_Hip_shared *basis_data; - code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; - CeedCallBackend(CeedBasisGetData(basis, &basis_data)); - data->W = basis_data->d_q_weight_1d; - code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; + code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; + CeedCallBackend(CeedBasisGetData(basis, &basis_data)); + data->W = basis_data->d_q_weight_1d; + code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; + } break; } // LCOV_EXCL_START @@ -423,12 +493,20 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp break; // No action case CEED_EVAL_INTERP: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + if (is_at_points) { + code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + } else { + code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + } break; case CEED_EVAL_GRAD: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; - if (use_3d_slices) { + if (is_at_points) { + code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + } else if (use_3d_slices) { code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; } else { @@ -446,16 +524,17 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp // LCOV_EXCL_STOP } } + CeedCallBackend(CeedBasisDestroy(&basis)); return CEED_ERROR_SUCCESS; } //------------------------------------------------------------------------------ // QFunction //------------------------------------------------------------------------------ -static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt dim, CeedInt num_input_fields, - CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, +static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt dim, CeedInt max_num_points, + CeedInt num_input_fields, CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, - CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, + CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, bool is_at_points, bool use_3d_slices) { std::string Q_name = "Q_1d"; CeedEvalMode eval_mode = CEED_EVAL_NONE; @@ -468,24 +547,128 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce code << " // ---- Output field " << i << "\n"; CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); - if (eval_mode == CEED_EVAL_NONE || eval_mode == CEED_EVAL_INTERP) { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + switch (eval_mode) { + case CEED_EVAL_NONE: + if (is_at_points) { + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n"; + } else { + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + } + break; + case CEED_EVAL_INTERP: + if (is_at_points) { + // Accumulator for point data + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " r_c" << var_suffix << "[i] = 0.0;\n"; + code << " }\n"; + } else { + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + } + break; + case CEED_EVAL_GRAD: + if (is_at_points) { + // Accumulator for point data + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "*dim];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " r_c" << var_suffix << "[i] = 0.0;\n"; + code << " }\n"; + } else if (use_3d_slices) { + // Accumulator for gradient slices + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) {\n"; + code << " r_q" << var_suffix << "[i] = 0.0;\n"; + code << " }\n"; + } else { + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; + } + break; + case CEED_EVAL_WEIGHT: + break; + // LCOV_EXCL_START + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + break; // TODO: Not implemented + // LCOV_EXCL_STOP } - if (eval_mode == CEED_EVAL_GRAD) { - if (use_3d_slices) { - // Accumulator for gradient slices - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) {\n"; - code << " r_q" << var_suffix << "[i] = 0.0;\n"; - code << " }\n"; - } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; + } + + if (is_at_points) { + // We need to handle batches of points + code << "\n // Note: Using batches of points\n"; + code << " const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));\n\n"; + code << " #pragma unroll\n"; + code << " for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {\n"; + code << " const CeedInt p = i % max_num_points;\n\n"; + + code << " // -- Coordinates\n"; + code << " CeedScalar r_x[dim];\n"; + code << " ReadPoint(data, elem, p, max_num_points, points.indices, points.coords, r_x);\n\n"; + + code << " // -- Input fields\n"; + for (CeedInt i = 0; i < num_input_fields; i++) { + std::string var_suffix = "_in_" + std::to_string(i); + + code << " // ---- Input field " << i << "\n"; + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + // Basis action + code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; + switch (eval_mode) { + case CEED_EVAL_NONE: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; + code << " ReadPoint(data, elem, p, max_num_points, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; + break; + case CEED_EVAL_INTERP: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; + code << " InterpAtPoints" << dim << "d(data, i, r_c" << var_suffix + << ", r_x, r_s" << var_suffix << ");\n"; + break; + case CEED_EVAL_GRAD: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; + code << " GradAtPoints" << dim << "d(data, i, r_c" << var_suffix + << ", r_x, r_s" << var_suffix << ");\n"; + break; + case CEED_EVAL_WEIGHT: + code << " CeedScalar r_s" << var_suffix << "[1];\n"; + code << " r_s" << var_suffix << "[0] = 1.0;\n"; + break; + // LCOV_EXCL_START + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + break; // TODO: Not implemented + // LCOV_EXCL_STOP } } - } + code << "\n // -- Output fields\n"; + for (CeedInt i = 0; i < num_output_fields; i++) { + std::string var_suffix = "_out_" + std::to_string(i); - // We treat quadrature points per slice in 3d to save registers - if (use_3d_slices) { + code << " // ---- Output field " << i << "\n"; + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); + // Basis action + switch (eval_mode) { + case CEED_EVAL_NONE: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; + break; + case CEED_EVAL_INTERP: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; + break; + case CEED_EVAL_GRAD: + code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; + break; + // LCOV_EXCL_START + case CEED_EVAL_WEIGHT: + break; // Should not occur + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + break; // TODO: Not implemented + // LCOV_EXCL_STOP + } + } + + } else if (use_3d_slices) { + // We treat quadrature points per slice in 3d to save registers code << "\n // Note: Using planes of 3D elements\n"; code << " #pragma unroll\n"; code << " for (CeedInt q = 0; q < " << Q_name << "; q++) {\n"; @@ -549,10 +732,9 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce case CEED_EVAL_WEIGHT: code << " CeedScalar r_s" << var_suffix << "[1];\n"; code << " r_s" << var_suffix << "[0] = r_q" << var_suffix << "[q];\n"; - break; // No action - // LCOV_EXCL_START + break; + // LCOV_EXCL_START case CEED_EVAL_DIV: - break; // TODO: Not implemented case CEED_EVAL_CURL: break; // TODO: Not implemented // LCOV_EXCL_STOP @@ -568,7 +750,7 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce switch (eval_mode) { case CEED_EVAL_NONE: code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; - break; // No action + break; case CEED_EVAL_INTERP: code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; break; @@ -579,7 +761,6 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce case CEED_EVAL_WEIGHT: break; // Should not occur case CEED_EVAL_DIV: - break; // TODO: Not implemented case CEED_EVAL_CURL: break; // TODO: Not implemented // LCOV_EXCL_STOP @@ -618,16 +799,64 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce // Apply QFunction code << "\n // -- Apply QFunction\n"; code << " " << qfunction_name << "(ctx, "; - if (dim != 3 || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices) { code << "1"; } else { code << "Q_1d"; } code << ", inputs, outputs);\n"; - // Copy or apply transpose grad, if needed - if (use_3d_slices) { - code << " // -- Output fields\n"; + if (is_at_points) { + // Map back to coefficients + code << "\n // -- Output fields\n"; + for (CeedInt i = 0; i < num_output_fields; i++) { + std::string var_suffix = "_out_" + std::to_string(i); + std::string P_name = "P_1d" + var_suffix; + + code << " // ---- Output field " << i << "\n"; + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); + // Basis action + code << " // EvalMode: " << CeedEvalModes[eval_mode] << "\n"; + switch (eval_mode) { + case CEED_EVAL_NONE: { + CeedInt comp_stride; + CeedElemRestriction elem_rstr; + + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); + code << " const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n"; + code << " WritePoint(data, elem, i, points.num_per_elem[elem], indices.outputs[" << i << "]" + << ", r_s" << var_suffix << ", d" << var_suffix << ");\n"; + break; + } + case CEED_EVAL_INTERP: + code << " if (i >= points.num_per_elem[elem]) {\n"; + code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n"; + code << " }\n"; + code << " InterpTransposeAtPoints" << dim << "d(data, i, r_s" + << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; + break; + case CEED_EVAL_GRAD: + code << " if (i >= points.num_per_elem[elem]) {\n"; + code << " for (CeedInt j = 0; j < num_comp" << var_suffix << "*dim; j++) r_s" << var_suffix << "[j] = 0.0;\n"; + code << " }\n"; + code << " GradTransposeAtPoints" << dim << "d(data, i, r_s" + << var_suffix << ", r_x, r_c" << var_suffix << ");\n"; + break; + // LCOV_EXCL_START + case CEED_EVAL_WEIGHT: + break; // Should not occur + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + break; // TODO: Not implemented + // LCOV_EXCL_STOP + } + } + } else if (use_3d_slices) { + // Copy or apply transpose grad, if needed + code << "\n // -- Output fields\n"; for (CeedInt i = 0; i < num_output_fields; i++) { std::string var_suffix = "_out_" + std::to_string(i); std::string P_name = "P_1d" + var_suffix; @@ -641,7 +870,7 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; code << " }\n"; - break; // No action + break; case CEED_EVAL_INTERP: code << " for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n"; code << " r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n"; @@ -655,7 +884,6 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce case CEED_EVAL_WEIGHT: break; // Should not occur case CEED_EVAL_DIV: - break; // TODO: Not implemented case CEED_EVAL_CURL: break; // TODO: Not implemented // LCOV_EXCL_STOP @@ -670,9 +898,9 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce // Build single operator kernel //------------------------------------------------------------------------------ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { - bool is_tensor = true, use_3d_slices = false; + bool is_tensor = true, is_at_points = false, use_3d_slices = false; Ceed ceed; - CeedInt Q_1d, num_input_fields, num_output_fields, dim = 1; + CeedInt Q_1d, num_input_fields, num_output_fields, dim = 1, max_num_points = 0, coords_comp_stride = 0; CeedQFunctionField *qf_input_fields, *qf_output_fields; CeedQFunction_Hip_gen *qf_data; CeedQFunction qf; @@ -695,15 +923,26 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); // Get operator data + CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); CeedCallBackend(CeedOperatorBuildKernelData_Hip_gen(ceed, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, &data->max_P_1d, &Q_1d, &dim, &is_tensor, &use_3d_slices)); if (dim == 0) dim = 1; data->dim = dim; + if (is_at_points) { + CeedElemRestriction_Hip *rstr_data; + CeedElemRestriction rstr_points = NULL; + + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); + CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); + CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride)); + CeedCallBackend(CeedElemRestrictionGetData(rstr_points, &rstr_data)); + data->points.indices = (CeedInt *)rstr_data->d_offsets; + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); + } + if (is_at_points) use_3d_slices = false; if (Q_1d == 0) { - CeedInt Q; - - CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); - Q_1d = Q; + if (is_at_points) Q_1d = max_num_points; + else CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q_1d)); } data->Q_1d = Q_1d; @@ -726,6 +965,8 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // TODO: Add non-tensor, AtPoints code << "// Tensor basis source\n"; code << "#include \n\n"; + code << "// AtPoints basis source\n"; + code << "#include \n\n"; code << "// CodeGen operator source\n"; code << "#include \n\n"; @@ -737,7 +978,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // Define CEED_Q_VLA code << "\n#undef CEED_Q_VLA\n"; - if (dim != 3 || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices) { code << "#define CEED_Q_VLA 1\n\n"; } else { code << "#define CEED_Q_VLA " << Q_1d << "\n\n"; @@ -761,14 +1002,15 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << "// d_[in,out]_i: CeedVector device array\n"; code << "// r_[in,out]_e_i: Element vector register\n"; code << "// r_[in,out]_q_i: Quadrature space vector register\n"; - code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; + code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficents register\n"; + code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; code << "// \n"; code << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; code << "// s_G_[in,out]_i: Gradient matrix, shared memory\n"; code << "// -----------------------------------------------------------------------------\n"; code << "\nextern \"C\" __launch_bounds__(BLOCK_SIZE)\n"; code << "__global__ void " << operator_name - << "(CeedInt num_elem, void* ctx, FieldsInt_Hip indices, Fields_Hip fields, Fields_Hip B, Fields_Hip G, CeedScalar* W) {\n"; + << "(CeedInt num_elem, void* ctx, FieldsInt_Hip indices, Fields_Hip fields, Fields_Hip B, Fields_Hip G, CeedScalar* W, Points_Hip points) {\n"; // Scratch buffers for (CeedInt i = 0; i < num_input_fields; i++) { @@ -785,6 +1027,10 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << " const CeedInt dim = " << dim << ";\n"; code << " const CeedInt Q_1d = " << Q_1d << ";\n"; + if (is_at_points) { + code << " const CeedInt max_num_points = " << max_num_points << ";\n"; + code << " const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; + } // Shared data code << " extern __shared__ CeedScalar slice[];\n"; @@ -798,11 +1044,13 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // Initialize constants, and matrices B and G code << "\n // Input field constants and basis data\n"; for (CeedInt i = 0; i < num_input_fields; i++) { - CeedCallBackend(CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, use_3d_slices)); + CeedCallBackend( + CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_at_points, use_3d_slices)); } code << "\n // Output field constants and basis data\n"; for (CeedInt i = 0; i < num_output_fields; i++) { - CeedCallBackend(CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); + CeedCallBackend( + CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); } // Loop over all elements @@ -811,7 +1059,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << " for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n"; // -- Compute minimum buffer space needed - CeedInt max_rstr_buffer_size = 0; + CeedInt max_rstr_buffer_size = 1; for (CeedInt i = 0; i < num_input_fields; i++) { CeedInt num_comp, elem_size; @@ -872,12 +1120,16 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { input_field_order[curr_index] = j; curr_index++; } + CeedCallBackend(CeedVectorDestroy(&vec_j)); + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j)); } + CeedCallBackend(CeedVectorDestroy(&vec_i)); + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i)); } } // -- Input restriction and basis - code << " // -- Input field restrictions and basis actions\n"; + code << "\n // -- Input field restrictions and basis actions\n"; for (CeedInt i = 0; i < num_input_fields; i++) { CeedInt f = input_field_order[i]; @@ -885,15 +1137,17 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // ---- Restriction CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, f, dim, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], Q_1d, - true, use_3d_slices)); + true, is_at_points, use_3d_slices)); // ---- Basis action - CeedCallBackend(CeedOperatorBuildKernelBasis_Hip_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, use_3d_slices)); + CeedCallBackend( + CeedOperatorBuildKernelBasis_Hip_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_at_points, use_3d_slices)); } // -- Q function - CeedCallBackend(CeedOperatorBuildKernelQFunction_Hip_gen(code, data, dim, num_input_fields, op_input_fields, qf_input_fields, num_output_fields, - op_output_fields, qf_output_fields, qfunction_name, Q_1d, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelQFunction_Hip_gen(code, data, dim, max_num_points, num_input_fields, op_input_fields, qf_input_fields, + num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_at_points, + use_3d_slices)); // -- Output basis and restriction code << "\n // -- Output field basis action and restrictions\n"; @@ -901,11 +1155,12 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << " // ---- Output field " << i << "\n"; // ---- Basis action - CeedCallBackend(CeedOperatorBuildKernelBasis_Hip_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); + CeedCallBackend( + CeedOperatorBuildKernelBasis_Hip_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); // ---- Restriction - CeedCallBackend( - CeedOperatorBuildKernelRestriction_Hip_gen(code, data, i, dim, NULL, op_output_fields[i], qf_output_fields[i], Q_1d, false, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, i, dim, NULL, op_output_fields[i], qf_output_fields[i], Q_1d, false, + is_at_points, use_3d_slices)); } // Close loop and function @@ -916,6 +1171,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { CeedInt block_sizes[3] = {0, 0, 0}; CeedInt num_elem; + // Compile CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); CeedCallBackend(BlockGridCalculate_Hip_gen(dim, num_elem, data->max_P_1d, Q_1d, block_sizes)); CeedCallBackend(CeedCompile_Hip(ceed, code.str().c_str(), &data->module, 2, "T_1D", block_sizes[0], "BLOCK_SIZE", diff --git a/backends/hip-gen/ceed-hip-gen-operator.c b/backends/hip-gen/ceed-hip-gen-operator.c index 11ee4943b8..59da8d5fbd 100644 --- a/backends/hip-gen/ceed-hip-gen-operator.c +++ b/backends/hip-gen/ceed-hip-gen-operator.c @@ -9,6 +9,7 @@ #include #include #include +#include #include "../hip/ceed-hip-common.h" #include "../hip/ceed-hip-compile.h" @@ -19,10 +20,14 @@ // Destroy operator //------------------------------------------------------------------------------ static int CeedOperatorDestroy_Hip_gen(CeedOperator op) { + Ceed ceed; CeedOperator_Hip_gen *impl; + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); CeedCallBackend(CeedOperatorGetData(op, &impl)); + if (impl->points.num_per_elem) CeedCallHip(ceed, hipFree((void **)impl->points.num_per_elem)); CeedCallBackend(CeedFree(&impl)); + CeedCallBackend(CeedDestroy(&ceed)); return CEED_ERROR_SUCCESS; } @@ -30,6 +35,7 @@ static int CeedOperatorDestroy_Hip_gen(CeedOperator op) { // Apply and add to output //------------------------------------------------------------------------------ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { + bool is_at_points; Ceed ceed; CeedInt num_elem, num_input_fields, num_output_fields; CeedEvalMode eval_mode; @@ -110,6 +116,39 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C } } + // Point coordinates, if needed + CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); + if (is_at_points) { + // Coords + CeedVector vec; + + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec)); + CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->points.coords)); + CeedCallBackend(CeedVectorDestroy(&vec)); + + // Points per elem + if (num_elem != data->points.num_elem) { + CeedInt *points_per_elem; + const CeedInt num_bytes = num_elem * sizeof(CeedInt); + CeedElemRestriction rstr_points = NULL; + + data->points.num_elem = num_elem; + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); + CeedCallBackend(CeedCalloc(num_elem, &points_per_elem)); + for (CeedInt e = 0; e < num_elem; e++) { + CeedInt num_points_elem; + + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points_elem)); + points_per_elem[e] = num_points_elem; + } + if (data->points.num_per_elem) CeedCallHip(ceed, hipFree((void **)data->points.num_per_elem)); + CeedCallHip(ceed, hipMalloc((void **)&data->points.num_per_elem, num_bytes)); + CeedCallHip(ceed, hipMemcpy((void *)data->points.num_per_elem, points_per_elem, num_bytes, hipMemcpyHostToDevice)); + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); + CeedCallBackend(CeedFree(&points_per_elem)); + } + } + // Get context data CeedCallBackend(CeedQFunctionGetInnerContextData(qf, CEED_MEM_DEVICE, &qf_data->d_c)); @@ -163,6 +202,7 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C if (vec == CEED_VECTOR_ACTIVE) vec = output_vec; // Check for multiple output modes CeedInt index = -1; + for (CeedInt j = 0; j < i; j++) { if (vec == output_vecs[j]) { index = j; @@ -175,6 +215,15 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C } } + // Restore point coordinates, if needed + if (is_at_points) { + CeedVector vec; + + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec)); + CeedCallBackend(CeedVectorRestoreArrayRead(vec, &data->points.coords)); + CeedCallBackend(CeedVectorDestroy(&vec)); + } + // Restore context data CeedCallBackend(CeedQFunctionRestoreInnerContextData(qf, &qf_data->d_c)); CeedCallBackend(CeedDestroy(&ceed)); diff --git a/backends/hip-gen/ceed-hip-gen.c b/backends/hip-gen/ceed-hip-gen.c index e867b4eb9e..b7a8c76d0a 100644 --- a/backends/hip-gen/ceed-hip-gen.c +++ b/backends/hip-gen/ceed-hip-gen.c @@ -39,6 +39,7 @@ static int CeedInit_Hip_gen(const char *resource, Ceed ceed) { CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "QFunctionCreate", CeedQFunctionCreate_Hip_gen)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "OperatorCreate", CeedOperatorCreate_Hip_gen)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "OperatorCreateAtPoints", CeedOperatorCreate_Hip_gen)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "Destroy", CeedDestroy_Hip)); return CEED_ERROR_SUCCESS; } diff --git a/backends/hip-gen/ceed-hip-gen.h b/backends/hip-gen/ceed-hip-gen.h index 139bab43bb..eb5dd0c893 100644 --- a/backends/hip-gen/ceed-hip-gen.h +++ b/backends/hip-gen/ceed-hip-gen.h @@ -22,6 +22,7 @@ typedef struct { Fields_Hip B; Fields_Hip G; CeedScalar *W; + Points_Hip points; } CeedOperator_Hip_gen; typedef struct { diff --git a/include/ceed/jit-source/hip/hip-gen-templates.h b/include/ceed/jit-source/hip/hip-gen-templates.h index 56521a3685..358583a79e 100644 --- a/include/ceed/jit-source/hip/hip-gen-templates.h +++ b/include/ceed/jit-source/hip/hip-gen-templates.h @@ -17,6 +17,38 @@ inline __device__ void LoadMatrix(SharedData_Hip &data, const CeedScalar *__rest for (CeedInt i = data.t_id; i < P * Q; i += blockDim.x * blockDim.y * blockDim.z) B[i] = d_B[i]; } +//------------------------------------------------------------------------------ +// AtPoints +//------------------------------------------------------------------------------ + +//------------------------------------------------------------------------------ +// L-vector -> single point +//------------------------------------------------------------------------------ +template +inline __device__ void ReadPoint(SharedData_Hip &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem, + const CeedInt *__restrict__ indices, const CeedScalar *__restrict__ d_u, CeedScalar *r_u) { + const CeedInt ind = indices[p + elem * NUM_PTS]; + + for (CeedInt comp = 0; comp < NUM_COMP; comp++) { + r_u[comp] = d_u[ind + comp * COMP_STRIDE]; + } +} + +//------------------------------------------------------------------------------ +// Single point -> L-vector +//------------------------------------------------------------------------------ +template +inline __device__ void WritePoint(SharedData_Hip &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem, + const CeedInt *__restrict__ indices, const CeedScalar *__restrict__ r_u, CeedScalar *d_u) { + if (p < points_in_elem) { + const CeedInt ind = indices[p + elem * NUM_PTS]; + + for (CeedInt comp = 0; comp < NUM_COMP; comp++) { + d_u[ind + comp * COMP_STRIDE] += r_u[comp]; + } + } +} + //------------------------------------------------------------------------------ // 1D //------------------------------------------------------------------------------ diff --git a/include/ceed/jit-source/hip/hip-types.h b/include/ceed/jit-source/hip/hip-types.h index 418e6fb02c..01fe82e08d 100644 --- a/include/ceed/jit-source/hip/hip-types.h +++ b/include/ceed/jit-source/hip/hip-types.h @@ -23,6 +23,13 @@ typedef struct { CeedInt *outputs[CEED_HIP_NUMBER_FIELDS]; } FieldsInt_Hip; +typedef struct { + CeedInt num_elem; + const CeedInt *num_per_elem; + const CeedInt *indices; + const CeedScalar *coords; +} Points_Hip; + typedef struct { CeedInt t_id_x; CeedInt t_id_y;