1 /// @file 2 /// Test creation, action, and destruction for mass matrix operator with multigrid level, tensor basis 3 /// \test Test creation, action, and destruction for mass matrix operator with multigrid level, tensor basis 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, bu; 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, Pfine, Q, CEED_GAUSS, &bu); 66 67 // QFunctions 68 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 69 CeedQFunctionAddInput(qf_setup, "weights", 1, CEED_EVAL_WEIGHT); 70 CeedQFunctionAddInput(qf_setup, "dx", 1*1, CEED_EVAL_GRAD); 71 CeedQFunctionAddOutput(qf_setup, "qdata", 1, CEED_EVAL_NONE); 72 73 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 74 CeedQFunctionAddInput(qf_mass, "qdata", 1, CEED_EVAL_NONE); 75 CeedQFunctionAddInput(qf_mass, "u", ncomp, CEED_EVAL_INTERP); 76 CeedQFunctionAddOutput(qf_mass, "v", ncomp, CEED_EVAL_INTERP); 77 78 // Operators 79 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 80 &op_setup); 81 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 82 &op_massFine); 83 84 CeedVectorCreate(ceed, Nx, &X); 85 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 86 CeedVectorCreate(ceed, nelem*Q, &qdata); 87 88 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, bx, 89 CEED_VECTOR_NONE); 90 CeedOperatorSetField(op_setup, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE); 91 CeedOperatorSetField(op_setup, "qdata", Erestrictui, CEED_BASIS_COLLOCATED, 92 CEED_VECTOR_ACTIVE); 93 94 CeedOperatorSetField(op_massFine, "qdata", Erestrictui, CEED_BASIS_COLLOCATED, 95 qdata); 96 CeedOperatorSetField(op_massFine, "u", ErestrictuFine, bu, CEED_VECTOR_ACTIVE); 97 CeedOperatorSetField(op_massFine, "v", ErestrictuFine, bu, CEED_VECTOR_ACTIVE); 98 99 CeedOperatorApply(op_setup, X, qdata, CEED_REQUEST_IMMEDIATE); 100 101 // Create multigrid level 102 CeedVectorCreate(ceed, ncomp*NuFine, &PMultFine); 103 CeedVectorSetValue(PMultFine, 1.0); 104 CeedBasis buCoarse, bCtoF; 105 CeedBasisCreateTensorH1Lagrange(ceed, 1, ncomp, Pcoarse, Q, CEED_GAUSS, 106 &buCoarse); 107 CeedBasisCreateTensorH1Lagrange(ceed, 1, ncomp, Pcoarse, Pfine, 108 CEED_GAUSS_LOBATTO, &bCtoF); 109 const CeedScalar *interpCtoF; 110 CeedBasisGetInterp1D(bCtoF, &interpCtoF); 111 CeedOperatorMultigridLevelCreateTensorH1(op_massFine, PMultFine, 112 ErestrictuCoarse, buCoarse, interpCtoF, &op_massCoarse, &op_prolong, 113 &op_restrict); 114 115 // Coarse problem 116 CeedVectorCreate(ceed, ncomp*NuCoarse, &Ucoarse); 117 CeedVectorSetValue(Ucoarse, 1.0); 118 CeedVectorCreate(ceed, ncomp*NuCoarse, &Vcoarse); 119 CeedOperatorApply(op_massCoarse, Ucoarse, Vcoarse, CEED_REQUEST_IMMEDIATE); 120 121 // Check output 122 CeedVectorGetArrayRead(Vcoarse, CEED_MEM_HOST, &hv); 123 sum = 0.; 124 for (CeedInt i=0; i<ncomp*NuCoarse; i++) { 125 sum += hv[i]; 126 } 127 if (fabs(sum-2.)>1e-10) 128 // LCOV_EXCL_START 129 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 130 // LCOV_EXCL_STOP 131 CeedVectorRestoreArrayRead(Vcoarse, &hv); 132 133 // Prolong coarse u 134 CeedVectorCreate(ceed, ncomp*NuFine, &Ufine); 135 CeedOperatorApply(op_prolong, Ucoarse, Ufine, CEED_REQUEST_IMMEDIATE); 136 137 // Fine problem 138 CeedVectorCreate(ceed, ncomp*NuFine, &Vfine); 139 CeedOperatorApply(op_massFine, Ufine, Vfine, CEED_REQUEST_IMMEDIATE); 140 141 // Check output 142 CeedVectorGetArrayRead(Vfine, CEED_MEM_HOST, &hv); 143 sum = 0.; 144 for (CeedInt i=0; i<ncomp*NuFine; i++) { 145 sum += hv[i]; 146 } 147 if (fabs(sum-2.)>1e-10) 148 // LCOV_EXCL_START 149 printf("Computed Area Fine Grid: %f != True Area: 1.0\n", sum); 150 // LCOV_EXCL_STOP 151 CeedVectorRestoreArrayRead(Vfine, &hv); 152 153 // Restrict state to coarse grid 154 CeedOperatorApply(op_restrict, Vfine, Vcoarse, CEED_REQUEST_IMMEDIATE); 155 156 // Check output 157 CeedVectorGetArrayRead(Vcoarse, CEED_MEM_HOST, &hv); 158 sum = 0.; 159 for (CeedInt i=0; i<ncomp*NuCoarse; i++) { 160 sum += hv[i]; 161 } 162 if (fabs(sum-2.)>1e-10) 163 // LCOV_EXCL_START 164 printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum); 165 // LCOV_EXCL_STOP 166 CeedVectorRestoreArrayRead(Vcoarse, &hv); 167 168 // Cleanup 169 CeedQFunctionDestroy(&qf_setup); 170 CeedQFunctionDestroy(&qf_mass); 171 CeedOperatorDestroy(&op_setup); 172 CeedOperatorDestroy(&op_massCoarse); 173 CeedOperatorDestroy(&op_massFine); 174 CeedOperatorDestroy(&op_prolong); 175 CeedOperatorDestroy(&op_restrict); 176 CeedElemRestrictionDestroy(&ErestrictuCoarse); 177 CeedElemRestrictionDestroy(&ErestrictuFine); 178 CeedElemRestrictionDestroy(&Erestrictx); 179 CeedElemRestrictionDestroy(&Erestrictui); 180 CeedBasisDestroy(&bu); 181 CeedBasisDestroy(&bx); 182 CeedVectorDestroy(&X); 183 CeedVectorDestroy(&Ucoarse); 184 CeedVectorDestroy(&Ufine); 185 CeedVectorDestroy(&Vcoarse); 186 CeedVectorDestroy(&Vfine); 187 CeedVectorDestroy(&PMultFine); 188 CeedVectorDestroy(&qdata); 189 CeedDestroy(&ceed); 190 return 0; 191 } 192