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