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