Actual source code: ex2.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: static char help[] = "Tests L2 projection with DMSwarm using delta function particles.\n";

  3:  #include <petscdmplex.h>
  4:  #include <petscfe.h>
  5:  #include <petscdmswarm.h>
  6:  #include <petscds.h>
  7:  #include <petscksp.h>
  8:  #include <petsc/private/petscfeimpl.h>
  9: typedef struct {
 10:   PetscInt       dim;                              /* The topological mesh dimension */
 11:   PetscBool      simplex;                          /* Flag for simplices or tensor cells */
 12:   char           meshFilename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
 13:   PetscInt       faces;                            /* Number of faces per edge if unit square/cube generated */
 14:   PetscReal      domain_lo[3], domain_hi[3];       /* Lower left and upper right mesh corners */
 15:   DMBoundaryType boundary[3];                      /* The domain boundary type, e.g. periodic */
 16:   PetscInt       particlesPerCell;                 /* The number of partices per cell */
 17:   PetscReal      particleRelDx;                    /* Relative particle position perturbation compared to average cell diameter h */
 18:   PetscReal      meshRelDx;                        /* Relative vertex position perturbation compared to average cell diameter h */
 19:   PetscInt       k;                                /* Mode number for test function */
 20:   PetscReal      momentTol;                        /* Tolerance for checking moment conservation */
 21:   PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
 22: } AppCtx;

 24: /* const char *const ex2FunctionTypes[] = {"linear","x2_x4","sin","ex2FunctionTypes","EX2_FUNCTION_",0}; */

 26: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 27: {
 28:   AppCtx  *ctx = (AppCtx *) a_ctx;
 29:   PetscInt d;

 31:   u[0] = 0.0;
 32:   for (d = 0; d < dim; ++d) u[0] += x[d]/(ctx->domain_hi[d] - ctx->domain_lo[d]);
 33:   return 0;
 34: }

 36: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 37: {
 38:   AppCtx  *ctx = (AppCtx *) a_ctx;
 39:   PetscInt d;

 41:   u[0] = 1;
 42:   for (d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d])*PetscSqr(ctx->domain_hi[d]) - PetscPowRealInt(x[d], 4);
 43:   return 0;
 44: }

 46: static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 47: {
 48:   AppCtx *ctx = (AppCtx *) a_ctx;

 50:   u[0] = PetscSinScalar(2*PETSC_PI*ctx->k*x[0]/(ctx->domain_hi[0] - ctx->domain_lo[0]));
 51:   return 0;
 52: }

 54: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 55: {
 56:   PetscInt       ii, bd;
 57:   char           fstring[PETSC_MAX_PATH_LEN] = "linear";
 58:   PetscBool      flag;

 62:   options->dim              = 2;
 63:   options->simplex          = PETSC_TRUE;
 64:   options->faces            = 1;
 65:   options->domain_lo[0]     = 0.0;
 66:   options->domain_lo[1]     = 0.0;
 67:   options->domain_lo[2]     = 0.0;
 68:   options->domain_hi[0]     = 1.0;
 69:   options->domain_hi[1]     = 1.0;
 70:   options->domain_hi[2]     = 1.0;
 71:   options->boundary[0]      = DM_BOUNDARY_NONE; /* PERIODIC (plotting does not work in parallel, moments not conserved) */
 72:   options->boundary[1]      = DM_BOUNDARY_NONE;
 73:   options->boundary[2]      = DM_BOUNDARY_NONE;
 74:   options->particlesPerCell = 1;
 75:   options->k                = 1;
 76:   options->particleRelDx    = 1.e-20;
 77:   options->meshRelDx        = 1.e-20;
 78:   options->momentTol        = 100.*PETSC_MACHINE_EPSILON;
 79:   PetscStrcpy(options->meshFilename, "");

 81:   PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
 82:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex2.c", options->dim, &options->dim, NULL);
 83:   PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex2.c", options->simplex, &options->simplex, NULL);
 84:   PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex2.c", options->meshFilename, options->meshFilename, PETSC_MAX_PATH_LEN, NULL);
 85:   PetscOptionsInt("-faces", "Number of faces per edge if unit square/cube generated", "ex2.c", options->faces, &options->faces, NULL);
 86:   PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
 87:   PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
 88:   PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
 89:   PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
 90:   ii = options->dim;
 91:   PetscOptionsRealArray("-domain_hi", "Domain size", "ex2.c", options->domain_hi, &ii, NULL);
 92:   ii = options->dim;
 93:   PetscOptionsRealArray("-domain_lo", "Domain size", "ex2.c", options->domain_lo, &ii, NULL);
 94:   bd = options->boundary[0];
 95:   PetscOptionsEList("-x_boundary", "The x-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[0]], &bd, NULL);
 96:   options->boundary[0] = (DMBoundaryType) bd;
 97:   bd = options->boundary[1];
 98:   PetscOptionsEList("-y_boundary", "The y-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[1]], &bd, NULL);
 99:   options->boundary[1] = (DMBoundaryType) bd;
100:   bd = options->boundary[2];
101:   PetscOptionsEList("-z_boundary", "The z-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[2]], &bd, NULL);
102:   options->boundary[2] = (DMBoundaryType) bd;
103:   PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, PETSC_MAX_PATH_LEN, NULL);
104:   PetscStrcmp(fstring, "linear", &flag);
105:   if (flag) {
106:     options->func = linear;
107:   } else {
108:     PetscStrcmp(fstring, "sin", &flag);
109:     if (flag) {
110:       options->func = sinx;
111:     } else {
112:       PetscStrcmp(fstring, "x2_x4", &flag);
113:       options->func = x2_x4;
114:       if (!flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown function %s",fstring);
115:     }
116:   }
117:   PetscOptionsEnd();

119:   return(0);
120: }

122: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
123: {
124:   PetscRandom    rnd;
125:   PetscReal      interval = user->meshRelDx;
126:   Vec            coordinates;
127:   PetscScalar   *coords;
128:   PetscReal      *hh;
129:   PetscInt       d, cdim, N, p, bs;

133:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
134:   PetscRandomSetInterval(rnd, -interval, interval);
135:   PetscRandomSetFromOptions(rnd);
136:   DMGetCoordinatesLocal(dm, &coordinates);
137:   DMGetCoordinateDim(dm, &cdim);
138:   PetscCalloc1(PetscMax(user->dim,cdim),&hh);
139:   for (d = 0; d < user->dim; ++d) hh[d] = (user->domain_hi[d] - user->domain_lo[d])/user->faces;
140:   VecGetLocalSize(coordinates, &N);
141:   VecGetBlockSize(coordinates, &bs);
142:   if (bs != cdim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %D != %D", bs, cdim);
143:   VecGetArray(coordinates, &coords);
144:   for (p = 0; p < N; p += cdim) {
145:     PetscScalar *coord = &coords[p], value;

147:     for (d = 0; d < cdim; ++d) {
148:       PetscRandomGetValue(rnd, &value);
149:       coord[d] = PetscMax(user->domain_lo[d], PetscMin(user->domain_hi[d], PetscRealPart(coord[d] + value*hh[d])));
150:     }
151:   }
152:   VecRestoreArray(coordinates, &coords);
153:   PetscRandomDestroy(&rnd);
154:   PetscFree(hh);
155:   return(0);
156: }

158: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
159: {
160:   PetscBool      flg;

164:   PetscStrcmp(user->meshFilename, "", &flg);
165:   if (flg) {
166:     PetscInt faces[3];

168:     faces[0] = user->faces; faces[1] = user->faces; faces[2] = user->faces;
169:     DMPlexCreateBoxMesh(comm, user->dim, user->simplex, faces, user->domain_lo, user->domain_hi, user->boundary, PETSC_TRUE, dm);
170:   } else {
171:     DMPlexCreateFromFile(comm, user->meshFilename, PETSC_TRUE, dm);
172:     DMGetDimension(*dm, &user->dim);
173:   }
174:   {
175:     DM distributedMesh = NULL;

177:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
178:     if (distributedMesh) {
179:       DMDestroy(dm);
180:       *dm  = distributedMesh;
181:     }
182:   }
183:   DMLocalizeCoordinates(*dm); /* needed for periodic */
184:   DMSetFromOptions(*dm);
185:   PerturbVertices(*dm, user);
186:   PetscObjectSetName((PetscObject) *dm, "Mesh");
187:   DMViewFromOptions(*dm, NULL, "-dm_view");
188:   return(0);
189: }

191: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
192:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
193:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
194:                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
195: {
196:   g0[0] = 1.0;
197: }

199: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
200: {
201:   PetscFE        fe;
202:   PetscDS        ds;
203:   PetscInt       dim;

207:   DMGetDimension(dm, &dim);
208:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, user->simplex, NULL, -1, &fe);
209:   PetscObjectSetName((PetscObject) fe, "fe");
210:   DMSetField(dm, 0, NULL, (PetscObject) fe);
211:   DMCreateDS(dm);
212:   PetscFEDestroy(&fe);
213:   /* Setup to form mass matrix */
214:   DMGetDS(dm, &ds);
215:   PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
216:   return(0);
217: }

219: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
220: {
221:   PetscRandom    rnd, rndp;
222:   PetscReal      interval = user->particleRelDx;
223:   PetscScalar    value, *vals;
224:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
225:   PetscInt      *cellid;
226:   PetscInt       Ncell, Np = user->particlesPerCell, p, c, dim, d;

230:   DMGetDimension(dm, &dim);
231:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
232:   DMSetType(*sw, DMSWARM);
233:   DMSetDimension(*sw, dim);

235:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
236:   PetscRandomSetInterval(rnd, -1.0, 1.0);
237:   PetscRandomSetFromOptions(rnd);
238:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
239:   PetscRandomSetInterval(rndp, -interval, interval);
240:   PetscRandomSetFromOptions(rndp);

242:   DMSwarmSetType(*sw, DMSWARM_PIC);
243:   DMSwarmSetCellDM(*sw, dm);
244:   DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
245:   DMSwarmFinalizeFieldRegister(*sw);
246:   DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
247:   DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
248:   DMSetFromOptions(*sw);
249:   DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
250:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
251:   DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **) &vals);

253:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
254:   for (c = 0; c < Ncell; ++c) {
255:     if (Np == 1) {
256:       DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
257:       cellid[c] = c;
258:       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
259:     } else {
260:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
261:       for (d = 0; d < dim; ++d) xi0[d] = -1.0;
262:       for (p = 0; p < Np; ++p) {
263:         const PetscInt n   = c*Np + p;
264:         PetscReal      sum = 0.0, refcoords[3];

266:         cellid[n] = c;
267:         for (d = 0; d < dim; ++d) {PetscRandomGetValue(rnd, &value); refcoords[d] = PetscRealPart(value); sum += refcoords[d];}
268:         if (user->simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
269:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
270:       }
271:     }
272:   }
273:   PetscFree5(centroid, xi0, v0, J, invJ);
274:   for (c = 0; c < Ncell; ++c) {
275:     for (p = 0; p < Np; ++p) {
276:       const PetscInt n = c*Np + p;

278:       for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value); coords[n*dim+d] += PetscRealPart(value);}
279:       user->func(dim, 0.0, &coords[n*dim], 1, &vals[c], user);
280:     }
281:   }
282:   DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
283:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
284:   DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **) &vals);
285:   PetscRandomDestroy(&rnd);
286:   PetscRandomDestroy(&rndp);
287:   PetscObjectSetName((PetscObject) *sw, "Particles");
288:   DMViewFromOptions(*sw, NULL, "-sw_view");
289:   return(0);
290: }

292: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
293: {
294:   DM                 dm;
295:   const PetscReal   *coords;
296:   const PetscScalar *w;
297:   PetscReal          mom[3] = {0.0, 0.0, 0.0};
298:   PetscInt           cell, cStart, cEnd, dim;
299:   PetscErrorCode     ierr;

302:   DMGetDimension(sw, &dim);
303:   DMSwarmGetCellDM(sw, &dm);
304:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
305:   DMSwarmSortGetAccess(sw);
306:   DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
307:   DMSwarmGetField(sw, "w_q", NULL, NULL, (void **) &w);
308:   for (cell = cStart; cell < cEnd; ++cell) {
309:     PetscInt *pidx;
310:     PetscInt  Np, p, d;

312:     DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
313:     for (p = 0; p < Np; ++p) {
314:       const PetscInt   idx = pidx[p];
315:       const PetscReal *c   = &coords[idx*dim];

317:       mom[0] += PetscRealPart(w[idx]);
318:       mom[1] += PetscRealPart(w[idx]) * c[0];
319:       for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d]*c[d];
320:     }
321:     PetscFree(pidx);
322:   }
323:   DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
324:   DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &w);
325:   DMSwarmSortRestoreAccess(sw);
326:   MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) sw));
327:   return(0);
328: }

330: static void f0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux,
331:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
332:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
333:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
334: {
335:   f0[0] = u[0];
336: }

338: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
339:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
340:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
341:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
342: {
343:   f0[0] = x[0]*u[0];
344: }

346: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
347:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
348:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
349:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
350: {
351:   PetscInt d;

353:   f0[0] = 0.0;
354:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
355: }

357: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
358: {
359:   PetscDS        prob;
361:   PetscScalar    mom;

364:   DMGetDS(dm, &prob);
365:   PetscDSSetObjective(prob, 0, &f0_1);
366:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
367:   moments[0] = PetscRealPart(mom);
368:   PetscDSSetObjective(prob, 0, &f0_x);
369:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
370:   moments[1] = PetscRealPart(mom);
371:   PetscDSSetObjective(prob, 0, &f0_r2);
372:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
373:   moments[2] = PetscRealPart(mom);
374:   return(0);
375: }

377: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
378: {
379:   MPI_Comm       comm;
380:   KSP            ksp;
381:   Mat            M;            /* FEM mass matrix */
382:   Mat            M_p;          /* Particle mass matrix */
383:   Vec            f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
384:   PetscReal      pmoments[3];  /* \int f, \int x f, \int r^2 f */
385:   PetscReal      fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
386:   PetscInt       m;

390:   PetscObjectGetComm((PetscObject) dm, &comm);
391:   KSPCreate(comm, &ksp);
392:   KSPSetOptionsPrefix(ksp, "ptof_");
393:   DMGetGlobalVector(dm, &fhat);
394:   DMGetGlobalVector(dm, &rhs);

396:   DMCreateMassMatrix(sw, dm, &M_p);
397:   MatViewFromOptions(M_p, NULL, "-M_p_view");

399:   /* make particle weight vector */
400:   DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);

402:   /* create matrix RHS vector */
403:   MatMultTranspose(M_p, f, rhs);
404:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
405:   PetscObjectSetName((PetscObject) rhs,"rhs");
406:   VecViewFromOptions(rhs, NULL, "-rhs_view");

408:   DMCreateMatrix(dm, &M);
409:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
410:   MatViewFromOptions(M, NULL, "-M_view");
411:   KSPSetOperators(ksp, M, M);
412:   KSPSetFromOptions(ksp);
413:   KSPSolve(ksp, rhs, fhat);
414:   PetscObjectSetName((PetscObject) fhat,"fhat");
415:   VecViewFromOptions(fhat, NULL, "-fhat_view");

417:   /* Check moments of field */
418:   computeParticleMoments(sw, pmoments, user);
419:   computeFEMMoments(dm, fhat, fmoments, user);
420:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
421:   for (m = 0; m < 3; ++m) {
422:     if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
423:   }

425:   KSPDestroy(&ksp);
426:   MatDestroy(&M);
427:   MatDestroy(&M_p);
428:   DMRestoreGlobalVector(dm, &fhat);
429:   DMRestoreGlobalVector(dm, &rhs);

431:   return(0);
432: }


435: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
436: {

438:   MPI_Comm       comm;
439:   KSP            ksp;
440:   Mat            M;            /* FEM mass matrix */
441:   Mat            M_p;          /* Particle mass matrix */
442:   Vec            f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
443:   PetscReal      pmoments[3];  /* \int f, \int x f, \int r^2 f */
444:   PetscReal      fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
445:   PetscInt       m;

449:   PetscObjectGetComm((PetscObject) dm, &comm);

451:   KSPCreate(comm, &ksp);
452:   KSPSetOptionsPrefix(ksp, "ftop_");
453:   KSPSetFromOptions(ksp);

455:   DMGetGlobalVector(dm, &fhat);
456:   DMGetGlobalVector(dm, &rhs);

458:   DMCreateMassMatrix(sw, dm, &M_p);
459:   MatViewFromOptions(M_p, NULL, "-M_p_view");

461:   /* make particle weight vector */
462:   DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);

464:   /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
465:   PetscObjectSetName((PetscObject) rhs,"rhs");
466:   VecViewFromOptions(rhs, NULL, "-rhs_view");
467:   DMCreateMatrix(dm, &M);
468:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
469:   MatViewFromOptions(M, NULL, "-M_view");
470:   MatMultTranspose(M, fhat, rhs);

472:   KSPSetOperators(ksp, M_p, M_p);
473:   KSPSolveTranspose(ksp, rhs, f);
474:   PetscObjectSetName((PetscObject) fhat,"fhat");
475:   VecViewFromOptions(fhat, NULL, "-fhat_view");

477:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);

479:   /* Check moments */
480:   computeParticleMoments(sw, pmoments, user);
481:   computeFEMMoments(dm, fhat, fmoments, user);
482:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
483:   for (m = 0; m < 3; ++m) {
484:     if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
485:   }
486:   KSPDestroy(&ksp);
487:   MatDestroy(&M);
488:   MatDestroy(&M_p);
489:   DMRestoreGlobalVector(dm, &fhat);
490:   DMRestoreGlobalVector(dm, &rhs);
491:   return(0);
492: }

494: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
495: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC){

497:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
498:   PetscInt         debug = mesh->printFEM;
499:   DM               dmC;
500:   PetscSection     section;
501:   PetscQuadrature  quad = NULL;
502:   PetscScalar     *interpolant, *gradsum;
503:   PetscFEGeom      fegeom;
504:   PetscReal       *coords;
505:   const PetscReal *quadPoints, *quadWeights;
506:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
507:   PetscErrorCode   ierr;

510:   VecGetDM(locC, &dmC);
511:   VecSet(locC, 0.0);
512:   DMGetDimension(dm, &dim);
513:   DMGetCoordinateDim(dm, &coordDim);
514:   fegeom.dimEmbed = coordDim;
515:   DMGetLocalSection(dm, &section);
516:   PetscSectionGetNumFields(section, &numFields);
517:   for (field = 0; field < numFields; ++field) {
518:     PetscObject  obj;
519:     PetscClassId id;
520:     PetscInt     Nc;

522:     DMGetField(dm, field, NULL, &obj);
523:     PetscObjectGetClassId(obj, &id);
524:     if (id == PETSCFE_CLASSID) {
525:       PetscFE fe = (PetscFE) obj;

527:       PetscFEGetQuadrature(fe, &quad);
528:       PetscFEGetNumComponents(fe, &Nc);
529:     } else if (id == PETSCFV_CLASSID) {
530:       PetscFV fv = (PetscFV) obj;

532:       PetscFVGetQuadrature(fv, &quad);
533:       PetscFVGetNumComponents(fv, &Nc);
534:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
535:     numComponents += Nc;
536:   }
537:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
538:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
539:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
540:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
541:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
542:   for (v = vStart; v < vEnd; ++v) {
543:     PetscInt   *star = NULL;
544:     PetscInt    starSize, st, d, fc;

546:     PetscArrayzero(gradsum, coordDim*numComponents);
547:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
548:     for (st = 0; st < starSize*2; st += 2) {
549:       const PetscInt cell = star[st];
550:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
551:       PetscScalar   *x    = NULL;

553:       if ((cell < cStart) || (cell >= cEnd)) continue;
554:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
555:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
556:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
557:         PetscObject  obj;
558:         PetscClassId id;
559:         PetscInt     Nb, Nc, q, qc = 0;

561:         PetscArrayzero(grad, coordDim*numComponents);
562:         DMGetField(dm, field, NULL, &obj);
563:         PetscObjectGetClassId(obj, &id);
564:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
565:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
566:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
567:         for (q = 0; q < Nq; ++q) {
568:           if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], cell, q);
569:           if (ierr) {
570:             PetscErrorCode ierr2;
571:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
572:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
573:             ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
574:             
575:           }
576:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &fegeom, q, interpolant);}
577:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
578:           for (fc = 0; fc < Nc; ++fc) {
579:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

581:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
582:           }
583:         }
584:         fieldOffset += Nb;
585:       }
586:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
587:       for (fc = 0; fc < numComponents; ++fc) {
588:         for (d = 0; d < coordDim; ++d) {
589:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
590:         }
591:       }
592:       if (debug) {
593:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
594:         for (fc = 0; fc < numComponents; ++fc) {
595:           for (d = 0; d < coordDim; ++d) {
596:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
597:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
598:           }
599:         }
600:         PetscPrintf(PETSC_COMM_SELF, "]\n");
601:       }
602:     }
603:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
604:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
605:   }
606:   PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
607:   return(0);
608: }

610: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
611: {

613:   MPI_Comm       comm;
614:   KSP            ksp;
615:   Mat            M;                   /* FEM mass matrix */
616:   Mat            M_p;                 /* Particle mass matrix */
617:   Vec            f, rhs, fhat, grad;  /* Particle field f, \int phi_i f, FEM field */
618:   PetscReal      pmoments[3];         /* \int f, \int x f, \int r^2 f */
619:   PetscReal      fmoments[3];         /* \int \hat f, \int x \hat f, \int r^2 \hat f */
620:   PetscInt       m;

624:   PetscObjectGetComm((PetscObject) dm, &comm);
625:   KSPCreate(comm, &ksp);
626:   KSPSetOptionsPrefix(ksp, "ptof_");
627:   DMGetGlobalVector(dm, &fhat);
628:   DMGetGlobalVector(dm, &rhs);
629:   DMGetGlobalVector(dm, &grad);

631:   DMCreateMassMatrix(sw, dm, &M_p);
632:   MatViewFromOptions(M_p, NULL, "-M_p_view");

634:   /* make particle weight vector */
635:   DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);

637:   /* create matrix RHS vector */
638:   MatMultTranspose(M_p, f, rhs);
639:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
640:   PetscObjectSetName((PetscObject) rhs,"rhs");
641:   VecViewFromOptions(rhs, NULL, "-rhs_view");

643:   DMCreateMatrix(dm, &M);
644:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);

646:   InterpolateGradient(dm, fhat, grad);

648:   MatViewFromOptions(M, NULL, "-M_view");
649:   KSPSetOperators(ksp, M, M);
650:   KSPSetFromOptions(ksp);
651:   KSPSolve(ksp, rhs, grad);
652:   PetscObjectSetName((PetscObject) fhat,"fhat");
653:   VecViewFromOptions(fhat, NULL, "-fhat_view");

655:   /* Check moments of field */
656:   computeParticleMoments(sw, pmoments, user);
657:   computeFEMMoments(dm, grad, fmoments, user);
658:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
659:   for (m = 0; m < 3; ++m) {
660:     if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
661:   }

663:   KSPDestroy(&ksp);
664:   MatDestroy(&M);
665:   MatDestroy(&M_p);
666:   DMRestoreGlobalVector(dm, &fhat);
667:   DMRestoreGlobalVector(dm, &rhs);
668:   DMRestoreGlobalVector(dm, &grad);

670:   return(0);
671: }

673: int main (int argc, char * argv[]) {
674:   MPI_Comm       comm;
675:   DM             dm, sw;
676:   AppCtx         user;

679:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
680:   comm = PETSC_COMM_WORLD;
681:   ProcessOptions(comm, &user);
682:   CreateMesh(comm, &dm, &user);
683:   CreateFEM(dm, &user);
684:   CreateParticles(dm, &sw, &user);
685:   TestL2ProjectionParticlesToField(dm, sw, &user);
686:   TestL2ProjectionFieldToParticles(dm, sw, &user);
687:   TestFieldGradientProjection(dm, sw, &user);
688:   DMDestroy(&dm);
689:   DMDestroy(&sw);
690:   PetscFinalize();
691:   return ierr;
692: }


695: /*TEST

697:   test:
698:     suffix: proj_tri_0
699:     requires: triangle !complex
700:     args: -dim 2 -faces 1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
701:     filter: grep -v marker | grep -v atomic | grep -v usage

703:   test:
704:     suffix: proj_tri_2_faces
705:     requires: triangle !complex
706:     args: -dim 2 -faces 2  -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
707:     filter: grep -v marker | grep -v atomic | grep -v usage

709:   test:
710:     suffix: proj_quad_0
711:     requires: triangle !complex
712:     args: -dim 2 -simplex 0 -faces 1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
713:     filter: grep -v marker | grep -v atomic | grep -v usage

715:   test:
716:     suffix: proj_quad_2_faces
717:     requires: triangle !complex
718:     args: -dim 2 -simplex 0 -faces 2 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
719:     filter: grep -v marker | grep -v atomic | grep -v usage

721:   test:
722:     suffix: proj_tri_5P
723:     requires: triangle !complex
724:     args: -dim 2 -faces 1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
725:     filter: grep -v marker | grep -v atomic | grep -v usage

727:   test:
728:     suffix: proj_quad_5P
729:     requires: triangle !complex
730:     args: -dim 2 -simplex 0 -faces 1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
731:     filter: grep -v marker | grep -v atomic | grep -v usage

733:   test:
734:     suffix: proj_tri_mdx
735:     requires: triangle !complex
736:     args: -dim 2 -faces 1 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
737:     filter: grep -v marker | grep -v atomic | grep -v usage

739:   test:
740:     suffix: proj_tri_mdx_5P
741:     requires: triangle !complex
742:     args: -dim 2 -faces 1 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
743:     filter: grep -v marker | grep -v atomic | grep -v usage

745:   test:
746:     suffix: proj_tri_3d
747:     requires: ctetgen !complex
748:     args: -dim 3 -faces 1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
749:     filter: grep -v marker | grep -v atomic | grep -v usage

751:   test:
752:     suffix: proj_tri_3d_2_faces
753:     requires: ctetgen !complex
754:     args: -dim 3 -faces 2 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
755:     filter: grep -v marker | grep -v atomic | grep -v usage

757:   test:
758:     suffix: proj_tri_3d_5P
759:     requires: ctetgen !complex
760:     args: -dim 3 -faces 1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
761:     filter: grep -v marker | grep -v atomic | grep -v usage

763:   test:
764:     suffix: proj_tri_3d_mdx
765:     requires: ctetgen !complex
766:     args: -dim 3 -faces 1 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
767:     filter: grep -v marker | grep -v atomic | grep -v usage

769:   test:
770:     suffix: proj_tri_3d_mdx_5P
771:     requires: ctetgen !complex
772:     args: -dim 3 -faces 1 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
773:     filter: grep -v marker | grep -v atomic | grep -v usage

775:   test:
776:     suffix: proj_tri_3d_mdx_2_faces
777:     requires: ctetgen !complex
778:     args: -dim 3 -faces 2 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
779:     filter: grep -v marker | grep -v atomic | grep -v usage

781:   test:
782:     suffix: proj_tri_3d_mdx_5P_2_faces
783:     requires: ctetgen !complex
784:     args: -dim 3 -faces 2 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu  -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
785:     filter: grep -v marker | grep -v atomic | grep -v usage

787: TEST*/