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