diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index f525460ea6..5347580d2c 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -1109,8 +1109,212 @@ static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op } //------------------------------------------------------------------------------ -// Assemble Operator +// Assemble Operator Diagonal AtPoints //------------------------------------------------------------------------------ +static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) { + bool is_active_at_points = true; + CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, elem_size_active = 1, num_comp_active; + CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; + Ceed ceed; + CeedVector point_coords = NULL, in_vec, out_vec; + CeedElemRestriction rstr_points = NULL; + CeedQFunctionField *qf_input_fields, *qf_output_fields; + CeedQFunction qf; + CeedOperatorField *op_input_fields, *op_output_fields; + CeedOperator_Ref *impl; + + CeedCallBackend(CeedOperatorGetData(op, &impl)); + CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); + CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); + CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); + + // Setup + CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); + + // Ceed + { + Ceed ceed_parent; + + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); + CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); + if (ceed_parent) ceed = ceed_parent; + } + + // Point coordinates + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); + + // Input and output vectors + { + CeedSize input_size, output_size; + + CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); + CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec)); + CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec)); + CeedCallBackend(CeedVectorSetValue(out_vec, 0.0)); + } + + // Input Evecs and Restriction + CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); + + // Check if active field is at points + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedRestrictionType rstr_type; + CeedVector vec; + CeedElemRestriction elem_rstr; + + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + // Skip non-active input + if (vec != CEED_VECTOR_ACTIVE) continue; + + // Get active restriction type + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); + CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); + is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; + if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); + } + + // Loop through elements + for (CeedInt e = 0; e < num_elem; e++) { + CeedInt num_points, e_vec_size = 0; + + // Setup points for element + CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); + + // Input basis apply for non-active bases + CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, + impl->point_coords_elem, true, e_data, impl, request)); + + // Loop over points on element + e_vec_size = (is_active_at_points ? num_points : elem_size_active) * num_comp_active; + for (CeedInt s = 0; s < e_vec_size; s++) { + for (CeedInt i = 0; i < num_input_fields; i++) { + bool is_active_input = false; + CeedInt size; + CeedRestrictionType rstr_type; + CeedEvalMode eval_mode; + CeedVector vec; + CeedElemRestriction elem_rstr; + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + // Skip non-active input + is_active_input = vec == CEED_VECTOR_ACTIVE; + if (!is_active_input) continue; + + // Get elem_size, eval_mode, size + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); + // Update unit vector + { + CeedScalar *array; + + if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); + CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); + array[s] = 1.0; + if (s > 0) array[s - 1] = 0.0; + CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); + } + // Basis action + switch (eval_mode) { + case CEED_EVAL_NONE: + break; + // Note - these basis eval modes require FEM fields + case CEED_EVAL_INTERP: + case CEED_EVAL_GRAD: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], + impl->q_vecs_in[i])); + break; + case CEED_EVAL_WEIGHT: + break; // No action + } + } + + // -- Q function + if (!impl->is_identity_qf) { + CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); + } + + // -- Output basis apply and restriction + CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, + num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); + + // -- Grab diagonal value + for (CeedInt i = 0; i < num_output_fields; i++) { + bool is_active_input = false; + CeedRestrictionType rstr_type; + CeedEvalMode eval_mode; + CeedVector vec; + CeedElemRestriction elem_rstr; + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + // ---- Skip non-active input + is_active_input = vec == CEED_VECTOR_ACTIVE; + if (!is_active_input) continue; + + // ---- Get elem_size, eval_mode, size + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); + // ---- Basis action + switch (eval_mode) { + case CEED_EVAL_NONE: + break; // No action + case CEED_EVAL_INTERP: + case CEED_EVAL_GRAD: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); + CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], + impl->e_vecs_out[i])); + break; + // LCOV_EXCL_START + case CEED_EVAL_WEIGHT: { + return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); + // LCOV_EXCL_STOP + } + } + // ---- Update output vector + { + CeedScalar *array, current_value = 0.0; + + CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); + current_value = array[s]; + CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); + CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[i], 0.0)); + CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); + array[s] = current_value; + CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); + } + // ---- Restrict output block + CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); + if (rstr_type == CEED_RESTRICTION_POINTS) { + CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); + } else { + CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); + } + } + } + num_points_offset += num_points; + } + + // Restore input arrays + CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); + + // Cleanup + CeedCallBackend(CeedVectorDestroy(&in_vec)); + CeedCallBackend(CeedVectorDestroy(&out_vec)); + CeedCallBackend(CeedVectorDestroy(&point_coords)); + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); + return CEED_ERROR_SUCCESS; +} //------------------------------------------------------------------------------ // Operator Destroy @@ -1180,6 +1384,7 @@ int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref)); CeedCallBackend( CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); return CEED_ERROR_SUCCESS; diff --git a/tests/t533-operator.c b/tests/t533-operator.c index a01dabda4a..2a19143ffa 100644 --- a/tests/t533-operator.c +++ b/tests/t533-operator.c @@ -89,7 +89,7 @@ int main(int argc, char **argv) { // Manually assemble diagonal CeedVectorSetValue(u, 0.0); - for (int i = 0; i < num_dofs; i++) { + for (CeedInt i = 0; i < num_dofs; i++) { CeedScalar *u_array; const CeedScalar *v_array; @@ -113,7 +113,7 @@ int main(int argc, char **argv) { const CeedScalar *assembled_array; CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); - for (int i = 0; i < num_dofs; i++) { + for (CeedInt i = 0; i < num_dofs; i++) { if (fabs(assembled_array[i] - assembled_true[i]) > 100. * CEED_EPSILON) { // LCOV_EXCL_START printf("[%" CeedInt_FMT "] Error in assembly: %f != %f\n", i, assembled_array[i], assembled_true[i]); diff --git a/tests/t592-operator.c b/tests/t592-operator.c index e0ccdace2e..91e519f3bb 100644 --- a/tests/t592-operator.c +++ b/tests/t592-operator.c @@ -1,6 +1,6 @@ /// @file /// Test assembly of mass matrix operator QFunction at points -/// \test Test assembly of mass matrix operator QFunction +/// \test Test assembly of mass matrix operator QFunction at points #include #include #include diff --git a/tests/t593-operator.c b/tests/t593-operator.c index 5b145d0884..2a2daceb88 100644 --- a/tests/t593-operator.c +++ b/tests/t593-operator.c @@ -1,5 +1,5 @@ /// @file -/// Bug reproducer for memcheck backends at points +/// Test 1D mass matrix operator at points with heterogeneous points per element /// \test Test 1D mass matrix operator at points with heterogeneous points per element #include #include diff --git a/tests/t594-operator.c b/tests/t594-operator.c new file mode 100644 index 0000000000..2d4e6d876c --- /dev/null +++ b/tests/t594-operator.c @@ -0,0 +1,179 @@ +/// @file +/// Test diagonal assembly of mass matrix operator at points +/// \test Test diagonal assembly of mass matrix operator at points +#include +#include +#include +#include + +#include "t500-operator.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedInt num_elem = 3, dim = 1, p = 3, q = 5; + CeedInt num_nodes_x = num_elem + 1, num_nodes_u = num_elem * (p - 1) + 1, num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; + CeedInt ind_x[num_elem * 2], ind_u[num_elem * p], ind_x_points[num_elem + 1 + num_points]; + CeedScalar x_array_mesh[num_nodes_x], x_array_points[num_points], assembled_true[num_nodes_u]; + CeedVector x_points = NULL, x_elem = NULL, q_data = NULL, u = NULL, v = NULL, assembled = NULL; + CeedElemRestriction elem_restriction_x_points, elem_restriction_q_data, elem_restriction_x, elem_restriction_u; + CeedBasis basis_x, basis_u; + CeedQFunction qf_setup, qf_mass; + CeedOperator op_setup, op_mass; + bool is_at_points; + + CeedInit(argv[1], &ceed); + + // Mesh coordinates + for (CeedInt i = 0; i < num_nodes_x; i++) x_array_mesh[i] = (CeedScalar)i / (num_nodes_x - 1); + for (CeedInt i = 0; i < num_elem; i++) { + ind_x[2 * i + 0] = i; + ind_x[2 * i + 1] = i + 1; + } + CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_nodes_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); + CeedVectorCreate(ceed, num_nodes_x, &x_elem); + CeedVectorSetArray(x_elem, CEED_MEM_HOST, CEED_USE_POINTER, x_array_mesh); + + // U mesh + for (CeedInt i = 0; i < num_elem; i++) { + for (CeedInt j = 0; j < p; j++) { + ind_u[p * i + j] = i * (p - 1) + j; + } + } + CeedElemRestrictionCreate(ceed, num_elem, p, 1, 1, num_nodes_u, CEED_MEM_HOST, CEED_USE_POINTER, ind_u, &elem_restriction_u); + + // Point reference coordinates + { + CeedScalar weight_tmp[num_points_per_elem + 1]; + CeedInt current_index = 0; + + // Use num_points_per_elem + 1 to test non-uniform quadrature + CeedGaussQuadrature(num_points_per_elem + 1, x_array_points, weight_tmp); + ind_x_points[0] = num_elem + 1; + for (CeedInt p = 0; p < num_points_per_elem + 1; p++, current_index++) { + ind_x_points[num_elem + 1 + current_index] = current_index; + } + // Use num_points_per_elem for middle elements + for (CeedInt e = 1; e < num_elem - 1; e++) { + CeedGaussQuadrature(num_points_per_elem, &x_array_points[current_index], weight_tmp); + ind_x_points[e] = num_elem + 1 + current_index; + for (CeedInt p = 0; p < num_points_per_elem; p++, current_index++) { + ind_x_points[num_elem + 1 + current_index] = current_index; + } + } + // Use num_points_per_elem - 1 to test non-uniform quadrature + CeedGaussQuadrature(num_points_per_elem - 1, &x_array_points[current_index], weight_tmp); + ind_x_points[num_elem - 1] = num_elem + 1 + current_index; + for (CeedInt p = 0; p < num_points_per_elem - 1; p++, current_index++) { + ind_x_points[num_elem + 1 + current_index] = current_index; + } + ind_x_points[num_elem] = num_elem + 1 + current_index; + + CeedVectorCreate(ceed, num_elem * num_points_per_elem, &x_points); + CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_USE_POINTER, x_array_points); + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x_points, + &elem_restriction_x_points); + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x_points, + &elem_restriction_q_data); + + // Q data + CeedVectorCreate(ceed, num_points, &q_data); + } + + // Basis creation + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, q, CEED_GAUSS, &basis_x); + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); + + // Setup geometric scaling + CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); + CeedQFunctionAddInput(qf_setup, "x", dim * dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); + + CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); + CeedOperatorSetField(op_setup, "x", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_setup, elem_restriction_x_points, x_points); + + CeedOperatorApply(op_setup, x_elem, q_data, CEED_REQUEST_IMMEDIATE); + + // Mass operator + CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); + CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); + CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP); + + CeedOperatorCreateAtPoints(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); + CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); + CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_mass, elem_restriction_x_points, x_points); + + CeedOperatorIsAtPoints(op_mass, &is_at_points); + if (!is_at_points) printf("Error: Operator should be at points\n"); + + CeedVectorCreate(ceed, num_nodes_u, &u); + CeedVectorSetValue(u, 0.0); + CeedVectorCreate(ceed, num_nodes_u, &v); + + // Assemble diagonal + CeedVectorCreate(ceed, num_nodes_u, &assembled); + CeedOperatorLinearAssembleDiagonal(op_mass, assembled, CEED_REQUEST_IMMEDIATE); + + // Manually assemble diagonal + CeedVectorSetValue(u, 0.0); + for (CeedInt i = 0; i < num_nodes_u; i++) { + CeedScalar *u_array; + const CeedScalar *v_array; + + // Set input + CeedVectorGetArray(u, CEED_MEM_HOST, &u_array); + u_array[i] = 1.0; + if (i) u_array[i - 1] = 0.0; + CeedVectorRestoreArray(u, &u_array); + + // Compute diag entry for DoF i + CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE); + + // Retrieve entry + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + assembled_true[i] = v_array[i]; + CeedVectorRestoreArrayRead(v, &v_array); + } + + // Check output + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); + for (CeedInt i = 0; i < num_nodes_u; i++) { + if (fabs(assembled_array[i] - assembled_true[i]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT "] Error in assembly: %f != %f\n", i, assembled_array[i], assembled_true[i]); + // LCOV_EXCL_STOP + } + } + CeedVectorRestoreArrayRead(assembled, &assembled_array); + } + + // Cleanup + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&x_elem); + CeedVectorDestroy(&assembled); + CeedElemRestrictionDestroy(&elem_restriction_q_data); + CeedElemRestrictionDestroy(&elem_restriction_x); + CeedElemRestrictionDestroy(&elem_restriction_x_points); + CeedElemRestrictionDestroy(&elem_restriction_u); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + CeedQFunctionDestroy(&qf_setup); + CeedQFunctionDestroy(&qf_mass); + CeedOperatorDestroy(&op_setup); + CeedOperatorDestroy(&op_mass); + CeedDestroy(&ceed); + return 0; +}