1 /// @file 2 /// Test creation, action, and destruction for mass matrix operator AtPoints 3 /// \test Test creation, action, and destruction for mass matrix operator AtPoints 4 #include "t591-operator.h" 5 6 #include <ceed.h> 7 #include <math.h> 8 #include <stdio.h> 9 #include <stdlib.h> 10 11 int main(int argc, char **argv) { 12 Ceed ceed; 13 CeedInt num_elem_1d = 3, num_elem = num_elem_1d * num_elem_1d, dim = 2, p_coarse = 2, p_fine = 3, q = 5; 14 CeedInt num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; 15 CeedInt num_nodes_coarse = (num_elem_1d * (p_coarse - 1) + 1) * (num_elem_1d * (p_coarse - 1) + 1); 16 CeedInt num_nodes_fine = (num_elem_1d * (p_fine - 1) + 1) * (num_elem_1d * (p_fine - 1) + 1); 17 CeedVector x_points, x_elem, q_data, u_coarse, u_fine, v_coarse, v_fine, p_mult_fine; 18 CeedElemRestriction elem_restriction_x_points, elem_restriction_q_data, elem_restriction_x, elem_restriction_u_coarse, elem_restriction_u_fine; 19 CeedBasis basis_x, basis_u_coarse, basis_u_fine; 20 CeedQFunction qf_setup, qf_mass; 21 CeedOperator op_setup, op_mass_coarse, op_mass_fine, op_prolong, op_restrict; 22 23 CeedInit(argv[1], &ceed); 24 25 // Point reference coordinates 26 CeedVectorCreate(ceed, dim * num_points, &x_points); 27 { 28 CeedScalar x_array[dim * num_points]; 29 30 for (CeedInt e = 0; e < num_elem; e++) { 31 for (CeedInt d = 0; d < dim; d++) { 32 x_array[num_points_per_elem * (e * dim + d) + 0] = 0.25; 33 x_array[num_points_per_elem * (e * dim + d) + 1] = d == 0 ? -0.25 : 0.25; 34 x_array[num_points_per_elem * (e * dim + d) + 2] = d == 0 ? 0.25 : -0.25; 35 x_array[num_points_per_elem * (e * dim + d) + 3] = 0.25; 36 } 37 } 38 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 39 } 40 { 41 CeedInt ind_x[num_elem + 1 + num_points]; 42 43 for (CeedInt i = 0; i <= num_elem; i++) ind_x[i] = num_elem + 1 + i * num_points_per_elem; 44 for (CeedInt i = 0; i < num_points; i++) ind_x[num_elem + 1 + i] = i; 45 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, 46 &elem_restriction_x_points); 47 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_q_data); 48 } 49 50 // Q data 51 CeedVectorCreate(ceed, num_points, &q_data); 52 53 // Cell coordinates 54 { 55 CeedInt p = 2, num_nodes = (num_elem_1d * (p - 1) + 1) * (num_elem_1d * (p - 1) + 1); 56 CeedInt ind_x[num_elem * p * p]; 57 58 for (CeedInt e = 0; e < num_elem; e++) { 59 CeedInt elem_xy[2] = {1, 1}, n_d[2] = {0, 0}; 60 61 for (CeedInt d = 0; d < dim; d++) n_d[d] = num_elem_1d * (p - 1) + 1; 62 { 63 CeedInt r_e = e; 64 65 for (CeedInt d = 0; d < dim; d++) { 66 elem_xy[d] = r_e % num_elem_1d; 67 r_e /= num_elem_1d; 68 } 69 } 70 CeedInt num_nodes_in_elem = p * p, *elem_nodes = ind_x + e * num_nodes_in_elem; 71 72 for (CeedInt n = 0; n < num_nodes_in_elem; n++) { 73 CeedInt g_node = 0, g_node_stride = 1, r_node = n; 74 75 for (CeedInt d = 0; d < dim; d++) { 76 g_node += (elem_xy[d] * (p - 1) + r_node % p) * g_node_stride; 77 g_node_stride *= n_d[d]; 78 r_node /= p; 79 } 80 elem_nodes[n] = p * g_node; 81 } 82 } 83 CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, 1, dim * num_nodes, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_x); 84 CeedVectorCreate(ceed, dim * num_nodes, &x_elem); 85 { 86 CeedScalar x_array[dim * num_nodes]; 87 88 for (CeedInt i = 0; i <= num_elem_1d; i++) { 89 for (CeedInt j = 0; j <= num_elem_1d; j++) { 90 x_array[(i * (num_elem_1d + 1) + j) * dim + 0] = j; 91 x_array[(i * (num_elem_1d + 1) + j) * dim + 1] = i; 92 } 93 } 94 CeedVectorSetArray(x_elem, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 95 } 96 } 97 98 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, q, CEED_GAUSS, &basis_x); 99 100 // Cell solution 101 { 102 CeedInt ind_u[num_elem * p_coarse * p_coarse]; 103 104 for (CeedInt e = 0; e < num_elem; e++) { 105 CeedInt elem_xy[2] = {1, 1}, n_d[2] = {0, 0}; 106 107 for (CeedInt d = 0; d < dim; d++) n_d[d] = num_elem_1d * (p_coarse - 1) + 1; 108 { 109 CeedInt r_e = e; 110 111 for (CeedInt d = 0; d < dim; d++) { 112 elem_xy[d] = r_e % num_elem_1d; 113 r_e /= num_elem_1d; 114 } 115 } 116 CeedInt num_nodes_in_elem = p_coarse * p_coarse, *elem_nodes = ind_u + e * num_nodes_in_elem; 117 118 for (CeedInt n = 0; n < num_nodes_in_elem; n++) { 119 CeedInt g_node = 0, g_node_stride = 1, r_node = n; 120 121 for (CeedInt d = 0; d < dim; d++) { 122 g_node += (elem_xy[d] * (p_coarse - 1) + r_node % p_coarse) * g_node_stride; 123 g_node_stride *= n_d[d]; 124 r_node /= p_coarse; 125 } 126 elem_nodes[n] = g_node; 127 } 128 } 129 CeedElemRestrictionCreate(ceed, num_elem, p_coarse * p_coarse, 1, 1, num_nodes_coarse, CEED_MEM_HOST, CEED_COPY_VALUES, ind_u, 130 &elem_restriction_u_coarse); 131 } 132 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_coarse, q, CEED_GAUSS, &basis_u_coarse); 133 { 134 CeedInt ind_u[num_elem * p_fine * p_fine]; 135 136 for (CeedInt e = 0; e < num_elem; e++) { 137 CeedInt elem_xy[2] = {1, 1}, n_d[2] = {0, 0}; 138 139 for (CeedInt d = 0; d < dim; d++) n_d[d] = num_elem_1d * (p_fine - 1) + 1; 140 { 141 CeedInt r_e = e; 142 143 for (CeedInt d = 0; d < dim; d++) { 144 elem_xy[d] = r_e % num_elem_1d; 145 r_e /= num_elem_1d; 146 } 147 } 148 CeedInt num_nodes_in_elem = p_fine * p_fine, *elem_nodes = ind_u + e * num_nodes_in_elem; 149 150 for (CeedInt n = 0; n < num_nodes_in_elem; n++) { 151 CeedInt g_node = 0, g_node_stride = 1, r_node = n; 152 153 for (CeedInt d = 0; d < dim; d++) { 154 g_node += (elem_xy[d] * (p_fine - 1) + r_node % p_fine) * g_node_stride; 155 g_node_stride *= n_d[d]; 156 r_node /= p_fine; 157 } 158 elem_nodes[n] = g_node; 159 } 160 } 161 CeedElemRestrictionCreate(ceed, num_elem, p_fine * p_fine, 1, 1, num_nodes_fine, CEED_MEM_HOST, CEED_COPY_VALUES, ind_u, 162 &elem_restriction_u_fine); 163 } 164 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_fine, q, CEED_GAUSS, &basis_u_fine); 165 166 // Setup geometric scaling 167 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 168 CeedQFunctionAddInput(qf_setup, "x", dim * dim, CEED_EVAL_GRAD); 169 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 170 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 171 172 CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 173 CeedOperatorSetField(op_setup, "x", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 174 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 175 CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 176 CeedOperatorAtPointsSetPoints(op_setup, elem_restriction_x_points, x_points); 177 178 CeedOperatorApply(op_setup, x_elem, q_data, CEED_REQUEST_IMMEDIATE); 179 180 // Mass operator 181 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 182 CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP); 183 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); 184 CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP); 185 186 CeedOperatorCreateAtPoints(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_fine); 187 CeedOperatorSetField(op_mass_fine, "u", elem_restriction_u_fine, basis_u_fine, CEED_VECTOR_ACTIVE); 188 CeedOperatorSetField(op_mass_fine, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); 189 CeedOperatorSetField(op_mass_fine, "v", elem_restriction_u_fine, basis_u_fine, CEED_VECTOR_ACTIVE); 190 CeedOperatorAtPointsSetPoints(op_mass_fine, elem_restriction_x_points, x_points); 191 192 CeedVectorCreate(ceed, num_nodes_fine, &u_fine); 193 CeedVectorCreate(ceed, num_nodes_fine, &v_fine); 194 CeedVectorCreate(ceed, num_nodes_fine, &p_mult_fine); 195 CeedVectorCreate(ceed, num_nodes_coarse, &u_coarse); 196 CeedVectorCreate(ceed, num_nodes_coarse, &v_coarse); 197 198 // Create multigrid level 199 CeedVectorSetValue(p_mult_fine, 1.0); 200 CeedOperatorMultigridLevelCreate(op_mass_fine, p_mult_fine, elem_restriction_u_coarse, basis_u_coarse, &op_mass_coarse, &op_prolong, &op_restrict); 201 202 // Coarse problem 203 CeedVectorSetValue(u_coarse, 1.0); 204 CeedOperatorApply(op_mass_coarse, u_coarse, v_coarse, CEED_REQUEST_IMMEDIATE); 205 206 // Check output 207 { 208 const CeedScalar *v_array; 209 CeedScalar sum = 0.; 210 211 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 212 for (CeedInt i = 0; i < num_nodes_coarse; i++) { 213 sum += v_array[i]; 214 } 215 CeedVectorRestoreArrayRead(v_coarse, &v_array); 216 if (fabs(sum - num_elem) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 2.0\n", sum); 217 } 218 219 // Prolong coarse u 220 CeedOperatorApply(op_prolong, u_coarse, u_fine, CEED_REQUEST_IMMEDIATE); 221 222 // Fine problem 223 CeedOperatorApply(op_mass_fine, u_fine, v_fine, CEED_REQUEST_IMMEDIATE); 224 225 // Check output 226 { 227 const CeedScalar *v_array; 228 CeedScalar sum = 0.; 229 230 CeedVectorGetArrayRead(v_fine, CEED_MEM_HOST, &v_array); 231 for (CeedInt i = 0; i < num_nodes_fine; i++) { 232 sum += v_array[i]; 233 } 234 CeedVectorRestoreArrayRead(v_fine, &v_array); 235 236 if (fabs(sum - num_elem) > 1000. * CEED_EPSILON) printf("Computed Area Fine Grid: %f != True Area: 2.0\n", sum); 237 } 238 // Restrict state to coarse grid 239 CeedOperatorApply(op_restrict, v_fine, v_coarse, CEED_REQUEST_IMMEDIATE); 240 241 // Check output 242 { 243 const CeedScalar *v_array; 244 CeedScalar sum = 0.; 245 246 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array); 247 for (CeedInt i = 0; i < num_nodes_coarse; i++) { 248 sum += v_array[i]; 249 } 250 CeedVectorRestoreArrayRead(v_coarse, &v_array); 251 if (fabs(sum - num_elem) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 2.0\n", sum); 252 } 253 254 CeedVectorDestroy(&x_points); 255 CeedVectorDestroy(&q_data); 256 CeedVectorDestroy(&x_elem); 257 CeedVectorDestroy(&u_coarse); 258 CeedVectorDestroy(&u_fine); 259 CeedVectorDestroy(&v_fine); 260 CeedVectorDestroy(&v_coarse); 261 CeedVectorDestroy(&p_mult_fine); 262 CeedElemRestrictionDestroy(&elem_restriction_x_points); 263 CeedElemRestrictionDestroy(&elem_restriction_q_data); 264 CeedElemRestrictionDestroy(&elem_restriction_x); 265 CeedElemRestrictionDestroy(&elem_restriction_u_coarse); 266 CeedElemRestrictionDestroy(&elem_restriction_u_fine); 267 CeedBasisDestroy(&basis_x); 268 CeedBasisDestroy(&basis_u_coarse); 269 CeedBasisDestroy(&basis_u_fine); 270 CeedQFunctionDestroy(&qf_setup); 271 CeedQFunctionDestroy(&qf_mass); 272 CeedOperatorDestroy(&op_setup); 273 CeedOperatorDestroy(&op_mass_coarse); 274 CeedOperatorDestroy(&op_mass_fine); 275 CeedOperatorDestroy(&op_prolong); 276 CeedOperatorDestroy(&op_restrict); 277 CeedDestroy(&ceed); 278 return 0; 279 } 280