1 /// @file 2 /// Test assembly of mass and Poisson operator diagonal 3 /// \test Test assembly of mass and Poisson operator diagonal 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t320-basis.h" 8 #include "t535-operator.h" 9 10 int main(int argc, char **argv) { 11 Ceed ceed; 12 CeedInterlaceMode imode = CEED_NONINTERLACED; 13 CeedElemRestriction Erestrictx, Erestrictu, 14 Erestrictxi, Erestrictui, 15 Erestrictqi; 16 CeedBasis bx, bu; 17 CeedQFunction qf_setup_mass, qf_setup_diff, qf_apply; 18 CeedOperator op_setup_mass, op_setup_diff, op_apply; 19 CeedVector qdata_mass, qdata_diff, X, A, U, V; 20 CeedInt nelem = 12, dim = 2, P = 6, Q = 4; 21 CeedInt nx = 3, ny = 2; 22 CeedInt row, col, offset; 23 CeedInt ndofs = (nx*2+1)*(ny*2+1), nqpts = nelem*Q*Q; 24 CeedInt indx[nelem*P*P]; 25 CeedScalar x[dim*ndofs], assembledTrue[ndofs]; 26 CeedScalar qref[dim*Q], qweight[Q]; 27 CeedScalar interp[P*Q], grad[dim*P*Q]; 28 CeedScalar *u; 29 const CeedScalar *a, *v; 30 31 CeedInit(argv[1], &ceed); 32 33 // DoF Coordinates 34 for (CeedInt i=0; i<ndofs; i++) { 35 x[i] = (1. / (nx*2)) * (CeedScalar) (i % (nx*2+1)); 36 x[i+ndofs] = (1. / (ny*2)) * (CeedScalar) (i / (nx*2+1)); 37 } 38 CeedVectorCreate(ceed, dim*ndofs, &X); 39 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 40 41 // Qdata Vectors 42 CeedVectorCreate(ceed, nqpts, &qdata_mass); 43 CeedVectorCreate(ceed, nqpts*dim*(dim+1)/2, &qdata_diff); 44 45 // Element Setup 46 for (CeedInt i=0; i<nelem/2; i++) { 47 col = i % nx; 48 row = i / nx; 49 offset = col*2 + row*(nx*2+1)*2; 50 51 indx[i*2*P + 0] = 2 + offset; 52 indx[i*2*P + 1] = 9 + offset; 53 indx[i*2*P + 2] = 16 + offset; 54 indx[i*2*P + 3] = 1 + offset; 55 indx[i*2*P + 4] = 8 + offset; 56 indx[i*2*P + 5] = 0 + offset; 57 58 indx[i*2*P + 6] = 14 + offset; 59 indx[i*2*P + 7] = 7 + offset; 60 indx[i*2*P + 8] = 0 + offset; 61 indx[i*2*P + 9] = 15 + offset; 62 indx[i*2*P + 10] = 8 + offset; 63 indx[i*2*P + 11] = 16 + offset; 64 } 65 66 // Restrictions 67 CeedElemRestrictionCreate(ceed, imode, nelem, P, ndofs, dim, CEED_MEM_HOST, 68 CEED_USE_POINTER, indx, &Erestrictx); 69 CeedInt stridesx[3] = {1, P, P*dim}; 70 CeedElemRestrictionCreateStrided(ceed, nelem, P, nelem*P, dim, stridesx, 71 &Erestrictxi); 72 73 CeedElemRestrictionCreate(ceed, imode, nelem, P, ndofs, 1, CEED_MEM_HOST, 74 CEED_USE_POINTER, indx, &Erestrictu); 75 CeedInt stridesu[3] = {1, Q, Q}; 76 CeedElemRestrictionCreateStrided(ceed, nelem, Q, nqpts, 1, stridesu, 77 &Erestrictui); 78 79 CeedInt stridesqd[3] = {1, Q, Q *dim *(dim+1)/2}; 80 CeedElemRestrictionCreateStrided(ceed, nelem, Q, nqpts, dim*(dim+1)/2, 81 stridesqd, &Erestrictqi); 82 83 // Bases 84 buildmats(qref, qweight, interp, grad); 85 CeedBasisCreateH1(ceed, CEED_TRIANGLE, dim, P, Q, interp, grad, qref, 86 qweight, &bx); 87 88 buildmats(qref, qweight, interp, grad); 89 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, P, Q, interp, grad, qref, 90 qweight, &bu); 91 92 // QFunction - setup mass 93 CeedQFunctionCreateInterior(ceed, 1, setup_mass, setup_mass_loc, 94 &qf_setup_mass); 95 CeedQFunctionAddInput(qf_setup_mass, "dx", dim*dim, CEED_EVAL_GRAD); 96 CeedQFunctionAddInput(qf_setup_mass, "_weight", 1, CEED_EVAL_WEIGHT); 97 CeedQFunctionAddOutput(qf_setup_mass, "qdata", 1, CEED_EVAL_NONE); 98 99 // Operator - setup mass 100 CeedOperatorCreate(ceed, qf_setup_mass, CEED_QFUNCTION_NONE, 101 CEED_QFUNCTION_NONE, &op_setup_mass); 102 CeedOperatorSetField(op_setup_mass, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE); 103 CeedOperatorSetField(op_setup_mass, "_weight", Erestrictxi, bx, 104 CEED_VECTOR_NONE); 105 CeedOperatorSetField(op_setup_mass, "qdata", Erestrictui, 106 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 107 108 // QFunction - setup diff 109 CeedQFunctionCreateInterior(ceed, 1, setup_diff, setup_diff_loc, 110 &qf_setup_diff); 111 CeedQFunctionAddInput(qf_setup_diff, "dx", dim*dim, CEED_EVAL_GRAD); 112 CeedQFunctionAddInput(qf_setup_diff, "_weight", 1, CEED_EVAL_WEIGHT); 113 CeedQFunctionAddOutput(qf_setup_diff, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 114 115 // Operator - setup diff 116 CeedOperatorCreate(ceed, qf_setup_diff, CEED_QFUNCTION_NONE, 117 CEED_QFUNCTION_NONE, &op_setup_diff); 118 CeedOperatorSetField(op_setup_diff, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE); 119 CeedOperatorSetField(op_setup_diff, "_weight", Erestrictxi, bx, 120 CEED_VECTOR_NONE); 121 CeedOperatorSetField(op_setup_diff, "qdata", Erestrictqi, 122 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 123 124 // Apply Setup Operators 125 CeedOperatorApply(op_setup_mass, X, qdata_mass, CEED_REQUEST_IMMEDIATE); 126 CeedOperatorApply(op_setup_diff, X, qdata_diff, CEED_REQUEST_IMMEDIATE); 127 128 // QFunction - apply 129 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 130 CeedQFunctionAddInput(qf_apply, "du", dim, CEED_EVAL_GRAD); 131 CeedQFunctionAddInput(qf_apply, "qdata_mass", 1, CEED_EVAL_NONE); 132 CeedQFunctionAddInput(qf_apply, "qdata_diff", dim*(dim+1)/2, CEED_EVAL_NONE); 133 CeedQFunctionAddInput(qf_apply, "u", 1, CEED_EVAL_INTERP); 134 CeedQFunctionAddOutput(qf_apply, "v", 1, CEED_EVAL_INTERP); 135 CeedQFunctionAddOutput(qf_apply, "dv", dim, CEED_EVAL_GRAD); 136 137 // Operator - apply 138 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 139 &op_apply); 140 CeedOperatorSetField(op_apply, "du", Erestrictu, bu, CEED_VECTOR_ACTIVE); 141 CeedOperatorSetField(op_apply, "qdata_mass", Erestrictui, 142 CEED_BASIS_COLLOCATED, qdata_mass); 143 CeedOperatorSetField(op_apply, "qdata_diff", Erestrictqi, 144 CEED_BASIS_COLLOCATED, qdata_diff); 145 CeedOperatorSetField(op_apply, "u", Erestrictu, bu, CEED_VECTOR_ACTIVE); 146 CeedOperatorSetField(op_apply, "v", Erestrictu, bu, CEED_VECTOR_ACTIVE); 147 CeedOperatorSetField(op_apply, "dv", Erestrictu, bu, CEED_VECTOR_ACTIVE); 148 149 // Assemble diagonal 150 CeedOperatorAssembleLinearDiagonal(op_apply, &A, CEED_REQUEST_IMMEDIATE); 151 152 // Manually assemble diagonal 153 CeedVectorCreate(ceed, ndofs, &U); 154 CeedVectorSetValue(U, 0.0); 155 CeedVectorCreate(ceed, ndofs, &V); 156 for (int i=0; i<ndofs; i++) { 157 // Set input 158 CeedVectorGetArray(U, CEED_MEM_HOST, &u); 159 u[i] = 1.0; 160 if (i) 161 u[i-1] = 0.0; 162 CeedVectorRestoreArray(U, &u); 163 164 // Compute diag entry for DoF i 165 CeedOperatorApply(op_apply, U, V, CEED_REQUEST_IMMEDIATE); 166 167 // Retrieve entry 168 CeedVectorGetArrayRead(V, CEED_MEM_HOST, &v); 169 assembledTrue[i] = v[i]; 170 CeedVectorRestoreArrayRead(V, &v); 171 } 172 173 // Check output 174 CeedVectorGetArrayRead(A, CEED_MEM_HOST, &a); 175 for (int i=0; i<ndofs; i++) 176 if (fabs(a[i] - assembledTrue[i]) > 1e-14) 177 // LCOV_EXCL_START 178 printf("[%d] Error in assembly: %f != %f\n", i, a[i], assembledTrue[i]); 179 // LCOV_EXCL_STOP 180 181 // Cleanup 182 CeedQFunctionDestroy(&qf_setup_mass); 183 CeedQFunctionDestroy(&qf_setup_diff); 184 CeedQFunctionDestroy(&qf_apply); 185 CeedOperatorDestroy(&op_setup_mass); 186 CeedOperatorDestroy(&op_setup_diff); 187 CeedOperatorDestroy(&op_apply); 188 CeedElemRestrictionDestroy(&Erestrictu); 189 CeedElemRestrictionDestroy(&Erestrictx); 190 CeedElemRestrictionDestroy(&Erestrictui); 191 CeedElemRestrictionDestroy(&Erestrictxi); 192 CeedElemRestrictionDestroy(&Erestrictqi); 193 CeedBasisDestroy(&bu); 194 CeedBasisDestroy(&bx); 195 CeedVectorDestroy(&X); 196 CeedVectorDestroy(&A); 197 CeedVectorDestroy(&qdata_mass); 198 CeedVectorDestroy(&qdata_diff); 199 CeedVectorDestroy(&U); 200 CeedVectorDestroy(&V); 201 CeedDestroy(&ceed); 202 return 0; 203 } 204