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