1 /// @file 2 /// Test creation reuse of the same QFunction for multiple operators 3 /// \test Test creation reuse of the same QFunction for multiple operators 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t502-operator.h" 8 9 int main(int argc, char **argv) { 10 Ceed ceed; 11 CeedInterlaceMode imode = CEED_NONINTERLACED, 12 imodeu = CEED_INTERLACED;; 13 CeedElemRestriction Erestrictx, Erestrictu, 14 Erestrictui_small, Erestrictui_large; 15 CeedBasis bx_small, bx_large, bu_small, bu_large; 16 CeedQFunction qf_setup, qf_mass; 17 CeedOperator op_setup_small, op_mass_small, 18 op_setup_large, op_mass_large; 19 CeedVector qdata_small, qdata_large, X, U, V; 20 CeedScalar *hu; 21 const CeedScalar *hv; 22 CeedInt nelem = 15, P = 5, Q = 8, scale = 3; 23 CeedInt Nx = nelem+1, Nu = nelem*(P-1)+1; 24 CeedInt indx[nelem*2], indu[nelem*P]; 25 CeedScalar x[Nx]; 26 CeedScalar sum1, sum2; 27 28 CeedInit(argv[1], &ceed); 29 for (CeedInt i=0; i<Nx; i++) x[i] = (CeedScalar) i / (Nx - 1); 30 for (CeedInt i=0; i<nelem; i++) { 31 indx[2*i+0] = i; 32 indx[2*i+1] = i+1; 33 } 34 // Restrictions 35 CeedElemRestrictionCreate(ceed, imode, nelem, 2, Nx, 1, CEED_MEM_HOST, 36 CEED_USE_POINTER, indx, &Erestrictx); 37 38 for (CeedInt i=0; i<nelem; i++) { 39 for (CeedInt j=0; j<P; j++) { 40 indu[P*i+j] = i*(P-1) + j; 41 } 42 } 43 CeedElemRestrictionCreate(ceed, imodeu, nelem, P, Nu, 2, CEED_MEM_HOST, 44 CEED_USE_POINTER, indu, &Erestrictu); 45 CeedInt stridesu_small[3] = {1, Q, Q}; 46 CeedElemRestrictionCreateStrided(ceed, nelem, Q, Q*nelem, 1, stridesu_small, 47 &Erestrictui_small); 48 CeedInt stridesu_large[3] = {1, Q*scale, Q*scale}; 49 CeedElemRestrictionCreateStrided(ceed, nelem, Q*scale, Q*nelem*scale, 1, 50 stridesu_large, &Erestrictui_large); 51 52 // Bases 53 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bx_small); 54 CeedBasisCreateTensorH1Lagrange(ceed, 1, 2, P, Q, CEED_GAUSS, &bu_small); 55 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q*scale, CEED_GAUSS, &bx_large); 56 CeedBasisCreateTensorH1Lagrange(ceed, 1, 2, P, Q*scale, CEED_GAUSS, &bu_large); 57 58 // QFunctions 59 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 60 CeedQFunctionAddInput(qf_setup, "_weight", 1, CEED_EVAL_WEIGHT); 61 CeedQFunctionAddInput(qf_setup, "x", 1, CEED_EVAL_GRAD); 62 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 63 64 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 65 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); 66 CeedQFunctionAddInput(qf_mass, "u", 2, CEED_EVAL_INTERP); 67 CeedQFunctionAddOutput(qf_mass, "v", 2, CEED_EVAL_INTERP); 68 69 // Input vector 70 CeedVectorCreate(ceed, Nx, &X); 71 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 72 73 // 'Small' Operators 74 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 75 &op_setup_small); 76 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 77 &op_mass_small); 78 79 CeedVectorCreate(ceed, nelem*Q, &qdata_small); 80 81 CeedOperatorSetField(op_setup_small, "_weight", CEED_ELEMRESTRICTION_NONE, 82 bx_small, CEED_VECTOR_NONE); 83 CeedOperatorSetField(op_setup_small, "x", Erestrictx, 84 bx_small, CEED_VECTOR_ACTIVE); 85 CeedOperatorSetField(op_setup_small, "rho", Erestrictui_small, 86 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 87 88 CeedOperatorSetField(op_mass_small, "rho", Erestrictui_small, 89 CEED_BASIS_COLLOCATED, qdata_small); 90 CeedOperatorSetField(op_mass_small, "u", Erestrictu, 91 bu_small, CEED_VECTOR_ACTIVE); 92 CeedOperatorSetField(op_mass_small, "v", Erestrictu, 93 bu_small, CEED_VECTOR_ACTIVE); 94 95 // 'Large' operators 96 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 97 &op_setup_large); 98 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 99 &op_mass_large); 100 101 CeedVectorCreate(ceed, nelem*Q*scale, &qdata_large); 102 103 CeedOperatorSetField(op_setup_large, "_weight", CEED_ELEMRESTRICTION_NONE, 104 bx_large, CEED_VECTOR_NONE); 105 CeedOperatorSetField(op_setup_large, "x", Erestrictx, 106 bx_large, CEED_VECTOR_ACTIVE); 107 CeedOperatorSetField(op_setup_large, "rho", Erestrictui_large, 108 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 109 110 CeedOperatorSetField(op_mass_large, "rho", Erestrictui_large, 111 CEED_BASIS_COLLOCATED, qdata_large); 112 CeedOperatorSetField(op_mass_large, "u", Erestrictu, 113 bu_large, CEED_VECTOR_ACTIVE); 114 CeedOperatorSetField(op_mass_large, "v", Erestrictu, 115 bu_large, CEED_VECTOR_ACTIVE); 116 117 // Setup 118 CeedOperatorApply(op_setup_small, X, qdata_small, CEED_REQUEST_IMMEDIATE); 119 CeedOperatorApply(op_setup_large, X, qdata_large, CEED_REQUEST_IMMEDIATE); 120 121 CeedVectorCreate(ceed, 2*Nu, &U); 122 CeedVectorGetArray(U, CEED_MEM_HOST, &hu); 123 for (int i = 0; i < Nu; i++) { 124 hu[2*i] = 1.0; 125 hu[2*i+1] = 2.0; 126 } 127 CeedVectorRestoreArray(U, &hu); 128 CeedVectorCreate(ceed, 2*Nu, &V); 129 130 // 'Small' operator 131 CeedOperatorApply(op_mass_small, U, V, CEED_REQUEST_IMMEDIATE); 132 133 // Check output 134 CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv); 135 sum1 = 0.; sum2 = 0.; 136 for (CeedInt i=0; i<Nu; i++) { 137 sum1 += hv[2*i]; 138 sum2 += hv[2*i+1]; 139 } 140 if (fabs(sum1-1.)>1e-10) printf("Computed Area: %f != True Area: 1.0\n", sum1); 141 if (fabs(sum2-2.)>1e-10) printf("Computed Area: %f != True Area: 2.0\n", sum2); 142 CeedVectorRestoreArrayRead(V, &hv); 143 144 // 'Large' operator 145 CeedOperatorApply(op_mass_large, U, V, CEED_REQUEST_IMMEDIATE); 146 147 // Check output 148 CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv); 149 sum1 = 0.; sum2 = 0.; 150 for (CeedInt i=0; i<Nu; i++) { 151 sum1 += hv[2*i]; 152 sum2 += hv[2*i+1]; 153 } 154 if (fabs(sum1-1.)>1e-10) printf("Computed Area: %f != True Area: 1.0\n", sum1); 155 if (fabs(sum2-2.)>1e-10) printf("Computed Area: %f != True Area: 2.0\n", sum2); 156 CeedVectorRestoreArrayRead(V, &hv); 157 158 CeedQFunctionDestroy(&qf_setup); 159 CeedQFunctionDestroy(&qf_mass); 160 CeedOperatorDestroy(&op_setup_small); 161 CeedOperatorDestroy(&op_mass_small); 162 CeedOperatorDestroy(&op_setup_large); 163 CeedOperatorDestroy(&op_mass_large); 164 CeedElemRestrictionDestroy(&Erestrictu); 165 CeedElemRestrictionDestroy(&Erestrictx); 166 CeedElemRestrictionDestroy(&Erestrictui_small); 167 CeedElemRestrictionDestroy(&Erestrictui_large); 168 CeedBasisDestroy(&bu_small); 169 CeedBasisDestroy(&bx_small); 170 CeedBasisDestroy(&bu_large); 171 CeedBasisDestroy(&bx_large); 172 CeedVectorDestroy(&X); 173 CeedVectorDestroy(&U); 174 CeedVectorDestroy(&V); 175 CeedVectorDestroy(&qdata_small); 176 CeedVectorDestroy(&qdata_large); 177 CeedDestroy(&ceed); 178 return 0; 179 } 180