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