1 /// @file 2 /// Test creation and use of FDM element inverse 3 /// \test Test creation and use of FDM element inverse 4 #include <ceed.h> 5 #include <math.h> 6 #include <stdio.h> 7 #include <stdlib.h> 8 9 #include "t540-operator.h" 10 11 int main(int argc, char **argv) { 12 Ceed ceed; 13 CeedElemRestriction elem_restriction_x, elem_restriction_u, elem_restriction_q_data; 14 CeedBasis basis_x, basis_u; 15 CeedQFunction qf_setup_mass, qf_apply; 16 CeedOperator op_setup_mass, op_apply, op_inverse; 17 CeedVector q_data_mass, x, u, v; 18 CeedInt num_elem = 1, p = 4, q = 5, dim = 2; 19 CeedInt num_dofs = p * p, num_qpts = num_elem * q * q; 20 21 CeedInit(argv[1], &ceed); 22 23 // Vectors 24 CeedVectorCreate(ceed, dim * num_elem * (2 * 2), &x); 25 { 26 CeedScalar x_array[dim * num_elem * (2 * 2)]; 27 28 for (CeedInt i = 0; i < 2; i++) { 29 for (CeedInt j = 0; j < 2; j++) { 30 x_array[i + j * 2 + 0 * 4] = i; 31 x_array[i + j * 2 + 1 * 4] = j; 32 } 33 } 34 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 35 } 36 CeedVectorCreate(ceed, num_dofs, &u); 37 CeedVectorCreate(ceed, num_dofs, &v); 38 CeedVectorCreate(ceed, num_qpts, &q_data_mass); 39 40 // Restrictions 41 CeedInt strides_x[3] = {1, 2 * 2, 2 * 2 * dim}; 42 CeedElemRestrictionCreateStrided(ceed, num_elem, 2 * 2, dim, dim * num_elem * 2 * 2, strides_x, &elem_restriction_x); 43 44 CeedInt strides_u[3] = {1, p * p, p * p}; 45 CeedElemRestrictionCreateStrided(ceed, num_elem, p * p, 1, num_dofs, strides_u, &elem_restriction_u); 46 47 CeedInt strides_q_data[3] = {1, q * q, q * q}; 48 CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts, strides_q_data, &elem_restriction_q_data); 49 50 // Bases 51 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, q, CEED_GAUSS, &basis_x); 52 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); 53 54 // QFunction - setup mass 55 CeedQFunctionCreateInterior(ceed, 1, setup_mass, setup_mass_loc, &qf_setup_mass); 56 CeedQFunctionAddInput(qf_setup_mass, "dx", dim * dim, CEED_EVAL_GRAD); 57 CeedQFunctionAddInput(qf_setup_mass, "weight", 1, CEED_EVAL_WEIGHT); 58 CeedQFunctionAddOutput(qf_setup_mass, "q data", 1, CEED_EVAL_NONE); 59 60 // Operator - setup mass 61 CeedOperatorCreate(ceed, qf_setup_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_mass); 62 CeedOperatorSetField(op_setup_mass, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 63 CeedOperatorSetField(op_setup_mass, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 64 CeedOperatorSetField(op_setup_mass, "q data", elem_restriction_q_data, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 65 66 // Apply Setup Operator 67 CeedOperatorApply(op_setup_mass, x, q_data_mass, CEED_REQUEST_IMMEDIATE); 68 69 // QFunction - apply 70 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 71 CeedQFunctionAddInput(qf_apply, "u", 1, CEED_EVAL_INTERP); 72 CeedQFunctionAddInput(qf_apply, "mass q data", 1, CEED_EVAL_NONE); 73 CeedQFunctionAddOutput(qf_apply, "v", 1, CEED_EVAL_INTERP); 74 75 // Operator - apply 76 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply); 77 CeedOperatorSetField(op_apply, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 78 CeedOperatorSetField(op_apply, "mass q data", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data_mass); 79 CeedOperatorSetField(op_apply, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 80 81 // Apply original operator 82 CeedVectorSetValue(u, 1.0); 83 CeedOperatorApply(op_apply, u, v, CEED_REQUEST_IMMEDIATE); 84 85 // Create FDM element inverse 86 CeedOperatorCreateFDMElementInverse(op_apply, &op_inverse, CEED_REQUEST_IMMEDIATE); 87 88 // Apply FDM element inverse 89 CeedOperatorApply(op_inverse, v, u, CEED_REQUEST_IMMEDIATE); 90 91 // Check output 92 { 93 const CeedScalar *u_array; 94 95 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 96 for (int i = 0; i < num_dofs; i++) { 97 if (fabs(u_array[i] - 1.0) > 500. * CEED_EPSILON) { 98 // LCOV_EXCL_START 99 printf("[%" CeedInt_FMT "] Error in inverse: %e - 1.0 = %e\n", i, u_array[i], u_array[i] - 1.); 100 // LCOV_EXCL_STOP 101 } 102 } 103 CeedVectorRestoreArrayRead(u, &u_array); 104 } 105 106 // Cleanup 107 CeedVectorDestroy(&x); 108 CeedVectorDestroy(&q_data_mass); 109 CeedVectorDestroy(&u); 110 CeedVectorDestroy(&v); 111 CeedElemRestrictionDestroy(&elem_restriction_u); 112 CeedElemRestrictionDestroy(&elem_restriction_x); 113 CeedElemRestrictionDestroy(&elem_restriction_q_data); 114 CeedBasisDestroy(&basis_u); 115 CeedBasisDestroy(&basis_x); 116 CeedQFunctionDestroy(&qf_setup_mass); 117 CeedQFunctionDestroy(&qf_apply); 118 CeedOperatorDestroy(&op_setup_mass); 119 CeedOperatorDestroy(&op_apply); 120 CeedOperatorDestroy(&op_inverse); 121 CeedDestroy(&ceed); 122 return 0; 123 } 124