1 /// @file 2 /// Test full assembly of mass matrix operator 3 /// \test Test full assembly of mass matrix operator AtPoints 4 #include <ceed.h> 5 #include <math.h> 6 #include <stdio.h> 7 #include <stdlib.h> 8 9 #include "t596-operator.h" 10 11 int main(int argc, char **argv) { 12 Ceed ceed; 13 14 CeedInit(argv[1], &ceed); 15 16 for (CeedInt num_comp = 1; num_comp <= 3; num_comp++) { 17 CeedElemRestriction elem_restriction_x, elem_restriction_x_points, elem_restriction_u, elem_restriction_q_data; 18 CeedBasis basis_x, basis_u; 19 CeedQFunction qf_setup, qf_mass; 20 CeedOperator op_setup, op_mass; 21 CeedVector q_data, x, x_points, u, v; 22 CeedInt p = 3, q = 4, dim = 2; 23 CeedInt n_x = 3, n_y = 2; 24 CeedInt num_elem = n_x * n_y; 25 CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; 26 CeedInt ind_x[num_elem * p * p]; 27 CeedScalar assembled_values[num_comp * num_comp * num_dofs * num_dofs]; 28 CeedScalar assembled_true[num_comp * num_comp * num_dofs * num_dofs]; 29 30 // Points 31 CeedVectorCreate(ceed, dim * num_points, &x_points); 32 { 33 CeedScalar x_array[dim * num_points]; 34 35 for (CeedInt e = 0; e < num_elem; e++) { 36 for (CeedInt d = 0; d < dim; d++) { 37 x_array[num_points_per_elem * (e * dim + d) + 0] = 0.25; 38 x_array[num_points_per_elem * (e * dim + d) + 1] = d == 0 ? -0.25 : 0.25; 39 x_array[num_points_per_elem * (e * dim + d) + 2] = d == 0 ? 0.25 : -0.25; 40 x_array[num_points_per_elem * (e * dim + d) + 3] = 0.25; 41 } 42 } 43 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 44 } 45 { 46 CeedInt ind_x[num_elem + 1 + num_points]; 47 48 for (CeedInt i = 0; i <= num_elem; i++) ind_x[i] = num_elem + 1 + i * num_points_per_elem; 49 for (CeedInt i = 0; i < num_points; i++) ind_x[num_elem + 1 + i] = i; 50 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, 51 &elem_restriction_x_points); 52 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_q_data); 53 } 54 55 // Vectors 56 CeedVectorCreate(ceed, dim * num_dofs, &x); 57 { 58 CeedScalar x_array[dim * num_dofs]; 59 60 for (CeedInt i = 0; i < n_x * 2 + 1; i++) { 61 for (CeedInt j = 0; j < n_y * 2 + 1; j++) { 62 x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x); 63 x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y); 64 } 65 } 66 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 67 } 68 CeedVectorCreate(ceed, num_comp * num_dofs, &u); 69 CeedVectorCreate(ceed, num_comp * num_dofs, &v); 70 CeedVectorCreate(ceed, num_points, &q_data); 71 72 // Restrictions 73 for (CeedInt i = 0; i < num_elem; i++) { 74 CeedInt col, row, offset; 75 col = i % n_x; 76 row = i / n_x; 77 offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1); 78 for (CeedInt j = 0; j < p; j++) { 79 for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j; 80 } 81 } 82 CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); 83 CeedElemRestrictionCreate(ceed, num_elem, p * p, num_comp, num_dofs, num_comp * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, 84 &elem_restriction_u); 85 86 // Bases 87 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x); 88 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u); 89 90 // QFunctions 91 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 92 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 93 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); 94 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 95 96 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 97 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); 98 CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP); 99 CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP); 100 { 101 CeedQFunctionContext qf_context; 102 103 CeedQFunctionContextCreate(ceed, &qf_context); 104 CeedQFunctionContextSetData(qf_context, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(CeedInt), &num_comp); 105 CeedQFunctionSetContext(qf_mass, qf_context); 106 CeedQFunctionContextDestroy(&qf_context); 107 } 108 109 // Operators 110 CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 111 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 112 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 113 CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 114 CeedOperatorAtPointsSetPoints(op_setup, elem_restriction_x_points, x_points); 115 116 CeedOperatorCreateAtPoints(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); 117 CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); 118 CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 119 CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 120 CeedOperatorAtPointsSetPoints(op_mass, elem_restriction_x_points, x_points); 121 122 // Apply Setup Operator 123 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); 124 125 // Fully assemble operator 126 CeedSize num_entries; 127 CeedInt *rows; 128 CeedInt *cols; 129 CeedVector assembled; 130 131 for (CeedInt k = 0; k < num_comp * num_comp * num_dofs * num_dofs; ++k) { 132 assembled_values[k] = 0.0; 133 assembled_true[k] = 0.0; 134 } 135 CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols); 136 CeedVectorCreate(ceed, num_entries, &assembled); 137 CeedOperatorLinearAssemble(op_mass, assembled); 138 { 139 const CeedScalar *assembled_array; 140 141 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); 142 for (CeedInt k = 0; k < num_entries; k++) { 143 assembled_values[rows[k] * num_comp * num_dofs + cols[k]] += assembled_array[k]; 144 } 145 CeedVectorRestoreArrayRead(assembled, &assembled_array); 146 } 147 148 // Manually assemble operator 149 CeedVectorSetValue(u, 0.0); 150 for (CeedInt j = 0; j < num_comp * num_dofs; j++) { 151 CeedScalar *u_array; 152 const CeedScalar *v_array; 153 154 // Set input 155 CeedVectorGetArray(u, CEED_MEM_HOST, &u_array); 156 u_array[j] = 1.0; 157 if (j) u_array[j - 1] = 0.0; 158 CeedVectorRestoreArray(u, &u_array); 159 160 // Compute entries for column j 161 CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE); 162 163 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 164 for (CeedInt i = 0; i < num_comp * num_dofs; i++) assembled_true[i * num_comp * num_dofs + j] = v_array[i]; 165 CeedVectorRestoreArrayRead(v, &v_array); 166 } 167 168 // Check output 169 for (CeedInt i = 0; i < num_comp * num_dofs; i++) { 170 for (CeedInt j = 0; j < num_comp * num_dofs; j++) { 171 if (fabs(assembled_values[i * num_dofs * num_comp + j] - assembled_true[i * num_dofs * num_comp + j]) > 100. * CEED_EPSILON) { 172 // LCOV_EXCL_START 173 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_dofs * num_comp + j], 174 assembled_true[i * num_dofs * num_comp + j]); 175 // LCOV_EXCL_STOP 176 } 177 } 178 } 179 180 // Cleanup 181 free(rows); 182 free(cols); 183 CeedVectorDestroy(&x); 184 CeedVectorDestroy(&x_points); 185 CeedVectorDestroy(&q_data); 186 CeedVectorDestroy(&u); 187 CeedVectorDestroy(&v); 188 CeedVectorDestroy(&assembled); 189 CeedElemRestrictionDestroy(&elem_restriction_u); 190 CeedElemRestrictionDestroy(&elem_restriction_x); 191 CeedElemRestrictionDestroy(&elem_restriction_x_points); 192 CeedElemRestrictionDestroy(&elem_restriction_q_data); 193 CeedBasisDestroy(&basis_u); 194 CeedBasisDestroy(&basis_x); 195 CeedQFunctionDestroy(&qf_setup); 196 CeedQFunctionDestroy(&qf_mass); 197 CeedOperatorDestroy(&op_setup); 198 CeedOperatorDestroy(&op_mass); 199 } 200 CeedDestroy(&ceed); 201 return 0; 202 } 203