1 /// @file 2 /// Test FLOP estimation for composite mass matrix operator 3 /// \test Test FLOP estimation for composite mass matrix operator 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t320-basis.h" 8 9 /* The mesh comprises of two rows of 3 quadralaterals followed by one row 10 of 6 triangles: 11 _ _ _ 12 |_|_|_| 13 |_|_|_| 14 |/|/|/| 15 16 */ 17 18 int main(int argc, char **argv) { 19 Ceed ceed; 20 CeedSize flop_estimate; 21 CeedElemRestriction elem_restr_x_tet, elem_restr_u_tet, 22 elem_restr_qd_i_tet, 23 elem_restr_x_hex, elem_restr_u_hex, 24 elem_restr_qd_i_hex; 25 CeedBasis basis_x_tet, basis_u_tet, 26 basis_x_hex, basis_u_hex; 27 CeedQFunction qf_mass; 28 CeedOperator op_mass_tet, op_mass_hex, op_mass; 29 CeedVector q_data_tet, q_data_hex; 30 CeedInt num_elem_tet = 6, P_tet = 6, Q_tet = 4, 31 num_elem_hex = 6, P_hex = 3, Q_hex = 4, dim = 2; 32 CeedInt n_x = 3, n_y = 3, 33 n_x_tet = 3, n_y_tet = 1, n_x_hex = 3; 34 CeedInt row, col, offset; 35 CeedInt num_dofs = (n_x*2+1)*(n_y*2+1), 36 num_qpts_tet = num_elem_tet*Q_tet, 37 num_qpts_hex = num_elem_hex*Q_hex*Q_hex; 38 CeedInt ind_x_tet[num_elem_tet*P_tet], 39 ind_x_hex[num_elem_hex*P_hex*P_hex]; 40 CeedScalar q_ref[dim*Q_tet], q_weight[Q_tet]; 41 CeedScalar interp[P_tet*Q_tet], grad[dim*P_tet*Q_tet]; 42 43 CeedInit(argv[1], &ceed); 44 45 // Qdata Vectors 46 CeedVectorCreate(ceed, num_qpts_tet, &q_data_tet); 47 CeedVectorCreate(ceed, num_qpts_hex, &q_data_hex); 48 49 // Set up Tet Elements 50 for (CeedInt i=0; i<num_elem_tet/2; i++) { 51 col = i % n_x_tet; 52 row = i / n_x_tet; 53 offset = col*2 + row*(n_x_tet*2+1)*2; 54 55 ind_x_tet[i*2*P_tet + 0] = 2 + offset; 56 ind_x_tet[i*2*P_tet + 1] = 9 + offset; 57 ind_x_tet[i*2*P_tet + 2] = 16 + offset; 58 ind_x_tet[i*2*P_tet + 3] = 1 + offset; 59 ind_x_tet[i*2*P_tet + 4] = 8 + offset; 60 ind_x_tet[i*2*P_tet + 5] = 0 + offset; 61 62 ind_x_tet[i*2*P_tet + 6] = 14 + offset; 63 ind_x_tet[i*2*P_tet + 7] = 7 + offset; 64 ind_x_tet[i*2*P_tet + 8] = 0 + offset; 65 ind_x_tet[i*2*P_tet + 9] = 15 + offset; 66 ind_x_tet[i*2*P_tet + 10] = 8 + offset; 67 ind_x_tet[i*2*P_tet + 11] = 16 + offset; 68 } 69 70 // -- Restrictions 71 CeedElemRestrictionCreate(ceed, num_elem_tet, P_tet, dim, num_dofs, 72 dim*num_dofs, 73 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 74 &elem_restr_x_tet); 75 76 CeedElemRestrictionCreate(ceed, num_elem_tet, P_tet, 1, 1, num_dofs, 77 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 78 &elem_restr_u_tet); 79 CeedInt strides_qd_tet[3] = {1, Q_tet, Q_tet}; 80 CeedElemRestrictionCreateStrided(ceed, num_elem_tet, Q_tet, 1, num_qpts_tet, 81 strides_qd_tet, &elem_restr_qd_i_tet); 82 83 // -- Bases 84 buildmats(q_ref, q_weight, interp, grad); 85 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, dim, P_tet, Q_tet, interp, grad, 86 q_ref, q_weight, &basis_x_tet); 87 88 buildmats(q_ref, q_weight, interp, grad); 89 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, P_tet, Q_tet, interp, grad, 90 q_ref, q_weight, &basis_u_tet); 91 92 // -- QFunction 93 CeedQFunctionCreateInteriorByName(ceed, "MassApply", &qf_mass); 94 95 // -- Operators 96 // ---- Mass Tet 97 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 98 &op_mass_tet); 99 CeedOperatorSetField(op_mass_tet, "u", elem_restr_u_tet, basis_u_tet, 100 CEED_VECTOR_ACTIVE); 101 CeedOperatorSetField(op_mass_tet, "qdata", elem_restr_qd_i_tet, 102 CEED_BASIS_COLLOCATED, 103 q_data_tet); 104 CeedOperatorSetField(op_mass_tet, "v", elem_restr_u_tet, basis_u_tet, 105 CEED_VECTOR_ACTIVE); 106 107 // Set up Hex Elements 108 for (CeedInt i=0; i<num_elem_hex; i++) { 109 col = i % n_x_hex; 110 row = i / n_x_hex; 111 offset = (n_x_tet*2+1)*(n_y_tet*2)*(1+row) + col*2; 112 for (CeedInt j=0; j<P_hex; j++) 113 for (CeedInt k=0; k<P_hex; k++) 114 ind_x_hex[P_hex*(P_hex*i+k)+j] = offset + k*(n_x_hex*2+1) + j; 115 } 116 117 // -- Restrictions 118 CeedElemRestrictionCreate(ceed, num_elem_hex, P_hex*P_hex, dim, num_dofs, 119 dim*num_dofs, 120 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 121 &elem_restr_x_hex); 122 123 CeedElemRestrictionCreate(ceed, num_elem_hex, P_hex*P_hex, 1, 1, num_dofs, 124 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 125 &elem_restr_u_hex); 126 CeedInt strides_qd_hex[3] = {1, Q_hex*Q_hex, Q_hex*Q_hex}; 127 CeedElemRestrictionCreateStrided(ceed, num_elem_hex, Q_hex*Q_hex, 1, 128 num_qpts_hex, 129 strides_qd_hex, &elem_restr_qd_i_hex); 130 131 // -- Bases 132 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, P_hex, Q_hex, CEED_GAUSS, 133 &basis_x_hex); 134 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_hex, Q_hex, CEED_GAUSS, 135 &basis_u_hex); 136 137 // -- Operators 138 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 139 &op_mass_hex); 140 CeedOperatorSetField(op_mass_hex, "u", elem_restr_u_hex, basis_u_hex, 141 CEED_VECTOR_ACTIVE); 142 CeedOperatorSetField(op_mass_hex, "qdata", elem_restr_qd_i_hex, 143 CEED_BASIS_COLLOCATED, 144 q_data_hex); 145 CeedOperatorSetField(op_mass_hex, "v", elem_restr_u_hex, basis_u_hex, 146 CEED_VECTOR_ACTIVE); 147 148 // Set up Composite Operator 149 // -- Create 150 CeedCompositeOperatorCreate(ceed, &op_mass); 151 // -- Add SubOperators 152 CeedCompositeOperatorAddSub(op_mass, op_mass_tet); 153 CeedCompositeOperatorAddSub(op_mass, op_mass_hex); 154 155 // Estimate FLOPs 156 CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); 157 CeedOperatorGetFlopsEstimate(op_mass, &flop_estimate); 158 159 // Check output 160 if (flop_estimate != 3042) 161 // LCOV_EXCL_START 162 printf("Incorrect FLOP estimate computed, %ld != 3042\n", flop_estimate); 163 // LCOV_EXCL_STOP 164 165 // Cleanup 166 CeedQFunctionDestroy(&qf_mass); 167 CeedOperatorDestroy(&op_mass_tet); 168 CeedOperatorDestroy(&op_mass_hex); 169 CeedOperatorDestroy(&op_mass); 170 CeedElemRestrictionDestroy(&elem_restr_u_tet); 171 CeedElemRestrictionDestroy(&elem_restr_x_tet); 172 CeedElemRestrictionDestroy(&elem_restr_qd_i_tet); 173 CeedElemRestrictionDestroy(&elem_restr_u_hex); 174 CeedElemRestrictionDestroy(&elem_restr_x_hex); 175 CeedElemRestrictionDestroy(&elem_restr_qd_i_hex); 176 CeedBasisDestroy(&basis_u_tet); 177 CeedBasisDestroy(&basis_x_tet); 178 CeedBasisDestroy(&basis_u_hex); 179 CeedBasisDestroy(&basis_x_hex); 180 CeedVectorDestroy(&q_data_tet); 181 CeedVectorDestroy(&q_data_hex); 182 CeedDestroy(&ceed); 183 return 0; 184 } 185