From b37f88251c8f416841969dcaf07a52601ea17f62 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Mon, 11 Mar 2024 15:41:30 -0600 Subject: [PATCH 1/4] cpu - initalize scratch vectors in AtPoints Operator impl --- backends/ref/ceed-ref-operator.c | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 666e3c6a16..8bbca38a6f 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -585,7 +585,10 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); CeedCallBackend(CeedOperatorGetData(op, &impl)); - if (is_input) CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); + if (is_input) { + CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); + CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0)); + } } // Loop over fields @@ -614,6 +617,7 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE || !is_input) { CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); + CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); } else { q_size = (CeedSize)max_num_points * size; CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); @@ -632,6 +636,7 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); q_size = (CeedSize)max_num_points * size; CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); + CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); break; case CEED_EVAL_WEIGHT: // Only on input fields CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); From bff3dd9cf7bc164f66dded8d46268babbba91b73 Mon Sep 17 00:00:00 2001 From: Zach Atkins Date: Mon, 11 Mar 2024 15:24:37 -0600 Subject: [PATCH 2/4] Add reproducer for memcheck bug in CeedOperatorInputBasisAtPoints_Ref --- tests/t593-operator.c | 152 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 152 insertions(+) create mode 100644 tests/t593-operator.c diff --git a/tests/t593-operator.c b/tests/t593-operator.c new file mode 100644 index 0000000000..5b145d0884 --- /dev/null +++ b/tests/t593-operator.c @@ -0,0 +1,152 @@ +/// @file +/// Bug reproducer for memcheck backends at points +/// \test Test 1D mass matrix operator at points with heterogeneous points per element +#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]; + CeedVector x_points = NULL, x_elem = NULL, q_data = NULL, u = NULL, v = 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 QFunction + CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE); + + // Check output + { + const CeedScalar *v_array; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_nodes_u; i++) { + if (fabs(v_array[i]) > 1e-14) printf("[%" CeedInt_FMT "] v %g != 0.0\n", i, v_array[i]); + } + CeedVectorRestoreArrayRead(v, &v_array); + } + + // Cleanup + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&x_elem); + 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; +} From 23629c235d530cd8e8aca9514a73f50bdf0ec8af Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Mon, 11 Mar 2024 16:22:42 -0600 Subject: [PATCH 3/4] cpu - padd AtPoints evec for memcheck --- backends/ref/ceed-ref-operator.c | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 8bbca38a6f..78344e5f6c 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -19,7 +19,8 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { Ceed ceed; - CeedSize e_size, q_size; + bool is_at_points; + CeedSize e_size, q_size, at_points_e_size_padding = 0; CeedInt num_comp, size, P; CeedQFunctionField *qf_fields; CeedOperatorField *op_fields; @@ -31,6 +32,7 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); if (ceed_parent) ceed = ceed_parent; } + CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); if (is_input) { CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); @@ -48,7 +50,24 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); if (eval_mode != CEED_EVAL_WEIGHT) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); - CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); + if (is_at_points) { + CeedSize e_size; + + CeedCallBackend(CeedElemRestrictionGetEVectorSize(elem_rstr, &e_size)); + if (at_points_e_size_padding == 0) { + CeedInt num_elem, num_points, max_points, num_comp; + + CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, num_elem - 1, &num_points)); + CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(elem_rstr, &max_points)); + CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); + at_points_e_size_padding = (max_points - num_points) * num_comp; + } + CeedCallBackend(CeedVectorCreate(ceed, e_size + at_points_e_size_padding, &e_vecs_full[i + start_e])); + CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); + } else { + CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); + } } switch (eval_mode) { From 6efa0d72b0004d844b67396ef4fa3300cbbcc0f4 Mon Sep 17 00:00:00 2001 From: Zach Atkins Date: Mon, 11 Mar 2024 16:49:02 -0600 Subject: [PATCH 4/4] Move fix to CeedOperatorSetupFieldsAtPoints_Ref --- backends/ref/ceed-ref-operator.c | 43 +++++++++++++++----------------- 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 78344e5f6c..3a3780bcf3 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -20,7 +20,7 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { Ceed ceed; bool is_at_points; - CeedSize e_size, q_size, at_points_e_size_padding = 0; + CeedSize e_size, q_size; CeedInt num_comp, size, P; CeedQFunctionField *qf_fields; CeedOperatorField *op_fields; @@ -50,24 +50,7 @@ static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool i CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); if (eval_mode != CEED_EVAL_WEIGHT) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); - if (is_at_points) { - CeedSize e_size; - - CeedCallBackend(CeedElemRestrictionGetEVectorSize(elem_rstr, &e_size)); - if (at_points_e_size_padding == 0) { - CeedInt num_elem, num_points, max_points, num_comp; - - CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); - CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, num_elem - 1, &num_points)); - CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(elem_rstr, &max_points)); - CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); - at_points_e_size_padding = (max_points - num_points) * num_comp; - } - CeedCallBackend(CeedVectorCreate(ceed, e_size + at_points_e_size_padding, &e_vecs_full[i + start_e])); - CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); - } else { - CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); - } + CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); } switch (eval_mode) { @@ -574,7 +557,7 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { Ceed ceed; CeedSize e_size, q_size; - CeedInt max_num_points, num_comp, size, P; + CeedInt e_size_padding = 0, max_num_points, num_comp, size, P; CeedQFunctionField *qf_fields; CeedOperatorField *op_fields; @@ -617,13 +600,27 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); if (eval_mode != CEED_EVAL_WEIGHT) { - CeedRestrictionType rstr_type; CeedElemRestriction elem_rstr; + CeedSize e_size; + bool is_at_points; CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); - CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); - CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); + CeedCallBackend(CeedElemRestrictionIsPoints(elem_rstr, &is_at_points)); + if (is_at_points) { + CeedCallBackend(CeedElemRestrictionGetEVectorSize(elem_rstr, &e_size)); + if (e_size_padding == 0) { + CeedInt num_points, num_elem; + + CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem)); + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, num_elem - 1, &num_points)); + e_size_padding = (max_num_points - num_points) * num_comp; + } + CeedCallBackend(CeedVectorCreate(ceed, e_size + e_size_padding, &e_vecs_full[i + start_e])); + CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); + } else { + CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); + } } switch (eval_mode) {