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