1 /// @file 2 /// Test creation and use of FDM element inverse 3 /// \test Test creation and use of FDM element inverse 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <string.h> 7 #include <math.h> 8 #include "t541-operator.h" 9 10 int main(int argc, char **argv) { 11 Ceed ceed; 12 CeedElemRestriction elem_restr_x_i, elem_restr_u_i, elem_restr_qd_i; 13 CeedBasis basis_x, basis_u; 14 CeedQFunction qf_setup_diff, qf_apply; 15 CeedOperator op_setup_diff, op_apply, op_inv; 16 CeedVector q_data_diff, X, U, V, W; 17 CeedInt num_elem = 1, P = 4, Q = 5, dim = 2; 18 CeedInt num_dofs = P*P, num_qpts = num_elem*Q*Q, q_data_size = dim*(dim+1)/2; 19 CeedScalar x[dim*num_elem*(2*2)]; 20 21 CeedInit(argv[1], &ceed); 22 23 // DoF Coordinates 24 for (CeedInt i=0; i<2; i++) 25 for (CeedInt j=0; j<2; j++) { 26 x[i+j*2+0*4] = i; 27 x[i+j*2+1*4] = j; 28 } 29 CeedVectorCreate(ceed, dim*num_elem*(2*2), &X); 30 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 31 32 // Qdata Vector 33 CeedVectorCreate(ceed, q_data_size*num_qpts, &q_data_diff); 34 35 // Element Setup 36 37 // Restrictions 38 CeedInt strides_x[3] = {1, 2*2, 2*2*dim}; 39 CeedElemRestrictionCreateStrided(ceed, num_elem, 2*2, dim, dim*num_elem*2*2, 40 strides_x, &elem_restr_x_i); 41 42 CeedInt strides_u[3] = {1, P*P, P*P}; 43 CeedElemRestrictionCreateStrided(ceed, num_elem, P*P, 1, num_dofs, strides_u, 44 &elem_restr_u_i); 45 46 CeedInt strides_qd[3] = {1, Q*Q, q_data_size *Q*Q}; 47 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q, q_data_size, 48 num_qpts*q_data_size, strides_qd, &elem_restr_qd_i); 49 50 // Bases 51 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &basis_x); 52 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_u); 53 54 // QFunction - setup diff 55 CeedQFunctionCreateInterior(ceed, 1, setup_diff, setup_diff_loc, 56 &qf_setup_diff); 57 CeedQFunctionAddInput(qf_setup_diff, "dx", dim*dim, CEED_EVAL_GRAD); 58 CeedQFunctionAddInput(qf_setup_diff, "weight", 1, CEED_EVAL_WEIGHT); 59 CeedQFunctionAddOutput(qf_setup_diff, "qdata", q_data_size, CEED_EVAL_NONE); 60 61 // Operator - setup diff 62 CeedOperatorCreate(ceed, qf_setup_diff, CEED_QFUNCTION_NONE, 63 CEED_QFUNCTION_NONE, &op_setup_diff); 64 CeedOperatorSetField(op_setup_diff, "dx", elem_restr_x_i, basis_x, 65 CEED_VECTOR_ACTIVE); 66 CeedOperatorSetField(op_setup_diff, "weight", CEED_ELEMRESTRICTION_NONE, 67 basis_x, CEED_VECTOR_NONE); 68 CeedOperatorSetField(op_setup_diff, "qdata", elem_restr_qd_i, 69 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 70 71 // Apply Setup Operator 72 CeedOperatorApply(op_setup_diff, X, q_data_diff, CEED_REQUEST_IMMEDIATE); 73 74 // QFunction - apply 75 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 76 CeedQFunctionAddInput(qf_apply, "u", dim, CEED_EVAL_GRAD); 77 CeedQFunctionAddInput(qf_apply, "qdata_diff", q_data_size, CEED_EVAL_NONE); 78 CeedQFunctionAddOutput(qf_apply, "v", dim, CEED_EVAL_GRAD); 79 80 // Operator - apply 81 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 82 &op_apply); 83 CeedOperatorSetField(op_apply, "u", elem_restr_u_i, basis_u, 84 CEED_VECTOR_ACTIVE); 85 CeedOperatorSetField(op_apply, "qdata_diff", elem_restr_qd_i, 86 CEED_BASIS_COLLOCATED, q_data_diff); 87 CeedOperatorSetField(op_apply, "v", elem_restr_u_i, basis_u, 88 CEED_VECTOR_ACTIVE); 89 90 // Create FDM element inverse 91 CeedOperatorCreateFDMElementInverse(op_apply, &op_inv, CEED_REQUEST_IMMEDIATE); 92 93 // Create vectors 94 CeedVectorCreate(ceed, num_dofs, &U); 95 CeedVectorCreate(ceed, num_dofs, &V); 96 CeedVectorCreate(ceed, num_dofs, &W); 97 98 // Create Schur complement for element corners 99 CeedScalar S[16]; 100 for (CeedInt i = 0; i < 4; i++) { 101 CeedScalar *u; 102 CeedVectorSetValue(U, 0.0); 103 CeedVectorGetArray(U, CEED_MEM_HOST, &u); 104 switch (i) { 105 case 0: u[0] = 1.0; break; 106 case 1: u[P-1] = 1.0; break; 107 case 2: u[P*P-P] = 1.0; break; 108 case 3: u[P*P-1] = 1.0; break; 109 } 110 CeedVectorRestoreArray(U, &u); 111 112 CeedOperatorApply(op_inv, U, V, CEED_REQUEST_IMMEDIATE); 113 114 const CeedScalar *v; 115 CeedVectorGetArrayRead(V, CEED_MEM_HOST, &v); 116 S[0*4+i] = -v[0]; 117 S[1*4+i] = -v[P-1]; 118 S[2*4+i] = -v[P*P-P]; 119 S[3*4+i] = -v[P*P-1]; 120 CeedVectorRestoreArrayRead(V, &v); 121 } 122 CeedScalar S_inv[16]; 123 { 124 CeedScalar det; 125 S_inv[0] = S[5] * S[10] * S[15] - 126 S[5] * S[11] * S[14] - 127 S[9] * S[6] * S[15] + 128 S[9] * S[7] * S[14] + 129 S[13] * S[6] * S[11] - 130 S[13] * S[7] * S[10]; 131 132 S_inv[4] = -S[4] * S[10] * S[15] + 133 S[4] * S[11] * S[14] + 134 S[8] * S[6] * S[15] - 135 S[8] * S[7] * S[14] - 136 S[12] * S[6] * S[11] + 137 S[12] * S[7] * S[10]; 138 139 S_inv[8] = S[4] * S[9] * S[15] - 140 S[4] * S[11] * S[13] - 141 S[8] * S[5] * S[15] + 142 S[8] * S[7] * S[13] + 143 S[12] * S[5] * S[11] - 144 S[12] * S[7] * S[9]; 145 146 S_inv[12] = -S[4] * S[9] * S[14] + 147 S[4] * S[10] * S[13] + 148 S[8] * S[5] * S[14] - 149 S[8] * S[6] * S[13] - 150 S[12] * S[5] * S[10] + 151 S[12] * S[6] * S[9]; 152 153 S_inv[1] = -S[1] * S[10] * S[15] + 154 S[1] * S[11] * S[14] + 155 S[9] * S[2] * S[15] - 156 S[9] * S[3] * S[14] - 157 S[13] * S[2] * S[11] + 158 S[13] * S[3] * S[10]; 159 160 S_inv[5] = S[0] * S[10] * S[15] - 161 S[0] * S[11] * S[14] - 162 S[8] * S[2] * S[15] + 163 S[8] * S[3] * S[14] + 164 S[12] * S[2] * S[11] - 165 S[12] * S[3] * S[10]; 166 167 S_inv[9] = -S[0] * S[9] * S[15] + 168 S[0] * S[11] * S[13] + 169 S[8] * S[1] * S[15] - 170 S[8] * S[3] * S[13] - 171 S[12] * S[1] * S[11] + 172 S[12] * S[3] * S[9]; 173 174 S_inv[13] = S[0] * S[9] * S[14] - 175 S[0] * S[10] * S[13] - 176 S[8] * S[1] * S[14] + 177 S[8] * S[2] * S[13] + 178 S[12] * S[1] * S[10] - 179 S[12] * S[2] * S[9]; 180 181 S_inv[2] = S[1] * S[6] * S[15] - 182 S[1] * S[7] * S[14] - 183 S[5] * S[2] * S[15] + 184 S[5] * S[3] * S[14] + 185 S[13] * S[2] * S[7] - 186 S[13] * S[3] * S[6]; 187 188 S_inv[6] = -S[0] * S[6] * S[15] + 189 S[0] * S[7] * S[14] + 190 S[4] * S[2] * S[15] - 191 S[4] * S[3] * S[14] - 192 S[12] * S[2] * S[7] + 193 S[12] * S[3] * S[6]; 194 195 S_inv[10] = S[0] * S[5] * S[15] - 196 S[0] * S[7] * S[13] - 197 S[4] * S[1] * S[15] + 198 S[4] * S[3] * S[13] + 199 S[12] * S[1] * S[7] - 200 S[12] * S[3] * S[5]; 201 202 S_inv[14] = -S[0] * S[5] * S[14] + 203 S[0] * S[6] * S[13] + 204 S[4] * S[1] * S[14] - 205 S[4] * S[2] * S[13] - 206 S[12] * S[1] * S[6] + 207 S[12] * S[2] * S[5]; 208 209 S_inv[3] = -S[1] * S[6] * S[11] + 210 S[1] * S[7] * S[10] + 211 S[5] * S[2] * S[11] - 212 S[5] * S[3] * S[10] - 213 S[9] * S[2] * S[7] + 214 S[9] * S[3] * S[6]; 215 216 S_inv[7] = S[0] * S[6] * S[11] - 217 S[0] * S[7] * S[10] - 218 S[4] * S[2] * S[11] + 219 S[4] * S[3] * S[10] + 220 S[8] * S[2] * S[7] - 221 S[8] * S[3] * S[6]; 222 223 S_inv[11] = -S[0] * S[5] * S[11] + 224 S[0] * S[7] * S[9] + 225 S[4] * S[1] * S[11] - 226 S[4] * S[3] * S[9] - 227 S[8] * S[1] * S[7] + 228 S[8] * S[3] * S[5]; 229 230 S_inv[15] = S[0] * S[5] * S[10] - 231 S[0] * S[6] * S[9] - 232 S[4] * S[1] * S[10] + 233 S[4] * S[2] * S[9] + 234 S[8] * S[1] * S[6] - 235 S[8] * S[2] * S[5]; 236 237 det = 1/(S[0]*S_inv[0] + S[1]*S_inv[4] + S[2]*S_inv[8] + S[3]*S_inv[12]); 238 239 for (CeedInt i = 0; i < 16; i++) 240 S_inv[i] *= det; 241 } 242 243 // Set initial values 244 { 245 CeedScalar nodes[P]; 246 CeedLobattoQuadrature(P, nodes, NULL); 247 CeedScalar *u; 248 CeedVectorGetArray(U, CEED_MEM_HOST, &u); 249 for (CeedInt i=0; i<P; i++) 250 for (CeedInt j=0; j<P; j++) 251 u[i*P+j] = -(nodes[i] - 1.0)*(nodes[i] + 1.0) - 252 (nodes[j] - 1.0)*(nodes[j] + 1.0); 253 CeedVectorRestoreArray(U, &u); 254 } 255 256 // Apply original operator 257 CeedOperatorApply(op_apply, U, V, CEED_REQUEST_IMMEDIATE); 258 259 // Apply FDM element inverse 260 { 261 // -- Zero corners 262 CeedScalar *v; 263 CeedVectorGetArray(V, CEED_MEM_HOST, &v); 264 v[0] = 0.0; 265 v[P-1] = 0.0; 266 v[P*P-P] = 0.0; 267 v[P*P-1] = 0.0; 268 CeedVectorRestoreArray(V, &v); 269 270 // -- Apply FDM inverse to interior 271 CeedOperatorApply(op_inv, V, W, CEED_REQUEST_IMMEDIATE); 272 273 // -- Pick off corners 274 const CeedScalar *w; 275 CeedScalar w_Pi[4]; 276 CeedVectorGetArrayRead(W, CEED_MEM_HOST, &w); 277 w_Pi[0] = w[0]; 278 w_Pi[1] = w[P-1]; 279 w_Pi[2] = w[P*P-P]; 280 w_Pi[3] = w[P*P-1]; 281 CeedVectorRestoreArrayRead(W, &w); 282 283 // -- Apply inverse of Schur complement 284 CeedScalar v_Pi[4]; 285 for (CeedInt i=0; i<4; i++) { 286 CeedScalar sum = 0.0; 287 for (CeedInt j=0; j<4; j++) { 288 sum += w_Pi[j] * S_inv[i*4+j]; 289 } 290 v_Pi[i] = sum; 291 } 292 293 // -- Set corners 294 CeedVectorGetArray(V, CEED_MEM_HOST, &v); 295 v[0] = v_Pi[0]; 296 v[P-1] = v_Pi[1]; 297 v[P*P-P] = v_Pi[2]; 298 v[P*P-1] = v_Pi[3]; 299 CeedVectorRestoreArray(V, &v); 300 301 // -- Apply full FDM inverse again 302 CeedOperatorApply(op_inv, V, W, CEED_REQUEST_IMMEDIATE); 303 } 304 305 306 // Check output 307 { 308 const CeedScalar *u, *w; 309 CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); 310 CeedVectorGetArrayRead(W, CEED_MEM_HOST, &w); 311 for (CeedInt i=0; i<P; i++) 312 for (CeedInt j=0; j<P; j++) 313 if (fabs(u[i*P+j] - w[i*P+j]) > 2e-3) 314 // LCOV_EXCL_START 315 printf("[%d, %d] Error in inverse: %e != %e\n", i, j, w[i*P+j], u[i*P+j]); 316 // LCOV_EXCL_STOP 317 CeedVectorRestoreArrayRead(U, &u); 318 CeedVectorRestoreArrayRead(W, &w); 319 } 320 321 // Cleanup 322 CeedQFunctionDestroy(&qf_setup_diff); 323 CeedQFunctionDestroy(&qf_apply); 324 CeedOperatorDestroy(&op_setup_diff); 325 CeedOperatorDestroy(&op_apply); 326 CeedOperatorDestroy(&op_inv); 327 CeedElemRestrictionDestroy(&elem_restr_u_i); 328 CeedElemRestrictionDestroy(&elem_restr_x_i); 329 CeedElemRestrictionDestroy(&elem_restr_qd_i); 330 CeedBasisDestroy(&basis_x); 331 CeedBasisDestroy(&basis_u); 332 CeedVectorDestroy(&X); 333 CeedVectorDestroy(&q_data_diff); 334 CeedVectorDestroy(&U); 335 CeedVectorDestroy(&V); 336 CeedVectorDestroy(&W); 337 CeedDestroy(&ceed); 338 return 0; 339 } 340