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_restriction_x, elem_restriction_q_data, elem_restriction_u_coarse, elem_restriction_u_fine; 11 CeedBasis basis_x, basis_u; 12 CeedQFunction qf_setup, qf_mass; 13 CeedOperator op_setup, op_mass_coarse, op_mass_fine, op_prolong, op_restrict; 14 CeedVector q_data, x, u_coarse, u_fine, v_coarse, v_fine, p_mult_fine; 15 CeedInt num_elem = 15, p_coarse = 3, p_fine = 5, q = 8; 16 CeedInt num_dofs_x = num_elem + 1, num_dofs_u_coarse = num_elem * (p_coarse - 1) + 1, num_dofs_u_fine = num_elem * (p_fine - 1) + 1; 17 CeedInt ind_u_coarse[num_elem * p_coarse], ind_u_fine[num_elem * p_fine], ind_x[num_elem * 2]; 18 19 CeedInit(argv[1], &ceed); 20 21 CeedVectorCreate(ceed, num_dofs_x, &x); 22 { 23 CeedScalar x_array[num_dofs_x]; 24 25 for (CeedInt i = 0; i < num_dofs_x; i++) x_array[i] = (CeedScalar)i / (num_dofs_x - 1); 26 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 27 } 28 CeedVectorCreate(ceed, num_dofs_u_fine, &p_mult_fine); 29 CeedVectorCreate(ceed, num_dofs_u_coarse, &u_coarse); 30 CeedVectorCreate(ceed, num_dofs_u_coarse, &v_coarse); 31 CeedVectorCreate(ceed, num_dofs_u_fine, &u_fine); 32 CeedVectorCreate(ceed, num_dofs_u_fine, &v_fine); 33 CeedVectorCreate(ceed, num_elem * q, &q_data); 34 35 // Restrictions 36 for (CeedInt i = 0; i < num_elem; i++) { 37 ind_x[2 * i + 0] = i; 38 ind_x[2 * i + 1] = i + 1; 39 } 40 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_dofs_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); 41 42 for (CeedInt i = 0; i < num_elem; i++) { 43 for (CeedInt j = 0; j < p_coarse; j++) { 44 ind_u_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j; 45 } 46 } 47 CeedElemRestrictionCreate(ceed, num_elem, p_coarse, 1, 1, num_dofs_u_coarse, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_coarse, 48 &elem_restriction_u_coarse); 49 50 for (CeedInt i = 0; i < num_elem; i++) { 51 for (CeedInt j = 0; j < p_fine; j++) { 52 ind_u_fine[p_fine * i + j] = i * (p_fine - 1) + j; 53 } 54 } 55 CeedElemRestrictionCreate(ceed, num_elem, p_fine, 1, 1, num_dofs_u_fine, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_fine, &elem_restriction_u_fine); 56 57 CeedInt strides_q_data[3] = {1, q, q}; 58 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data); 59 60 // Bases 61 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x); 62 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_fine, q, CEED_GAUSS, &basis_u); 63 64 // QFunctions 65 CeedQFunctionCreateInteriorByName(ceed, "Mass1DBuild", &qf_setup); 66 CeedQFunctionCreateInteriorByName(ceed, "MassApply", &qf_mass); 67 68 // Operators 69 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 70 71 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 72 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 73 CeedOperatorSetField(op_setup, "qdata", elem_restriction_q_data, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 74 75 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_fine); 76 CeedOperatorSetField(op_mass_fine, "qdata", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data); 77 CeedOperatorSetField(op_mass_fine, "u", elem_restriction_u_fine, basis_u, CEED_VECTOR_ACTIVE); 78 CeedOperatorSetField(op_mass_fine, "v", elem_restriction_u_fine, basis_u, CEED_VECTOR_ACTIVE); 79 80 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); 81 82 // Create multigrid level 83 CeedVectorSetValue(p_mult_fine, 1.0); 84 CeedBasis basis_u_coarse, basis_coarse_to_fine; 85 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_coarse, q, CEED_GAUSS, &basis_u_coarse); 86 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_coarse, p_fine, CEED_GAUSS_LOBATTO, &basis_coarse_to_fine); 87 const CeedScalar *interp_coarse_to_fine; 88 CeedBasisGetInterp1D(basis_coarse_to_fine, &interp_coarse_to_fine); 89 CeedOperatorMultigridLevelCreateH1(op_mass_fine, p_mult_fine, elem_restriction_u_coarse, basis_u_coarse, interp_coarse_to_fine, &op_mass_coarse, 90 &op_prolong, &op_restrict); 91 92 // Coarse problem 93 CeedVectorSetValue(u_coarse, 1.0); 94 CeedOperatorApply(op_mass_coarse, u_coarse, v_coarse, CEED_REQUEST_IMMEDIATE); 95 96 // Check output 97 { 98 const CeedScalar *v_array; 99 CeedScalar sum = 0.; 100 101 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 102 for (CeedInt i = 0; i < num_dofs_u_coarse; i++) { 103 sum += v_array[i]; 104 } 105 CeedVectorRestoreArrayRead(v_coarse, &v_array); 106 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 107 } 108 109 // Prolong coarse u 110 CeedOperatorApply(op_prolong, u_coarse, u_fine, CEED_REQUEST_IMMEDIATE); 111 112 // Fine problem 113 CeedOperatorApply(op_mass_fine, u_fine, v_fine, CEED_REQUEST_IMMEDIATE); 114 115 // Check output 116 { 117 const CeedScalar *v_array; 118 CeedScalar sum = 0.; 119 120 CeedVectorGetArrayRead(v_fine, CEED_MEM_HOST, &v_array); 121 for (CeedInt i = 0; i < num_dofs_u_fine; i++) { 122 sum += v_array[i]; 123 } 124 CeedVectorRestoreArrayRead(v_fine, &v_array); 125 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Fine Grid: %f != True Area: 1.0\n", sum); 126 } 127 128 // Restrict state to coarse grid 129 CeedOperatorApply(op_restrict, v_fine, v_coarse, CEED_REQUEST_IMMEDIATE); 130 131 // Check output 132 { 133 const CeedScalar *v_array; 134 CeedScalar sum = 0.; 135 136 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 137 for (CeedInt i = 0; i < num_dofs_u_coarse; i++) { 138 sum += v_array[i]; 139 } 140 CeedVectorRestoreArrayRead(v_coarse, &v_array); 141 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 142 } 143 144 // Cleanup 145 CeedVectorDestroy(&x); 146 CeedVectorDestroy(&u_coarse); 147 CeedVectorDestroy(&u_fine); 148 CeedVectorDestroy(&v_coarse); 149 CeedVectorDestroy(&v_fine); 150 CeedVectorDestroy(&p_mult_fine); 151 CeedVectorDestroy(&q_data); 152 CeedElemRestrictionDestroy(&elem_restriction_u_coarse); 153 CeedElemRestrictionDestroy(&elem_restriction_u_fine); 154 CeedElemRestrictionDestroy(&elem_restriction_x); 155 CeedElemRestrictionDestroy(&elem_restriction_q_data); 156 CeedBasisDestroy(&basis_u); 157 CeedBasisDestroy(&basis_u_coarse); 158 CeedBasisDestroy(&basis_coarse_to_fine); 159 CeedBasisDestroy(&basis_x); 160 CeedQFunctionDestroy(&qf_setup); 161 CeedQFunctionDestroy(&qf_mass); 162 CeedOperatorDestroy(&op_setup); 163 CeedOperatorDestroy(&op_mass_coarse); 164 CeedOperatorDestroy(&op_mass_fine); 165 CeedOperatorDestroy(&op_prolong); 166 CeedOperatorDestroy(&op_restrict); 167 CeedDestroy(&ceed); 168 return 0; 169 } 170