1 #include "../include/libceedsetup.h" 2 #include "../include/petscutils.h" 3 #include <stdio.h> 4 5 // ----------------------------------------------------------------------------- 6 // Destroy libCEED operator objects 7 // ----------------------------------------------------------------------------- 8 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) { 9 int ierr; 10 11 CeedVectorDestroy(&data->qdata); 12 CeedVectorDestroy(&data->Xceed); 13 CeedVectorDestroy(&data->Yceed); 14 CeedBasisDestroy(&data->basisx); 15 CeedBasisDestroy(&data->basisu); 16 CeedElemRestrictionDestroy(&data->Erestrictu); 17 CeedElemRestrictionDestroy(&data->Erestrictx); 18 CeedElemRestrictionDestroy(&data->Erestrictui); 19 CeedElemRestrictionDestroy(&data->Erestrictqdi); 20 CeedQFunctionDestroy(&data->qfApply); 21 CeedOperatorDestroy(&data->opApply); 22 if (i > 0) { 23 CeedOperatorDestroy(&data->opProlong); 24 CeedBasisDestroy(&data->basisctof); 25 CeedOperatorDestroy(&data->opRestrict); 26 } 27 ierr = PetscFree(data); CHKERRQ(ierr); 28 29 PetscFunctionReturn(0); 30 }; 31 32 // ----------------------------------------------------------------------------- 33 // Set up libCEED for a given degree 34 // ----------------------------------------------------------------------------- 35 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, 36 CeedInt topodim, CeedInt qextra, 37 PetscInt ncompx, PetscInt ncompu, 38 PetscInt gsize, PetscInt xlsize, 39 bpData bpData, CeedData data, 40 PetscBool setup_rhs, CeedVector rhsceed, 41 CeedVector *target) { 42 int ierr; 43 DM dmcoord; 44 Vec coords; 45 const PetscScalar *coordArray; 46 CeedBasis basisx, basisu; 47 CeedElemRestriction Erestrictx, Erestrictu, Erestrictui, Erestrictqdi; 48 CeedQFunction qfSetupGeo, qfApply; 49 CeedOperator opSetupGeo, opApply; 50 CeedVector xcoord, qdata, Xceed, Yceed; 51 CeedInt P, Q, nqpts, cStart, cEnd, nelem, qdatasize = bpData.qdatasize; 52 CeedScalar R = 1, // radius of the sphere 53 l = 1.0/PetscSqrtReal(3.0); // half edge of the inscribed cube 54 55 // CEED bases 56 P = degree + 1; 57 Q = P + qextra; 58 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompu, P, Q, bpData.qmode, 59 &basisu); 60 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompx, 2, Q, bpData.qmode, 61 &basisx); 62 CeedBasisGetNumQuadraturePoints(basisu, &nqpts); 63 64 // CEED restrictions 65 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 66 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 67 CHKERRQ(ierr); 68 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, topodim, 0, 0, 0, 69 &Erestrictx); 70 CHKERRQ(ierr); 71 ierr = CreateRestrictionFromPlex(ceed, dm, P, topodim, 0, 0, 0, &Erestrictu); 72 CHKERRQ(ierr); 73 74 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); 75 nelem = cEnd - cStart; 76 77 CeedElemRestrictionCreateStrided(ceed, nelem, nqpts, ncompu, ncompu*nelem*nqpts, 78 CEED_STRIDES_BACKEND, &Erestrictui); 79 CeedElemRestrictionCreateStrided(ceed, nelem, nqpts, qdatasize, 80 qdatasize*nelem*nqpts, 81 CEED_STRIDES_BACKEND, &Erestrictqdi); 82 83 // Element coordinates 84 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 85 ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr); 86 87 CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL); 88 CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, 89 (PetscScalar *)coordArray); 90 ierr = VecRestoreArrayRead(coords, &coordArray); 91 92 // Create the persistent vectors that will be needed in setup and apply 93 CeedVectorCreate(ceed, qdatasize*nelem*nqpts, &qdata); 94 CeedVectorCreate(ceed, xlsize, &Xceed); 95 CeedVectorCreate(ceed, xlsize, &Yceed); 96 97 // Create the QFunction that builds the context data 98 CeedQFunctionCreateInterior(ceed, 1, bpData.setupgeo, bpData.setupgeofname, 99 &qfSetupGeo); 100 CeedQFunctionAddInput(qfSetupGeo, "x", ncompx, CEED_EVAL_INTERP); 101 CeedQFunctionAddInput(qfSetupGeo, "dx", ncompx*topodim, CEED_EVAL_GRAD); 102 CeedQFunctionAddInput(qfSetupGeo, "weight", 1, CEED_EVAL_WEIGHT); 103 CeedQFunctionAddOutput(qfSetupGeo, "qdata", qdatasize, CEED_EVAL_NONE); 104 105 // Create the operator that builds the quadrature data 106 CeedOperatorCreate(ceed, qfSetupGeo, NULL, NULL, &opSetupGeo); 107 CeedOperatorSetField(opSetupGeo, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE); 108 CeedOperatorSetField(opSetupGeo, "dx", Erestrictx, basisx, CEED_VECTOR_ACTIVE); 109 CeedOperatorSetField(opSetupGeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx, 110 CEED_VECTOR_NONE); 111 CeedOperatorSetField(opSetupGeo, "qdata", Erestrictqdi, 112 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 113 114 // Setup qdata 115 CeedOperatorApply(opSetupGeo, xcoord, qdata, CEED_REQUEST_IMMEDIATE); 116 117 // Set up PDE operator 118 CeedInt inscale = bpData.inmode == CEED_EVAL_GRAD ? topodim : 1; 119 CeedInt outscale = bpData.outmode == CEED_EVAL_GRAD ? topodim : 1; 120 CeedQFunctionCreateInterior(ceed, 1, bpData.apply, bpData.applyfname, &qfApply); 121 CeedQFunctionAddInput(qfApply, "u", ncompu*inscale, bpData.inmode); 122 CeedQFunctionAddInput(qfApply, "qdata", qdatasize, CEED_EVAL_NONE); 123 CeedQFunctionAddOutput(qfApply, "v", ncompu*outscale, bpData.outmode); 124 125 // Create the mass or diff operator 126 CeedOperatorCreate(ceed, qfApply, NULL, NULL, &opApply); 127 CeedOperatorSetField(opApply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE); 128 CeedOperatorSetField(opApply, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED, 129 qdata); 130 CeedOperatorSetField(opApply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE); 131 132 // Set up RHS if needed 133 if (setup_rhs) { 134 CeedQFunction qfSetupRHS; 135 CeedOperator opSetupRHS; 136 CeedVectorCreate(ceed, nelem*nqpts*ncompu, target); 137 138 // Create the q-function that sets up the RHS and true solution 139 CeedQFunctionCreateInterior(ceed, 1, bpData.setuprhs, bpData.setuprhsfname, 140 &qfSetupRHS); 141 CeedQFunctionAddInput(qfSetupRHS, "x", ncompx, CEED_EVAL_INTERP); 142 CeedQFunctionAddInput(qfSetupRHS, "qdata", qdatasize, CEED_EVAL_NONE); 143 CeedQFunctionAddOutput(qfSetupRHS, "true_soln", ncompu, CEED_EVAL_NONE); 144 CeedQFunctionAddOutput(qfSetupRHS, "rhs", ncompu, CEED_EVAL_INTERP); 145 146 // Create the operator that builds the RHS and true solution 147 CeedOperatorCreate(ceed, qfSetupRHS, NULL, NULL, &opSetupRHS); 148 CeedOperatorSetField(opSetupRHS, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE); 149 CeedOperatorSetField(opSetupRHS, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED, 150 qdata); 151 CeedOperatorSetField(opSetupRHS, "true_soln", Erestrictui, 152 CEED_BASIS_COLLOCATED, *target); 153 CeedOperatorSetField(opSetupRHS, "rhs", Erestrictu, basisu, 154 CEED_VECTOR_ACTIVE); 155 156 // Set up the libCEED context 157 CeedQFunctionContext rhsSetupCtx; 158 CeedQFunctionContextCreate(ceed, &rhsSetupCtx); 159 CeedScalar rhsSetupData[2] = {R, l}; 160 CeedQFunctionContextSetData(rhsSetupCtx, CEED_MEM_HOST, CEED_COPY_VALUES, 161 sizeof rhsSetupData, &rhsSetupData); 162 CeedQFunctionSetContext(qfSetupRHS, rhsSetupCtx); 163 CeedQFunctionContextDestroy(&rhsSetupCtx); 164 165 // Setup RHS and target 166 CeedOperatorApply(opSetupRHS, xcoord, rhsceed, CEED_REQUEST_IMMEDIATE); 167 168 // Cleanup 169 CeedQFunctionDestroy(&qfSetupRHS); 170 CeedOperatorDestroy(&opSetupRHS); 171 } 172 173 // Cleanup 174 CeedQFunctionDestroy(&qfSetupGeo); 175 CeedOperatorDestroy(&opSetupGeo); 176 CeedVectorDestroy(&xcoord); 177 178 // Save libCEED data required for level 179 data->basisx = basisx; data->basisu = basisu; 180 data->Erestrictx = Erestrictx; 181 data->Erestrictu = Erestrictu; 182 data->Erestrictui = Erestrictui; 183 data->Erestrictqdi = Erestrictqdi; 184 data->qfApply = qfApply; 185 data->opApply = opApply; 186 data->qdata = qdata; 187 data->Xceed = Xceed; 188 data->Yceed = Yceed; 189 190 PetscFunctionReturn(0); 191 }; 192 193 // ----------------------------------------------------------------------------- 194 // Setup libCEED level transfer operator objects 195 // ----------------------------------------------------------------------------- 196 PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt numlevels, 197 CeedInt ncompu, 198 CeedData *data, CeedInt *leveldegrees, 199 CeedQFunction qfrestrict, CeedQFunction qfprolong) { 200 // Return early if numlevels=1 201 if (numlevels==1) 202 PetscFunctionReturn(0); 203 204 // Set up each level 205 for (CeedInt i=1; i<numlevels; i++) { 206 // P coarse and P fine 207 CeedInt Pc = leveldegrees[i-1] + 1; 208 CeedInt Pf = leveldegrees[i] + 1; 209 210 // Restriction - Fine to corse 211 CeedBasis basisctof; 212 CeedOperator opRestrict; 213 214 // Basis 215 CeedBasisCreateTensorH1Lagrange(ceed, 3, ncompu, Pc, Pf, 216 CEED_GAUSS_LOBATTO, &basisctof); 217 218 // Create the restriction operator 219 CeedOperatorCreate(ceed, qfrestrict, CEED_QFUNCTION_NONE, 220 CEED_QFUNCTION_NONE, &opRestrict); 221 CeedOperatorSetField(opRestrict, "input", data[i]->Erestrictu, 222 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 223 CeedOperatorSetField(opRestrict, "output", data[i-1]->Erestrictu, 224 basisctof, CEED_VECTOR_ACTIVE); 225 226 // Save libCEED data required for level 227 data[i]->basisctof = basisctof; 228 data[i]->opRestrict = opRestrict; 229 230 // Interpolation - Corse to fine 231 CeedOperator opProlong; 232 233 // Create the prolongation operator 234 CeedOperatorCreate(ceed, qfprolong, CEED_QFUNCTION_NONE, 235 CEED_QFUNCTION_NONE, &opProlong); 236 CeedOperatorSetField(opProlong, "input", data[i-1]->Erestrictu, 237 basisctof, CEED_VECTOR_ACTIVE); 238 CeedOperatorSetField(opProlong, "output", data[i]->Erestrictu, 239 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 240 241 // Save libCEED data required for level 242 data[i]->opProlong = opProlong; 243 } 244 245 PetscFunctionReturn(0); 246 }; 247 248 // ----------------------------------------------------------------------------- 249