diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 9adc6e4c3c..89fd47df76 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -1524,6 +1524,212 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, Ce return CEED_ERROR_SUCCESS; } +//------------------------------------------------------------------------------ +// Assemble Operator AtPoints +//------------------------------------------------------------------------------ +static int CeedSingleOperatorAssembleAtPoints_Ref(CeedOperator op, CeedInt offset, CeedVector values) { + CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1; + CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}, *assembled; + 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)); + CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed)); + CeedCallBackend(CeedDestroy(&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)); + } + + // Assembled array + CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_HOST, &assembled)); + + // Clear input Evecs + for (CeedInt i = 0; i < num_input_fields; i++) { + bool is_active; + CeedVector vec; + + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + is_active = vec == CEED_VECTOR_ACTIVE; + CeedCallBackend(CeedVectorDestroy(&vec)); + if (!is_active || impl->skip_rstr_in[i]) continue; + CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); + } + + // Input Evecs and Restriction + CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, CEED_REQUEST_IMMEDIATE)); + + // 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, CEED_REQUEST_IMMEDIATE)); + 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, false, e_data, impl, CEED_REQUEST_IMMEDIATE)); + + // Loop over points on element + for (CeedInt i = 0; i < num_input_fields; i++) { + bool is_active_at_points = true, is_active; + CeedInt elem_size_active = 1; + CeedRestrictionType rstr_type; + CeedVector vec; + CeedElemRestriction elem_rstr; + + // -- Skip non-active input + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + is_active = vec == CEED_VECTOR_ACTIVE; + CeedCallBackend(CeedVectorDestroy(&vec)); + if (!is_active || impl->skip_rstr_in[i]) continue; + + // -- Get active restriction type + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); + is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; + if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); + else elem_size_active = num_points; + CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); + + e_vec_size = elem_size_active * num_comp_active; + for (CeedInt s = 0; s < e_vec_size; s++) { + const CeedInt comp_in = s / elem_size_active; + const CeedInt node_in = s % elem_size_active; + + // -- 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)); + } + // Input basis apply for 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, false, true, e_data, impl, CEED_REQUEST_IMMEDIATE)); + + // -- 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, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec, + impl->point_coords_elem, true, impl, CEED_REQUEST_IMMEDIATE)); + + // -- Build element matrix + for (CeedInt j = 0; j < num_output_fields; j++) { + bool is_active; + CeedInt elem_size = 0; + CeedRestrictionType rstr_type; + CeedVector vec; + CeedElemRestriction elem_rstr; + + // ---- Skip non-active output + CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec)); + is_active = vec == CEED_VECTOR_ACTIVE; + CeedCallBackend(CeedVectorDestroy(&vec)); + if (!is_active || impl->skip_rstr_out[j]) continue; + + // ---- Check if elem size matches + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); + if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) { + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); + continue; + } + if (rstr_type == CEED_RESTRICTION_POINTS) { + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size)); + } else { + CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); + } + { + CeedInt num_comp = 0; + + CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); + if (e_vec_size != num_comp * elem_size) continue; + } + // ---- Copy output + { + const CeedScalar *output; + + CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_out[j], CEED_MEM_HOST, &output)); + for (CeedInt k = 0; k < e_vec_size; k++) { + const CeedInt comp_out = k / elem_size_active; + const CeedInt node_out = k % elem_size_active; + + assembled[offset + e * e_vec_size * e_vec_size + (comp_in * num_comp_active + comp_out) * elem_size_active * elem_size_active + + node_out * elem_size_active + node_in] = output[k]; + } + CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_out[j], &output)); + } + } + // -- Reset unit vector + if (s == e_vec_size - 1) { + CeedScalar *array; + + CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); + array[s] = 0.0; + CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); + } + } + } + num_points_offset += num_points; + } + + // Restore input arrays + CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); + + // Restore assembled values + CeedCallBackend(CeedVectorRestoreArray(values, &assembled)); + + // Cleanup + CeedCallBackend(CeedDestroy(&ceed)); + CeedCallBackend(CeedVectorDestroy(&in_vec)); + CeedCallBackend(CeedVectorDestroy(&out_vec)); + CeedCallBackend(CeedVectorDestroy(&point_coords)); + CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); + CeedCallBackend(CeedQFunctionDestroy(&qf)); + return CEED_ERROR_SUCCESS; +} + //------------------------------------------------------------------------------ // Operator Destroy //------------------------------------------------------------------------------ @@ -1592,6 +1798,7 @@ int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { CeedCallBackend( CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssembleAtPoints_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); CeedCallBackend(CeedDestroy(&ceed)); diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index 970f3549b5..5de71fed38 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -567,7 +567,7 @@ static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, C @ref Developer **/ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVector values) { - bool is_composite; + bool is_composite, is_at_points; CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); @@ -595,6 +595,10 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVecto } } + CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); + CeedCheck(!is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, + "Backend does not implement CeedOperatorLinearAssemble for AtPoints operator"); + // Assemble QFunction CeedInt layout_qf[3]; const CeedScalar *assembled_qf_array; diff --git a/tests/t596-operator.c b/tests/t596-operator.c new file mode 100644 index 0000000000..81ca865ebd --- /dev/null +++ b/tests/t596-operator.c @@ -0,0 +1,202 @@ +/// @file +/// Test full assembly of mass matrix operator +/// \test Test full assembly of mass matrix operator AtPoints +#include +#include +#include +#include + +#include "t596-operator.h" + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt num_comp = 1; num_comp <= 3; num_comp++) { + CeedElemRestriction elem_restriction_x, elem_restriction_x_points, elem_restriction_u, elem_restriction_q_data; + CeedBasis basis_x, basis_u; + CeedQFunction qf_setup, qf_mass; + CeedOperator op_setup, op_mass; + CeedVector q_data, x, x_points, u, v; + CeedInt p = 3, q = 4, dim = 2; + CeedInt n_x = 3, n_y = 2; + CeedInt num_elem = n_x * n_y; + CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; + CeedInt ind_x[num_elem * p * p]; + CeedScalar assembled_values[num_comp * num_comp * num_dofs * num_dofs]; + CeedScalar assembled_true[num_comp * num_comp * num_dofs * num_dofs]; + + // Points + CeedVectorCreate(ceed, dim * num_points, &x_points); + { + CeedScalar x_array[dim * num_points]; + + for (CeedInt e = 0; e < num_elem; e++) { + for (CeedInt d = 0; d < dim; d++) { + x_array[num_points_per_elem * (e * dim + d) + 0] = 0.25; + x_array[num_points_per_elem * (e * dim + d) + 1] = d == 0 ? -0.25 : 0.25; + x_array[num_points_per_elem * (e * dim + d) + 2] = d == 0 ? 0.25 : -0.25; + x_array[num_points_per_elem * (e * dim + d) + 3] = 0.25; + } + } + CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + { + CeedInt ind_x[num_elem + 1 + num_points]; + + for (CeedInt i = 0; i <= num_elem; i++) ind_x[i] = num_elem + 1 + i * num_points_per_elem; + for (CeedInt i = 0; i < num_points; i++) ind_x[num_elem + 1 + i] = i; + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, + &elem_restriction_x_points); + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_q_data); + } + + // Vectors + CeedVectorCreate(ceed, dim * num_dofs, &x); + { + CeedScalar x_array[dim * num_dofs]; + + for (CeedInt i = 0; i < n_x * 2 + 1; i++) { + for (CeedInt j = 0; j < n_y * 2 + 1; j++) { + x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x); + x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y); + } + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedVectorCreate(ceed, num_comp * num_dofs, &u); + CeedVectorCreate(ceed, num_comp * num_dofs, &v); + CeedVectorCreate(ceed, num_points, &q_data); + + // Restrictions + for (CeedInt i = 0; i < num_elem; i++) { + CeedInt col, row, offset; + col = i % n_x; + row = i / n_x; + offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1); + for (CeedInt j = 0; j < p; j++) { + for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j; + } + } + CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); + CeedElemRestrictionCreate(ceed, num_elem, p * p, num_comp, num_dofs, num_comp * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, + &elem_restriction_u); + + // Bases + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x); + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u); + + // QFunctions + CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); + CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); + CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); + + CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); + CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); + CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP); + CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP); + { + CeedQFunctionContext qf_context; + + CeedQFunctionContextCreate(ceed, &qf_context); + CeedQFunctionContextSetData(qf_context, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(CeedInt), &num_comp); + CeedQFunctionSetContext(qf_mass, qf_context); + CeedQFunctionContextDestroy(&qf_context); + } + + // Operators + CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); + CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_setup, elem_restriction_x_points, x_points); + + CeedOperatorCreateAtPoints(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); + CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); + CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_mass, elem_restriction_x_points, x_points); + + // Apply Setup Operator + CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); + + // Fully assemble operator + CeedSize num_entries; + CeedInt *rows; + CeedInt *cols; + CeedVector assembled; + + for (CeedInt k = 0; k < num_comp * num_comp * num_dofs * num_dofs; ++k) { + assembled_values[k] = 0.0; + assembled_true[k] = 0.0; + } + CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols); + CeedVectorCreate(ceed, num_entries, &assembled); + CeedOperatorLinearAssemble(op_mass, assembled); + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); + for (CeedInt k = 0; k < num_entries; k++) { + assembled_values[rows[k] * num_comp * num_dofs + cols[k]] += assembled_array[k]; + } + CeedVectorRestoreArrayRead(assembled, &assembled_array); + } + + // Manually assemble operator + CeedVectorSetValue(u, 0.0); + for (CeedInt j = 0; j < num_comp * num_dofs; j++) { + CeedScalar *u_array; + const CeedScalar *v_array; + + // Set input + CeedVectorGetArray(u, CEED_MEM_HOST, &u_array); + u_array[j] = 1.0; + if (j) u_array[j - 1] = 0.0; + CeedVectorRestoreArray(u, &u_array); + + // Compute entries for column j + CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE); + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_comp * num_dofs; i++) assembled_true[i * num_comp * num_dofs + j] = v_array[i]; + CeedVectorRestoreArrayRead(v, &v_array); + } + + // Check output + for (CeedInt i = 0; i < num_comp * num_dofs; i++) { + for (CeedInt j = 0; j < num_comp * num_dofs; j++) { + if (fabs(assembled_values[i * num_dofs * num_comp + j] - assembled_true[i * num_dofs * num_comp + j]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_dofs * num_comp + j], + assembled_true[i * num_dofs * num_comp + j]); + // LCOV_EXCL_STOP + } + } + } + + // Cleanup + free(rows); + free(cols); + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&assembled); + CeedElemRestrictionDestroy(&elem_restriction_u); + CeedElemRestrictionDestroy(&elem_restriction_x); + CeedElemRestrictionDestroy(&elem_restriction_x_points); + CeedElemRestrictionDestroy(&elem_restriction_q_data); + CeedBasisDestroy(&basis_u); + CeedBasisDestroy(&basis_x); + CeedQFunctionDestroy(&qf_setup); + CeedQFunctionDestroy(&qf_mass); + CeedOperatorDestroy(&op_setup); + CeedOperatorDestroy(&op_mass); + } + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t596-operator.h b/tests/t596-operator.h new file mode 100644 index 0000000000..4fe3d700f5 --- /dev/null +++ b/tests/t596-operator.h @@ -0,0 +1,29 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +#include + +CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + const CeedScalar *weight = in[0], *J = in[1]; + CeedScalar *rho = out[0]; + + for (CeedInt i = 0; i < Q; i++) { + rho[i] = weight[i] * (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]); + } + return 0; +} + +CEED_QFUNCTION(mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + CeedInt num_comp = *(CeedInt *)ctx; + const CeedScalar *rho = in[0], *u = in[1]; + CeedScalar *v = out[0]; + + for (CeedInt i = 0; i < Q; i++) { + for (CeedInt c = 0; c < num_comp; c++) v[i + c * Q] = rho[i] * c * u[i + c * Q]; + } + return 0; +} diff --git a/tests/t597-operator.c b/tests/t597-operator.c new file mode 100644 index 0000000000..25d6b3cf3f --- /dev/null +++ b/tests/t597-operator.c @@ -0,0 +1,203 @@ +/// @file +/// Test full assembly of Poisson operator AtPoints +/// \test Test full assembly of Poisson operator AtPoints +#include +#include +#include +#include + +#include "t597-operator.h" + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt num_comp = 1; num_comp <= 3; num_comp++) { + CeedElemRestriction elem_restriction_x, elem_restriction_x_points, elem_restriction_u, elem_restriction_q_data; + CeedBasis basis_x, basis_u; + CeedQFunction qf_setup, qf_diff; + CeedOperator op_setup, op_diff; + CeedVector q_data, x, x_points, u, v; + CeedInt p = 3, q = 4, dim = 2; + CeedInt n_x = 3, n_y = 2; + CeedInt num_elem = n_x * n_y; + CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; + CeedInt ind_x[num_elem * p * p]; + CeedScalar assembled_values[num_comp * num_comp * num_dofs * num_dofs]; + CeedScalar assembled_true[num_comp * num_comp * num_dofs * num_dofs]; + + // Points + CeedVectorCreate(ceed, dim * num_points, &x_points); + { + CeedScalar x_array[dim * num_points]; + + for (CeedInt e = 0; e < num_elem; e++) { + for (CeedInt d = 0; d < dim; d++) { + x_array[num_points_per_elem * (e * dim + d) + 0] = 0.25; + x_array[num_points_per_elem * (e * dim + d) + 1] = d == 0 ? -0.25 : 0.25; + x_array[num_points_per_elem * (e * dim + d) + 2] = d == 0 ? 0.25 : -0.25; + x_array[num_points_per_elem * (e * dim + d) + 3] = 0.25; + } + } + CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + { + CeedInt ind_x[num_elem + 1 + num_points]; + + for (CeedInt i = 0; i <= num_elem; i++) ind_x[i] = num_elem + 1 + i * num_points_per_elem; + for (CeedInt i = 0; i < num_points; i++) ind_x[num_elem + 1 + i] = i; + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, + &elem_restriction_x_points); + CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim * (dim + 1) / 2, num_points * dim * (dim + 1) / 2, CEED_MEM_HOST, + CEED_COPY_VALUES, ind_x, &elem_restriction_q_data); + } + + // Vectors + CeedVectorCreate(ceed, dim * num_dofs, &x); + { + CeedScalar x_array[dim * num_dofs]; + + for (CeedInt i = 0; i < n_x * 2 + 1; i++) { + for (CeedInt j = 0; j < n_y * 2 + 1; j++) { + x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x); + x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y); + } + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedVectorCreate(ceed, num_comp * num_dofs, &u); + CeedVectorCreate(ceed, num_comp * num_dofs, &v); + CeedVectorCreate(ceed, num_points * dim * (dim + 1) / 2, &q_data); + + // Restrictions + for (CeedInt i = 0; i < num_elem; i++) { + CeedInt col, row, offset; + col = i % n_x; + row = i / n_x; + offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1); + for (CeedInt j = 0; j < p; j++) { + for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j; + } + } + CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); + CeedElemRestrictionCreate(ceed, num_elem, p * p, num_comp, num_dofs, num_comp * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, + &elem_restriction_u); + + // Bases + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x); + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u); + + // QFunction - setup + CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); + CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddOutput(qf_setup, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE); + + // Operator - setup + CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); + CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_setup, "q data", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_setup, elem_restriction_x_points, x_points); + + // Apply Setup Operator + CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); + + // QFunction - apply + CeedQFunctionCreateInterior(ceed, 1, diff, diff_loc, &qf_diff); + CeedQFunctionAddInput(qf_diff, "du", num_comp * dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_diff, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE); + CeedQFunctionAddOutput(qf_diff, "dv", num_comp * dim, CEED_EVAL_GRAD); + { + CeedQFunctionContext qf_context; + + CeedQFunctionContextCreate(ceed, &qf_context); + CeedQFunctionContextSetData(qf_context, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(CeedInt), &num_comp); + CeedQFunctionSetContext(qf_diff, qf_context); + CeedQFunctionContextDestroy(&qf_context); + } + + // Operator - apply + CeedOperatorCreateAtPoints(ceed, qf_diff, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff); + CeedOperatorSetField(op_diff, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_diff, "q data", elem_restriction_q_data, CEED_BASIS_NONE, q_data); + CeedOperatorSetField(op_diff, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorAtPointsSetPoints(op_diff, elem_restriction_x_points, x_points); + + // Fully assemble operator + CeedSize num_entries; + CeedInt *rows; + CeedInt *cols; + CeedVector assembled; + + for (CeedInt k = 0; k < num_comp * num_comp * num_dofs * num_dofs; ++k) { + assembled_values[k] = 0.0; + assembled_true[k] = 0.0; + } + CeedOperatorLinearAssembleSymbolic(op_diff, &num_entries, &rows, &cols); + CeedVectorCreate(ceed, num_entries, &assembled); + CeedOperatorLinearAssemble(op_diff, assembled); + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); + for (CeedInt k = 0; k < num_entries; k++) assembled_values[rows[k] * num_comp * num_dofs + cols[k]] += assembled_array[k]; + CeedVectorRestoreArrayRead(assembled, &assembled_array); + } + + // Manually assemble operator + CeedVectorSetValue(u, 0.0); + for (CeedInt j = 0; j < num_comp * num_dofs; j++) { + CeedScalar *u_array; + const CeedScalar *v_array; + + // Set input + CeedVectorGetArray(u, CEED_MEM_HOST, &u_array); + u_array[j] = 1.0; + if (j) u_array[j - 1] = 0.0; + CeedVectorRestoreArray(u, &u_array); + + // Compute entries for column j + CeedOperatorApply(op_diff, u, v, CEED_REQUEST_IMMEDIATE); + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_comp * num_dofs; i++) assembled_true[i * num_comp * num_dofs + j] = v_array[i]; + CeedVectorRestoreArrayRead(v, &v_array); + } + + // Check output + for (CeedInt i = 0; i < num_comp * num_dofs; i++) { + for (CeedInt j = 0; j < num_comp * num_dofs; j++) { + if (fabs(assembled_values[i * num_comp * num_dofs + j] - assembled_true[i * num_comp * num_dofs + j]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_comp * num_dofs + j], + assembled_true[i * num_comp * num_dofs + j]); + // LCOV_EXCL_STOP + } + } + } + + // Cleanup + free(rows); + free(cols); + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&q_data); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&assembled); + CeedElemRestrictionDestroy(&elem_restriction_u); + CeedElemRestrictionDestroy(&elem_restriction_x); + CeedElemRestrictionDestroy(&elem_restriction_x_points); + CeedElemRestrictionDestroy(&elem_restriction_q_data); + CeedBasisDestroy(&basis_u); + CeedBasisDestroy(&basis_x); + CeedQFunctionDestroy(&qf_setup); + CeedQFunctionDestroy(&qf_diff); + CeedOperatorDestroy(&op_setup); + CeedOperatorDestroy(&op_diff); + } + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t597-operator.h b/tests/t597-operator.h new file mode 100644 index 0000000000..b68854b6fe --- /dev/null +++ b/tests/t597-operator.h @@ -0,0 +1,59 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +#include + +CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store + // the symmetric part of the result. + + // in[0] is Jacobians with shape [2, nc=2, Q] + // in[1] is quadrature weights, size (Q) + const CeedScalar *J = in[0], *qw = in[1]; + + // out[0] is qdata, size (Q) + CeedScalar *qd = out[0]; + + // Quadrature point loop + for (CeedInt i = 0; i < Q; i++) { + // J: 0 2 qd: 0 2 adj(J): J22 -J12 + // 1 3 2 1 -J21 J11 + const CeedScalar J11 = J[i + Q * 0]; + const CeedScalar J21 = J[i + Q * 1]; + const CeedScalar J12 = J[i + Q * 2]; + const CeedScalar J22 = J[i + Q * 3]; + const CeedScalar w = qw[i] / (J11 * J22 - J21 * J12); + qd[i + Q * 0] = w * (J12 * J12 + J22 * J22); + qd[i + Q * 2] = w * (J11 * J11 + J21 * J21); + qd[i + Q * 1] = -w * (J11 * J12 + J21 * J22); + } + + return 0; +} + +CEED_QFUNCTION(diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + CeedInt num_comp = *(CeedInt *)ctx; + // in[0] is gradient u, shape [2, nc=1, Q] + // in[1] is quadrature data, size (3*Q) + const CeedScalar *du = in[0], *qd = in[1]; + + // out[0] is output to multiply against gradient v, shape [2, nc=1, Q] + CeedScalar *dv = out[0]; + + // Quadrature point loop + for (CeedInt i = 0; i < Q; i++) { + for (CeedInt c = 0; c < num_comp; c++) { + const CeedScalar du0 = du[i + Q * (2 * c + 0)]; + const CeedScalar du1 = du[i + Q * (2 * c + 1)]; + + dv[i + Q * (2 * c + 0)] = qd[i + Q * 0] * du0 + qd[i + Q * 2] * du1; + dv[i + Q * (2 * c + 1)] = qd[i + Q * 2] * du0 + qd[i + Q * 1] * du1; + } + } + + return 0; +}