1 /// @file 2 /// Test creation, action, and destruction for mass matrix operator with multigrid level, tensor basis and interpolation basis generation 3 /// \test Test creation, action, and destruction for mass matrix operator with multigrid level, tensor basis and interpolation basis generation 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 8 #include "t502-operator.h" 9 10 int main(int argc, char **argv) { 11 Ceed ceed; 12 CeedElemRestriction Erestrictx, Erestrictui, 13 ErestrictuCoarse, ErestrictuFine; 14 CeedBasis bx, bCoarse, bFine; 15 CeedQFunction qf_setup, qf_mass; 16 CeedOperator op_setup, op_massCoarse, op_massFine, 17 op_prolong, op_restrict; 18 CeedVector qdata, X, Ucoarse, Ufine, 19 Vcoarse, Vfine, PMultFine; 20 const CeedScalar *hv; 21 CeedInt nelem = 15, Pcoarse = 3, Pfine = 5, Q = 8, ncomp = 2; 22 CeedInt Nx = nelem+1, NuCoarse = nelem*(Pcoarse-1)+1, 23 NuFine = nelem*(Pfine-1)+1; 24 CeedInt induCoarse[nelem*Pcoarse], induFine[nelem*Pfine], 25 indx[nelem*2]; 26 CeedScalar x[Nx]; 27 CeedScalar sum; 28 29 CeedInit(argv[1], &ceed); 30 31 for (CeedInt i=0; i<Nx; i++) 32 x[i] = (CeedScalar) i / (Nx - 1); 33 for (CeedInt i=0; i<nelem; i++) { 34 indx[2*i+0] = i; 35 indx[2*i+1] = i+1; 36 } 37 // Restrictions 38 CeedElemRestrictionCreate(ceed, nelem, 2, 1, 1, Nx, CEED_MEM_HOST, 39 CEED_USE_POINTER, indx, &Erestrictx); 40 41 for (CeedInt i=0; i<nelem; i++) { 42 for (CeedInt j=0; j<Pcoarse; j++) { 43 induCoarse[Pcoarse*i+j] = i*(Pcoarse-1) + j; 44 } 45 } 46 CeedElemRestrictionCreate(ceed, nelem, Pcoarse, ncomp, NuCoarse, 47 ncomp*NuCoarse, CEED_MEM_HOST, CEED_USE_POINTER, 48 induCoarse, &ErestrictuCoarse); 49 50 for (CeedInt i=0; i<nelem; i++) { 51 for (CeedInt j=0; j<Pfine; j++) { 52 induFine[Pfine*i+j] = i*(Pfine-1) + j; 53 } 54 } 55 CeedElemRestrictionCreate(ceed, nelem, Pfine, ncomp, NuFine, 56 ncomp*NuFine, CEED_MEM_HOST, CEED_USE_POINTER, 57 induFine, &ErestrictuFine); 58 59 CeedInt stridesu[3] = {1, Q, Q}; 60 CeedElemRestrictionCreateStrided(ceed, nelem, Q, 1, Q*nelem, stridesu, 61 &Erestrictui); 62 63 // Bases 64 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bx); 65 CeedBasisCreateTensorH1Lagrange(ceed, 1, ncomp, Pcoarse, Q, CEED_GAUSS, 66 &bCoarse); 67 CeedBasisCreateTensorH1Lagrange(ceed, 1, ncomp, Pfine, Q, CEED_GAUSS, &bFine); 68 69 // QFunctions 70 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 71 CeedQFunctionAddInput(qf_setup, "weights", 1, CEED_EVAL_WEIGHT); 72 CeedQFunctionAddInput(qf_setup, "dx", 1*1, CEED_EVAL_GRAD); 73 CeedQFunctionAddOutput(qf_setup, "qdata", 1, CEED_EVAL_NONE); 74 75 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 76 CeedQFunctionAddInput(qf_mass, "qdata", 1, CEED_EVAL_NONE); 77 CeedQFunctionAddInput(qf_mass, "u", ncomp, CEED_EVAL_INTERP); 78 CeedQFunctionAddOutput(qf_mass, "v", ncomp, CEED_EVAL_INTERP); 79 80 // Operators 81 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 82 &op_setup); 83 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 84 &op_massFine); 85 86 CeedVectorCreate(ceed, Nx, &X); 87 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 88 CeedVectorCreate(ceed, nelem*Q, &qdata); 89 90 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, bx, 91 CEED_VECTOR_NONE); 92 CeedOperatorSetField(op_setup, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE); 93 CeedOperatorSetField(op_setup, "qdata", Erestrictui, CEED_BASIS_COLLOCATED, 94 CEED_VECTOR_ACTIVE); 95 96 CeedOperatorSetField(op_massFine, "qdata", Erestrictui, CEED_BASIS_COLLOCATED, 97 qdata); 98 CeedOperatorSetField(op_massFine, "u", ErestrictuFine, bFine, 99 CEED_VECTOR_ACTIVE); 100 CeedOperatorSetField(op_massFine, "v", ErestrictuFine, bFine, 101 CEED_VECTOR_ACTIVE); 102 103 CeedOperatorApply(op_setup, X, qdata, CEED_REQUEST_IMMEDIATE); 104 105 // Create multigrid level 106 CeedVectorCreate(ceed, ncomp*NuFine, &PMultFine); 107 CeedVectorSetValue(PMultFine, 1.0); 108 CeedOperatorMultigridLevelCreate(op_massFine, PMultFine, ErestrictuCoarse, 109 bCoarse, &op_massCoarse, &op_prolong, &op_restrict); 110 111 // Coarse problem 112 CeedVectorCreate(ceed, ncomp*NuCoarse, &Ucoarse); 113 CeedVectorSetValue(Ucoarse, 1.0); 114 CeedVectorCreate(ceed, ncomp*NuCoarse, &Vcoarse); 115 CeedOperatorApply(op_massCoarse, Ucoarse, Vcoarse, CEED_REQUEST_IMMEDIATE); 116 117 // Check output 118 CeedVectorGetArrayRead(Vcoarse, CEED_MEM_HOST, &hv); 119 sum = 0.; 120 for (CeedInt i=0; i<ncomp*NuCoarse; i++) { 121 sum += hv[i]; 122 } 123 if (fabs(sum-2.)>1e-10) 124 // LCOV_EXCL_START 125 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 126 // LCOV_EXCL_STOP 127 CeedVectorRestoreArrayRead(Vcoarse, &hv); 128 129 // Prolong coarse u 130 CeedVectorCreate(ceed, ncomp*NuFine, &Ufine); 131 CeedOperatorApply(op_prolong, Ucoarse, Ufine, CEED_REQUEST_IMMEDIATE); 132 133 // Fine problem 134 CeedVectorCreate(ceed, ncomp*NuFine, &Vfine); 135 CeedOperatorApply(op_massFine, Ufine, Vfine, CEED_REQUEST_IMMEDIATE); 136 137 // Check output 138 CeedVectorGetArrayRead(Vfine, CEED_MEM_HOST, &hv); 139 sum = 0.; 140 for (CeedInt i=0; i<ncomp*NuFine; i++) { 141 sum += hv[i]; 142 } 143 if (fabs(sum-2.)>1e-10) 144 // LCOV_EXCL_START 145 printf("Computed Area Fine Grid: %f != True Area: 1.0\n", sum); 146 // LCOV_EXCL_STOP 147 CeedVectorRestoreArrayRead(Vfine, &hv); 148 149 // Restrict state to coarse grid 150 CeedOperatorApply(op_restrict, Vfine, Vcoarse, CEED_REQUEST_IMMEDIATE); 151 152 // Check output 153 CeedVectorGetArrayRead(Vcoarse, CEED_MEM_HOST, &hv); 154 sum = 0.; 155 for (CeedInt i=0; i<ncomp*NuCoarse; i++) { 156 sum += hv[i]; 157 } 158 if (fabs(sum-2.)>1e-10) 159 // LCOV_EXCL_START 160 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 161 // LCOV_EXCL_STOP 162 CeedVectorRestoreArrayRead(Vcoarse, &hv); 163 164 // Cleanup 165 CeedQFunctionDestroy(&qf_setup); 166 CeedQFunctionDestroy(&qf_mass); 167 CeedOperatorDestroy(&op_setup); 168 CeedOperatorDestroy(&op_massCoarse); 169 CeedOperatorDestroy(&op_massFine); 170 CeedOperatorDestroy(&op_prolong); 171 CeedOperatorDestroy(&op_restrict); 172 CeedElemRestrictionDestroy(&ErestrictuCoarse); 173 CeedElemRestrictionDestroy(&ErestrictuFine); 174 CeedElemRestrictionDestroy(&Erestrictx); 175 CeedElemRestrictionDestroy(&Erestrictui); 176 CeedBasisDestroy(&bCoarse); 177 CeedBasisDestroy(&bFine); 178 CeedBasisDestroy(&bx); 179 CeedVectorDestroy(&X); 180 CeedVectorDestroy(&Ucoarse); 181 CeedVectorDestroy(&Ufine); 182 CeedVectorDestroy(&Vcoarse); 183 CeedVectorDestroy(&Vfine); 184 CeedVectorDestroy(&PMultFine); 185 CeedVectorDestroy(&qdata); 186 CeedDestroy(&ceed); 187 return 0; 188 } 189