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_coarse, op_mass_fine, op_prolong, op_restrict; 14 CeedVector x, u_coarse, u_fine, v_coarse, v_fine, p_mult_fine; 15 CeedInt num_elem = 15, num_elem_sub = 5, num_sub_ops = 3, p_coarse = 3, p_fine = 5, q = 8, num_comp = 2; 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_sub * p_coarse], ind_u_fine[num_elem_sub * p_fine], ind_x[num_elem_sub * 2]; 18 19 CeedInit(argv[1], &ceed); 20 21 // Vectors 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_comp * num_dofs_u_coarse, &u_coarse); 30 CeedVectorCreate(ceed, num_comp * num_dofs_u_coarse, &v_coarse); 31 CeedVectorCreate(ceed, num_comp * num_dofs_u_fine, &u_fine); 32 CeedVectorCreate(ceed, num_comp * num_dofs_u_fine, &v_fine); 33 34 // Composite operators 35 CeedCompositeOperatorCreate(ceed, &op_mass_coarse); 36 CeedCompositeOperatorCreate(ceed, &op_mass_fine); 37 CeedCompositeOperatorCreate(ceed, &op_prolong); 38 CeedCompositeOperatorCreate(ceed, &op_restrict); 39 40 // Setup fine suboperators 41 for (CeedInt i = 0; i < num_sub_ops; i++) { 42 CeedVector q_data; 43 CeedElemRestriction elem_restriction_x, elem_restriction_q_data, elem_restriction_u_fine; 44 CeedBasis basis_x, basis_u_fine; 45 CeedQFunction qf_setup, qf_mass; 46 CeedOperator sub_op_setup, sub_op_mass_fine; 47 48 // -- QData 49 CeedVectorCreate(ceed, num_elem * q, &q_data); 50 51 // -- Restrictions 52 CeedInt offset = num_elem_sub * i; 53 for (CeedInt j = 0; j < num_elem_sub; j++) { 54 ind_x[2 * j + 0] = offset + j; 55 ind_x[2 * j + 1] = offset + j + 1; 56 } 57 CeedElemRestrictionCreate(ceed, num_elem_sub, 2, 1, 1, num_dofs_x, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_x); 58 59 offset = num_elem_sub * i * (p_fine - 1); 60 for (CeedInt j = 0; j < num_elem_sub; j++) { 61 for (CeedInt k = 0; k < p_fine; k++) { 62 ind_u_fine[p_fine * j + k] = offset + j * (p_fine - 1) + k; 63 } 64 } 65 CeedElemRestrictionCreate(ceed, num_elem_sub, p_fine, num_comp, num_dofs_u_fine, num_comp * num_dofs_u_fine, CEED_MEM_HOST, CEED_COPY_VALUES, 66 ind_u_fine, &elem_restriction_u_fine); 67 68 CeedInt strides_q_data[3] = {1, q, q}; 69 CeedElemRestrictionCreateStrided(ceed, num_elem_sub, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data); 70 71 // -- Bases 72 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x); 73 CeedBasisCreateTensorH1Lagrange(ceed, 1, num_comp, p_fine, q, CEED_GAUSS, &basis_u_fine); 74 75 // -- QFunctions 76 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 77 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 78 CeedQFunctionAddInput(qf_setup, "dx", 1 * 1, CEED_EVAL_GRAD); 79 CeedQFunctionAddOutput(qf_setup, "q data", 1, CEED_EVAL_NONE); 80 81 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 82 CeedQFunctionAddInput(qf_mass, "q data", 1, CEED_EVAL_NONE); 83 CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP); 84 CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP); 85 86 // -- Operators 87 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &sub_op_setup); 88 CeedOperatorSetField(sub_op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 89 CeedOperatorSetField(sub_op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 90 CeedOperatorSetField(sub_op_setup, "q data", elem_restriction_q_data, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 91 92 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &sub_op_mass_fine); 93 CeedOperatorSetField(sub_op_mass_fine, "q data", elem_restriction_q_data, CEED_BASIS_COLLOCATED, q_data); 94 CeedOperatorSetField(sub_op_mass_fine, "u", elem_restriction_u_fine, basis_u_fine, CEED_VECTOR_ACTIVE); 95 CeedOperatorSetField(sub_op_mass_fine, "v", elem_restriction_u_fine, basis_u_fine, CEED_VECTOR_ACTIVE); 96 97 // -- Create qdata 98 CeedOperatorApply(sub_op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); 99 100 // -- Composite operators 101 CeedCompositeOperatorAddSub(op_mass_fine, sub_op_mass_fine); 102 103 // -- Cleanup 104 CeedVectorDestroy(&q_data); 105 CeedElemRestrictionDestroy(&elem_restriction_u_fine); 106 CeedElemRestrictionDestroy(&elem_restriction_x); 107 CeedElemRestrictionDestroy(&elem_restriction_q_data); 108 CeedBasisDestroy(&basis_u_fine); 109 CeedBasisDestroy(&basis_x); 110 CeedQFunctionDestroy(&qf_setup); 111 CeedQFunctionDestroy(&qf_mass); 112 CeedOperatorDestroy(&sub_op_setup); 113 CeedOperatorDestroy(&sub_op_mass_fine); 114 } 115 116 // Scale for suboperator multiplicity 117 CeedVectorCreate(ceed, num_comp * num_dofs_u_fine, &p_mult_fine); 118 CeedCompositeOperatorGetMultiplicity(op_mass_fine, 0, NULL, p_mult_fine); 119 120 // Setup coarse and prolong/restriction suboperators 121 for (CeedInt i = 0; i < num_sub_ops; i++) { 122 CeedElemRestriction elem_restriction_u_coarse; 123 CeedBasis basis_u_coarse; 124 CeedOperator *sub_ops_mass_fine, sub_op_mass_coarse, sub_op_prolong, sub_op_restrict; 125 126 // -- Fine grid operator 127 CeedCompositeOperatorGetSubList(op_mass_fine, &sub_ops_mass_fine); 128 129 // -- Restrictions 130 CeedInt offset = num_elem_sub * i * (p_coarse - 1); 131 for (CeedInt j = 0; j < num_elem_sub; j++) { 132 for (CeedInt k = 0; k < p_coarse; k++) { 133 ind_u_coarse[p_coarse * j + k] = offset + j * (p_coarse - 1) + k; 134 } 135 } 136 CeedElemRestrictionCreate(ceed, num_elem_sub, p_coarse, num_comp, num_dofs_u_coarse, num_comp * num_dofs_u_coarse, CEED_MEM_HOST, 137 CEED_COPY_VALUES, ind_u_coarse, &elem_restriction_u_coarse); 138 139 // -- Bases 140 CeedBasisCreateTensorH1Lagrange(ceed, 1, num_comp, p_coarse, q, CEED_GAUSS, &basis_u_coarse); 141 142 // -- Create multigrid level 143 CeedOperatorMultigridLevelCreate(sub_ops_mass_fine[i], p_mult_fine, elem_restriction_u_coarse, basis_u_coarse, &sub_op_mass_coarse, 144 &sub_op_prolong, &sub_op_restrict); 145 146 // -- Composite operators 147 CeedCompositeOperatorAddSub(op_mass_coarse, sub_op_mass_coarse); 148 CeedCompositeOperatorAddSub(op_prolong, sub_op_prolong); 149 CeedCompositeOperatorAddSub(op_restrict, sub_op_restrict); 150 151 // -- Cleanup 152 CeedElemRestrictionDestroy(&elem_restriction_u_coarse); 153 CeedBasisDestroy(&basis_u_coarse); 154 CeedOperatorDestroy(&sub_op_mass_coarse); 155 CeedOperatorDestroy(&sub_op_prolong); 156 CeedOperatorDestroy(&sub_op_restrict); 157 } 158 159 // Coarse problem 160 CeedVectorSetValue(u_coarse, 1.0); 161 CeedOperatorApply(op_mass_coarse, u_coarse, v_coarse, CEED_REQUEST_IMMEDIATE); 162 163 // Check output 164 { 165 const CeedScalar *v_array; 166 CeedScalar sum = 0.; 167 168 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 169 for (CeedInt i = 0; i < num_comp * num_dofs_u_coarse; i++) { 170 sum += v_array[i]; 171 } 172 CeedVectorRestoreArrayRead(v_coarse, &v_array); 173 if (fabs(sum - 2.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 2.0\n", sum); 174 } 175 176 // Prolong coarse u 177 CeedOperatorApply(op_prolong, u_coarse, u_fine, CEED_REQUEST_IMMEDIATE); 178 179 // Fine problem 180 CeedOperatorApply(op_mass_fine, u_fine, v_fine, CEED_REQUEST_IMMEDIATE); 181 182 // Check output 183 { 184 const CeedScalar *v_array; 185 CeedScalar sum = 0.; 186 CeedVectorGetArrayRead(v_fine, CEED_MEM_HOST, &v_array); 187 for (CeedInt i = 0; i < num_comp * num_dofs_u_fine; i++) { 188 sum += v_array[i]; 189 } 190 CeedVectorRestoreArrayRead(v_fine, &v_array); 191 192 if (fabs(sum - 2.) > 1000. * CEED_EPSILON) printf("Computed Area Fine Grid: %f != True Area: 2.0\n", sum); 193 } 194 // Restrict state to coarse grid 195 CeedOperatorApply(op_restrict, v_fine, v_coarse, CEED_REQUEST_IMMEDIATE); 196 197 // Check output 198 { 199 const CeedScalar *v_array; 200 CeedScalar sum = 0.; 201 202 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 203 for (CeedInt i = 0; i < num_comp * num_dofs_u_coarse; i++) { 204 sum += v_array[i]; 205 } 206 CeedVectorRestoreArrayRead(v_coarse, &v_array); 207 if (fabs(sum - 2.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 2.0\n", sum); 208 } 209 210 // Cleanup 211 CeedVectorDestroy(&x); 212 CeedVectorDestroy(&u_coarse); 213 CeedVectorDestroy(&u_fine); 214 CeedVectorDestroy(&v_coarse); 215 CeedVectorDestroy(&v_fine); 216 CeedVectorDestroy(&p_mult_fine); 217 CeedOperatorDestroy(&op_mass_coarse); 218 CeedOperatorDestroy(&op_mass_fine); 219 CeedOperatorDestroy(&op_prolong); 220 CeedOperatorDestroy(&op_restrict); 221 CeedDestroy(&ceed); 222 return 0; 223 } 224