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