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