1 /// @file 2 /// Test creation, action, and destruction for mass matrix operator with multigrid level, nontensor basis 3 /// \test Test creation, action, and destruction for mass matrix operator with multigrid level, nontensor basis 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 8 int main(int argc, char **argv) { 9 Ceed ceed; 10 CeedElemRestriction elem_restr_x, elem_restr_qd, 11 elem_restr_u_c, elem_restr_u_f; 12 CeedBasis basis_x, basis_u; 13 CeedQFunction qf_setup, qf_mass; 14 CeedOperator op_setup, op_mass_c, op_mass_f, 15 op_prolong, op_restrict; 16 CeedVector q_data, X, U_c, U_f, 17 V_c, V_f, p_mult_f; 18 const CeedScalar *hv; 19 CeedInt num_elem = 15, P_c = 3, P_f = 5, Q = 8; 20 CeedInt num_dofs_x = num_elem+1, num_dofs_u_c = num_elem*(P_c-1)+1, 21 num_dofs_u_f = num_elem*(P_f-1)+1; 22 CeedInt ind_u_c[num_elem*P_c], ind_u_f[num_elem*P_f], 23 ind_x[num_elem*2]; 24 CeedScalar x[num_dofs_x]; 25 CeedScalar sum; 26 27 CeedInit(argv[1], &ceed); 28 29 for (CeedInt i=0; i<num_dofs_x; i++) 30 x[i] = (CeedScalar) i / (num_dofs_x - 1); 31 for (CeedInt i=0; i<num_elem; i++) { 32 ind_x[2*i+0] = i; 33 ind_x[2*i+1] = i+1; 34 } 35 // Restrictions 36 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_dofs_x, CEED_MEM_HOST, 37 CEED_USE_POINTER, ind_x, &elem_restr_x); 38 39 for (CeedInt i=0; i<num_elem; i++) { 40 for (CeedInt j=0; j<P_c; j++) { 41 ind_u_c[P_c*i+j] = i*(P_c-1) + j; 42 } 43 } 44 CeedElemRestrictionCreate(ceed, num_elem, P_c, 1, 1, num_dofs_u_c, 45 CEED_MEM_HOST, 46 CEED_USE_POINTER, ind_u_c, &elem_restr_u_c); 47 48 for (CeedInt i=0; i<num_elem; i++) { 49 for (CeedInt j=0; j<P_f; j++) { 50 ind_u_f[P_f*i+j] = i*(P_f-1) + j; 51 } 52 } 53 CeedElemRestrictionCreate(ceed, num_elem, P_f, 1, 1, num_dofs_u_f, 54 CEED_MEM_HOST, 55 CEED_USE_POINTER, ind_u_f, &elem_restr_u_f); 56 57 CeedInt strides_qd[3] = {1, Q, Q}; 58 CeedElemRestrictionCreateStrided(ceed, num_elem, Q, 1, Q*num_elem, strides_qd, 59 &elem_restr_qd); 60 61 // Bases 62 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &basis_x); 63 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P_f, Q, CEED_GAUSS, &basis_u); 64 65 // QFunctions 66 CeedQFunctionCreateInteriorByName(ceed, "Mass1DBuild", &qf_setup); 67 CeedQFunctionCreateInteriorByName(ceed, "MassApply", &qf_mass); 68 69 // Operators 70 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 71 &op_setup); 72 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 73 &op_mass_f); 74 75 CeedVectorCreate(ceed, num_dofs_x, &X); 76 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 77 CeedVectorCreate(ceed, num_elem*Q, &q_data); 78 79 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, 80 CEED_VECTOR_NONE); 81 CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 82 CeedOperatorSetField(op_setup, "qdata", elem_restr_qd, CEED_BASIS_COLLOCATED, 83 CEED_VECTOR_ACTIVE); 84 85 CeedOperatorSetField(op_mass_f, "qdata", elem_restr_qd, CEED_BASIS_COLLOCATED, 86 q_data); 87 CeedOperatorSetField(op_mass_f, "u", elem_restr_u_f, basis_u, 88 CEED_VECTOR_ACTIVE); 89 CeedOperatorSetField(op_mass_f, "v", elem_restr_u_f, basis_u, 90 CEED_VECTOR_ACTIVE); 91 92 CeedOperatorApply(op_setup, X, q_data, CEED_REQUEST_IMMEDIATE); 93 94 // Create multigrid level 95 CeedVectorCreate(ceed, num_dofs_u_f, &p_mult_f); 96 CeedVectorSetValue(p_mult_f, 1.0); 97 CeedBasis basis_u_c, basis_c_to_f; 98 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P_c, Q, CEED_GAUSS, &basis_u_c); 99 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P_c, P_f, CEED_GAUSS_LOBATTO, 100 &basis_c_to_f); 101 const CeedScalar *interp_c_to_f; 102 CeedBasisGetInterp1D(basis_c_to_f, &interp_c_to_f); 103 CeedOperatorMultigridLevelCreateH1(op_mass_f, p_mult_f, elem_restr_u_c, 104 basis_u_c, interp_c_to_f, &op_mass_c, 105 &op_prolong, &op_restrict); 106 107 // Coarse problem 108 CeedVectorCreate(ceed, num_dofs_u_c, &U_c); 109 CeedVectorSetValue(U_c, 1.0); 110 CeedVectorCreate(ceed, num_dofs_u_c, &V_c); 111 CeedOperatorApply(op_mass_c, U_c, V_c, CEED_REQUEST_IMMEDIATE); 112 113 // Check output 114 CeedVectorGetArrayRead(V_c, CEED_MEM_HOST, &hv); 115 sum = 0.; 116 for (CeedInt i=0; i<num_dofs_u_c; i++) { 117 sum += hv[i]; 118 } 119 if (fabs(sum-1.)>1e-10) 120 // LCOV_EXCL_START 121 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 122 // LCOV_EXCL_STOP 123 CeedVectorRestoreArrayRead(V_c, &hv); 124 125 // Prolong coarse u 126 CeedVectorCreate(ceed, num_dofs_u_f, &U_f); 127 CeedOperatorApply(op_prolong, U_c, U_f, CEED_REQUEST_IMMEDIATE); 128 129 // Fine problem 130 CeedVectorCreate(ceed, num_dofs_u_f, &V_f); 131 CeedOperatorApply(op_mass_f, U_f, V_f, CEED_REQUEST_IMMEDIATE); 132 133 // Check output 134 CeedVectorGetArrayRead(V_f, CEED_MEM_HOST, &hv); 135 sum = 0.; 136 for (CeedInt i=0; i<num_dofs_u_f; i++) { 137 sum += hv[i]; 138 } 139 if (fabs(sum-1.)>1e-10) 140 // LCOV_EXCL_START 141 printf("Computed Area Fine Grid: %f != True Area: 1.0\n", sum); 142 // LCOV_EXCL_STOP 143 CeedVectorRestoreArrayRead(V_f, &hv); 144 145 // Restrict state to coarse grid 146 CeedOperatorApply(op_restrict, V_f, V_c, CEED_REQUEST_IMMEDIATE); 147 148 // Check output 149 CeedVectorGetArrayRead(V_c, CEED_MEM_HOST, &hv); 150 sum = 0.; 151 for (CeedInt i=0; i<num_dofs_u_c; i++) { 152 sum += hv[i]; 153 } 154 if (fabs(sum-1.)>1e-10) 155 // LCOV_EXCL_START 156 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 157 // LCOV_EXCL_STOP 158 CeedVectorRestoreArrayRead(V_c, &hv); 159 160 // Cleanup 161 CeedQFunctionDestroy(&qf_setup); 162 CeedQFunctionDestroy(&qf_mass); 163 CeedOperatorDestroy(&op_setup); 164 CeedOperatorDestroy(&op_mass_c); 165 CeedOperatorDestroy(&op_mass_f); 166 CeedOperatorDestroy(&op_prolong); 167 CeedOperatorDestroy(&op_restrict); 168 CeedElemRestrictionDestroy(&elem_restr_u_c); 169 CeedElemRestrictionDestroy(&elem_restr_u_f); 170 CeedElemRestrictionDestroy(&elem_restr_x); 171 CeedElemRestrictionDestroy(&elem_restr_qd); 172 CeedBasisDestroy(&basis_u); 173 CeedBasisDestroy(&basis_u_c); 174 CeedBasisDestroy(&basis_c_to_f); 175 CeedBasisDestroy(&basis_x); 176 CeedVectorDestroy(&X); 177 CeedVectorDestroy(&U_c); 178 CeedVectorDestroy(&U_f); 179 CeedVectorDestroy(&V_c); 180 CeedVectorDestroy(&V_f); 181 CeedVectorDestroy(&p_mult_f); 182 CeedVectorDestroy(&q_data); 183 CeedDestroy(&ceed); 184 return 0; 185 } 186