Actual source code: ex2.c

  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:   PetscReal L[3];                             /* Dimensions of the mesh bounding box */
 11:   PetscInt  particlesPerCell;                 /* The number of partices per cell */
 12:   PetscReal particleRelDx;                    /* Relative particle position perturbation compared to average cell diameter h */
 13:   PetscReal meshRelDx;                        /* Relative vertex position perturbation compared to average cell diameter h */
 14:   PetscInt  k;                                /* Mode number for test function */
 15:   PetscReal momentTol;                        /* Tolerance for checking moment conservation */
 16:   PetscBool useBlockDiagPrec;                 /* Use the block diagonal of the normal equations as a preconditioner */
 17:   PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
 18: } AppCtx;

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

 22: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 23: {
 24:   AppCtx  *ctx = (AppCtx *) a_ctx;
 25:   PetscInt d;

 27:   u[0] = 0.0;
 28:   for (d = 0; d < dim; ++d) u[0] += x[d]/(ctx->L[d]);
 29:   return 0;
 30: }

 32: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 33: {
 34:   AppCtx  *ctx = (AppCtx *) a_ctx;
 35:   PetscInt d;

 37:   u[0] = 1;
 38:   for (d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d])*PetscSqr(ctx->L[d]) - PetscPowRealInt(x[d], 4);
 39:   return 0;
 40: }

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

 46:   u[0] = PetscSinScalar(2*PETSC_PI*ctx->k*x[0]/(ctx->L[0]));
 47:   return 0;
 48: }

 50: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 51: {
 52:   char           fstring[PETSC_MAX_PATH_LEN] = "linear";
 53:   PetscBool      flag;

 57:   options->particlesPerCell = 1;
 58:   options->k                = 1;
 59:   options->particleRelDx    = 1.e-20;
 60:   options->meshRelDx        = 1.e-20;
 61:   options->momentTol        = 100.*PETSC_MACHINE_EPSILON;
 62:   options->useBlockDiagPrec = PETSC_FALSE;

 64:   PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
 65:   PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
 66:   PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
 67:   PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
 68:   PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
 69:   PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL);
 70:   PetscStrcmp(fstring, "linear", &flag);
 71:   if (flag) {
 72:     options->func = linear;
 73:   } else {
 74:     PetscStrcmp(fstring, "sin", &flag);
 75:     if (flag) {
 76:       options->func = sinx;
 77:     } else {
 78:       PetscStrcmp(fstring, "x2_x4", &flag);
 79:       options->func = x2_x4;
 81:     }
 82:   }
 83:   PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL);
 84:   PetscOptionsBool("-block_diag_prec", "Use the block diagonal of the normal equations to precondition the particle projection", "ex2.c", options->useBlockDiagPrec, &options->useBlockDiagPrec, NULL);
 85:   PetscOptionsEnd();

 87:   return 0;
 88: }

 90: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
 91: {
 92:   PetscRandom    rnd;
 93:   PetscReal      interval = user->meshRelDx;
 94:   Vec            coordinates;
 95:   PetscScalar   *coords;
 96:   PetscReal     *hh, low[3], high[3];
 97:   PetscInt       d, cdim, cEnd, N, p, bs;

100:   DMGetBoundingBox(dm, low, high);
101:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
102:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
103:   PetscRandomSetInterval(rnd, -interval, interval);
104:   PetscRandomSetFromOptions(rnd);
105:   DMGetCoordinatesLocal(dm, &coordinates);
106:   DMGetCoordinateDim(dm, &cdim);
107:   PetscCalloc1(cdim,&hh);
108:   for (d = 0; d < cdim; ++d) hh[d] = (user->L[d])/PetscPowReal(cEnd, 1./cdim);
109:   VecGetLocalSize(coordinates, &N);
110:   VecGetBlockSize(coordinates, &bs);
112:   VecGetArray(coordinates, &coords);
113:   for (p = 0; p < N; p += cdim) {
114:     PetscScalar *coord = &coords[p], value;

116:     for (d = 0; d < cdim; ++d) {
117:       PetscRandomGetValue(rnd, &value);
118:       coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value*hh[d])));
119:     }
120:   }
121:   VecRestoreArray(coordinates, &coords);
122:   PetscRandomDestroy(&rnd);
123:   PetscFree(hh);
124:   return 0;
125: }

127: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
128: {
129:   PetscReal      low[3], high[3];
130:   PetscInt       cdim, d;

133:   DMCreate(comm, dm);
134:   DMSetType(*dm, DMPLEX);
135:   DMSetFromOptions(*dm);
136:   DMViewFromOptions(*dm, NULL, "-dm_view");

138:   DMGetCoordinateDim(*dm, &cdim);
139:   DMGetBoundingBox(*dm, low, high);
140:   for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
141:   PerturbVertices(*dm, user);
142:   return 0;
143: }

145: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
146:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
147:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
148:                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
149: {
150:   g0[0] = 1.0;
151: }

153: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
154: {
155:   PetscFE        fe;
156:   PetscDS        ds;
157:   DMPolytopeType ct;
158:   PetscBool      simplex;
159:   PetscInt       dim, cStart;

162:   DMGetDimension(dm, &dim);
163:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
164:   DMPlexGetCellType(dm, cStart, &ct);
165:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
166:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, NULL, -1, &fe);
167:   PetscObjectSetName((PetscObject) fe, "fe");
168:   DMSetField(dm, 0, NULL, (PetscObject) fe);
169:   DMCreateDS(dm);
170:   PetscFEDestroy(&fe);
171:   /* Setup to form mass matrix */
172:   DMGetDS(dm, &ds);
173:   PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
174:   return 0;
175: }

177: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
178: {
179:   PetscRandom    rnd, rndp;
180:   PetscReal      interval = user->particleRelDx;
181:   DMPolytopeType ct;
182:   PetscBool      simplex;
183:   PetscScalar    value, *vals;
184:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
185:   PetscInt      *cellid;
186:   PetscInt       Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;

189:   DMGetDimension(dm, &dim);
190:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
191:   DMPlexGetCellType(dm, cStart, &ct);
192:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
193:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
194:   DMSetType(*sw, DMSWARM);
195:   DMSetDimension(*sw, dim);

197:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
198:   PetscRandomSetInterval(rnd, -1.0, 1.0);
199:   PetscRandomSetFromOptions(rnd);
200:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
201:   PetscRandomSetInterval(rndp, -interval, interval);
202:   PetscRandomSetFromOptions(rndp);

204:   DMSwarmSetType(*sw, DMSWARM_PIC);
205:   DMSwarmSetCellDM(*sw, dm);
206:   DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
207:   DMSwarmFinalizeFieldRegister(*sw);
208:   DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
209:   DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
210:   DMSetFromOptions(*sw);
211:   DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
212:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
213:   DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **) &vals);

215:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
216:   for (c = 0; c < Ncell; ++c) {
217:     if (Np == 1) {
218:       DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
219:       cellid[c] = c;
220:       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
221:     } else {
222:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
223:       for (d = 0; d < dim; ++d) xi0[d] = -1.0;
224:       for (p = 0; p < Np; ++p) {
225:         const PetscInt n   = c*Np + p;
226:         PetscReal      sum = 0.0, refcoords[3];

228:         cellid[n] = c;
229:         for (d = 0; d < dim; ++d) {PetscRandomGetValue(rnd, &value)); refcoords[d] = PetscRealPart(value; sum += refcoords[d];}
230:         if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
231:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
232:       }
233:     }
234:   }
235:   PetscFree5(centroid, xi0, v0, J, invJ);
236:   for (c = 0; c < Ncell; ++c) {
237:     for (p = 0; p < Np; ++p) {
238:       const PetscInt n = c*Np + p;

240:       for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value)); coords[n*dim+d] += PetscRealPart(value;}
241:       user->func(dim, 0.0, &coords[n*dim], 1, &vals[c], user);
242:     }
243:   }
244:   DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
245:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
246:   DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **) &vals);
247:   PetscRandomDestroy(&rnd);
248:   PetscRandomDestroy(&rndp);
249:   PetscObjectSetName((PetscObject) *sw, "Particles");
250:   DMViewFromOptions(*sw, NULL, "-sw_view");
251:   return 0;
252: }

254: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
255: {
256:   DM                 dm;
257:   const PetscReal   *coords;
258:   const PetscScalar *w;
259:   PetscReal          mom[3] = {0.0, 0.0, 0.0};
260:   PetscInt           cell, cStart, cEnd, dim;

263:   DMGetDimension(sw, &dim);
264:   DMSwarmGetCellDM(sw, &dm);
265:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
266:   DMSwarmSortGetAccess(sw);
267:   DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
268:   DMSwarmGetField(sw, "w_q", NULL, NULL, (void **) &w);
269:   for (cell = cStart; cell < cEnd; ++cell) {
270:     PetscInt *pidx;
271:     PetscInt  Np, p, d;

273:     DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
274:     for (p = 0; p < Np; ++p) {
275:       const PetscInt   idx = pidx[p];
276:       const PetscReal *c   = &coords[idx*dim];

278:       mom[0] += PetscRealPart(w[idx]);
279:       mom[1] += PetscRealPart(w[idx]) * c[0];
280:       for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d]*c[d];
281:     }
282:     PetscFree(pidx);
283:   }
284:   DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
285:   DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &w);
286:   DMSwarmSortRestoreAccess(sw);
287:   MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) sw));
288:   return 0;
289: }

291: static void f0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux,
292:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
293:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
294:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
295: {
296:   f0[0] = u[0];
297: }

299: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
300:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
301:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
302:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
303: {
304:   f0[0] = x[0]*u[0];
305: }

307: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
308:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
309:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
310:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
311: {
312:   PetscInt d;

314:   f0[0] = 0.0;
315:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
316: }

318: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
319: {
320:   PetscDS        prob;
321:   PetscScalar    mom;

324:   DMGetDS(dm, &prob);
325:   PetscDSSetObjective(prob, 0, &f0_1);
326:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
327:   moments[0] = PetscRealPart(mom);
328:   PetscDSSetObjective(prob, 0, &f0_x);
329:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
330:   moments[1] = PetscRealPart(mom);
331:   PetscDSSetObjective(prob, 0, &f0_r2);
332:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
333:   moments[2] = PetscRealPart(mom);
334:   return 0;
335: }

337: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
338: {
339:   MPI_Comm       comm;
340:   KSP            ksp;
341:   Mat            M;            /* FEM mass matrix */
342:   Mat            M_p;          /* Particle mass matrix */
343:   Vec            f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
344:   PetscReal      pmoments[3];  /* \int f, \int x f, \int r^2 f */
345:   PetscReal      fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
346:   PetscInt       m;

349:   PetscObjectGetComm((PetscObject) dm, &comm);
350:   KSPCreate(comm, &ksp);
351:   KSPSetOptionsPrefix(ksp, "ptof_");
352:   DMGetGlobalVector(dm, &fhat);
353:   DMGetGlobalVector(dm, &rhs);

355:   DMCreateMassMatrix(sw, dm, &M_p);
356:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

361:   /* create matrix RHS vector */
362:   MatMultTranspose(M_p, f, rhs);
363:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
364:   PetscObjectSetName((PetscObject) rhs,"rhs");
365:   VecViewFromOptions(rhs, NULL, "-rhs_view");

367:   DMCreateMatrix(dm, &M);
368:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
369:   MatViewFromOptions(M, NULL, "-M_view");
370:   KSPSetOperators(ksp, M, M);
371:   KSPSetFromOptions(ksp);
372:   KSPSolve(ksp, rhs, fhat);
373:   PetscObjectSetName((PetscObject) fhat,"fhat");
374:   VecViewFromOptions(fhat, NULL, "-fhat_view");

376:   /* Check moments of field */
377:   computeParticleMoments(sw, pmoments, user);
378:   computeFEMMoments(dm, fhat, fmoments, user);
379:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
380:   for (m = 0; m < 3; ++m) {
382:   }

384:   KSPDestroy(&ksp);
385:   MatDestroy(&M);
386:   MatDestroy(&M_p);
387:   DMRestoreGlobalVector(dm, &fhat);
388:   DMRestoreGlobalVector(dm, &rhs);

390:   return 0;
391: }

393: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
394: {

396:   MPI_Comm       comm;
397:   KSP            ksp;
398:   Mat            M;            /* FEM mass matrix */
399:   Mat            M_p, PM_p;    /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
400:   Vec            f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
401:   PetscReal      pmoments[3];  /* \int f, \int x f, \int r^2 f */
402:   PetscReal      fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
403:   PetscInt       m;

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

408:   KSPCreate(comm, &ksp);
409:   KSPSetOptionsPrefix(ksp, "ftop_");
410:   KSPSetFromOptions(ksp);

412:   DMGetGlobalVector(dm, &fhat);
413:   DMGetGlobalVector(dm, &rhs);

415:   DMCreateMassMatrix(sw, dm, &M_p);
416:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

421:   /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
422:   PetscObjectSetName((PetscObject) rhs,"rhs");
423:   VecViewFromOptions(rhs, NULL, "-rhs_view");
424:   DMCreateMatrix(dm, &M);
425:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
426:   MatViewFromOptions(M, NULL, "-M_view");
427:   MatMultTranspose(M, fhat, rhs);
428:   if (user->useBlockDiagPrec) DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);
429:   else                        {PetscObjectReference((PetscObject) M_p); PM_p = M_p;}

431:   KSPSetOperators(ksp, M_p, PM_p);
432:   KSPSolveTranspose(ksp, rhs, f);
433:   PetscObjectSetName((PetscObject) fhat,"fhat");
434:   VecViewFromOptions(fhat, NULL, "-fhat_view");

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

438:   /* Check moments */
439:   computeParticleMoments(sw, pmoments, user);
440:   computeFEMMoments(dm, fhat, fmoments, user);
441:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
442:   for (m = 0; m < 3; ++m) {
444:   }
445:   KSPDestroy(&ksp);
446:   MatDestroy(&M);
447:   MatDestroy(&M_p);
448:   MatDestroy(&PM_p);
449:   DMRestoreGlobalVector(dm, &fhat);
450:   DMRestoreGlobalVector(dm, &rhs);
451:   return 0;
452: }

454: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
455: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
456: {
457:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
458:   PetscInt         debug = mesh->printFEM;
459:   DM               dmC;
460:   PetscSection     section;
461:   PetscQuadrature  quad = NULL;
462:   PetscScalar     *interpolant, *gradsum;
463:   PetscFEGeom      fegeom;
464:   PetscReal       *coords;
465:   const PetscReal *quadPoints, *quadWeights;
466:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;

468:   VecGetDM(locC, &dmC);
469:   VecSet(locC, 0.0);
470:   DMGetDimension(dm, &dim);
471:   DMGetCoordinateDim(dm, &coordDim);
472:   fegeom.dimEmbed = coordDim;
473:   DMGetLocalSection(dm, &section);
474:   PetscSectionGetNumFields(section, &numFields);
475:   for (field = 0; field < numFields; ++field) {
476:     PetscObject  obj;
477:     PetscClassId id;
478:     PetscInt     Nc;

480:     DMGetField(dm, field, NULL, &obj);
481:     PetscObjectGetClassId(obj, &id);
482:     if (id == PETSCFE_CLASSID) {
483:       PetscFE fe = (PetscFE) obj;

485:       PetscFEGetQuadrature(fe, &quad);
486:       PetscFEGetNumComponents(fe, &Nc);
487:     } else if (id == PETSCFV_CLASSID) {
488:       PetscFV fv = (PetscFV) obj;

490:       PetscFVGetQuadrature(fv, &quad);
491:       PetscFVGetNumComponents(fv, &Nc);
492:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
493:     numComponents += Nc;
494:   }
495:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
497:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
498:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
499:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
500:   for (v = vStart; v < vEnd; ++v) {
501:     PetscInt *star = NULL;
502:     PetscInt starSize, st, d, fc;

504:     PetscArrayzero(gradsum, coordDim*numComponents);
505:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
506:     for (st = 0; st < starSize*2; st += 2) {
507:       const PetscInt cell = star[st];
508:       PetscScalar    *grad = &gradsum[coordDim*numComponents];
509:       PetscScalar    *x    = NULL;

511:       if ((cell < cStart) || (cell >= cEnd)) continue;
512:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
513:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
514:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
515:         PetscObject  obj;
516:         PetscClassId id;
517:         PetscInt     Nb, Nc, q, qc = 0;

519:         PetscArrayzero(grad, coordDim*numComponents);
520:         DMGetField(dm, field, NULL, &obj);
521:         PetscObjectGetClassId(obj, &id);
522:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc));PetscCall(PetscFEGetDimension((PetscFE) obj, &Nb);}
523:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
524:         else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
525:         for (q = 0; q < Nq; ++q) {
527:           if (id == PETSCFE_CLASSID)      PetscFEInterpolateGradient_Static((PetscFE) obj, 1, &x[fieldOffset], &fegeom, q, interpolant);
528:           else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
529:           for (fc = 0; fc < Nc; ++fc) {
530:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

532:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
533:           }
534:         }
535:         fieldOffset += Nb;
536:       }
537:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
538:       for (fc = 0; fc < numComponents; ++fc) {
539:         for (d = 0; d < coordDim; ++d) {
540:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
541:         }
542:       }
543:       if (debug) {
544:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
545:         for (fc = 0; fc < numComponents; ++fc) {
546:           for (d = 0; d < coordDim; ++d) {
547:             if (fc || d > 0) PetscPrintf(PETSC_COMM_SELF, ", ");
548:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
549:           }
550:         }
551:         PetscPrintf(PETSC_COMM_SELF, "]\n");
552:       }
553:     }
554:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
555:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
556:   }
557:   PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
558:   return 0;
559: }

561: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
562: {

564:   MPI_Comm       comm;
565:   KSP            ksp;
566:   Mat            M;                   /* FEM mass matrix */
567:   Mat            M_p;                 /* Particle mass matrix */
568:   Vec            f, rhs, fhat, grad;  /* Particle field f, \int phi_i f, FEM field */
569:   PetscReal      pmoments[3];         /* \int f, \int x f, \int r^2 f */
570:   PetscReal      fmoments[3];         /* \int \hat f, \int x \hat f, \int r^2 \hat f */
571:   PetscInt       m;

574:   PetscObjectGetComm((PetscObject) dm, &comm);
575:   KSPCreate(comm, &ksp);
576:   KSPSetOptionsPrefix(ksp, "ptof_");
577:   DMGetGlobalVector(dm, &fhat);
578:   DMGetGlobalVector(dm, &rhs);
579:   DMGetGlobalVector(dm, &grad);

581:   DMCreateMassMatrix(sw, dm, &M_p);
582:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

587:   /* create matrix RHS vector */
588:   MatMultTranspose(M_p, f, rhs);
589:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
590:   PetscObjectSetName((PetscObject) rhs,"rhs");
591:   VecViewFromOptions(rhs, NULL, "-rhs_view");

593:   DMCreateMatrix(dm, &M);
594:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);

596:   InterpolateGradient(dm, fhat, grad);

598:   MatViewFromOptions(M, NULL, "-M_view");
599:   KSPSetOperators(ksp, M, M);
600:   KSPSetFromOptions(ksp);
601:   KSPSolve(ksp, rhs, grad);
602:   PetscObjectSetName((PetscObject) fhat,"fhat");
603:   VecViewFromOptions(fhat, NULL, "-fhat_view");

605:   /* Check moments of field */
606:   computeParticleMoments(sw, pmoments, user);
607:   computeFEMMoments(dm, grad, fmoments, user);
608:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
609:   for (m = 0; m < 3; ++m) {
611:   }

613:   KSPDestroy(&ksp);
614:   MatDestroy(&M);
615:   MatDestroy(&M_p);
616:   DMRestoreGlobalVector(dm, &fhat);
617:   DMRestoreGlobalVector(dm, &rhs);
618:   DMRestoreGlobalVector(dm, &grad);

620:   return 0;
621: }

623: int main (int argc, char *argv[])
624: {
625:   MPI_Comm       comm;
626:   DM             dm, sw;
627:   AppCtx         user;

629:   PetscInitialize(&argc, &argv, NULL, help);
630:   comm = PETSC_COMM_WORLD;
631:   ProcessOptions(comm, &user);
632:   CreateMesh(comm, &dm, &user);
633:   CreateFEM(dm, &user);
634:   CreateParticles(dm, &sw, &user);
635:   TestL2ProjectionParticlesToField(dm, sw, &user);
636:   TestL2ProjectionFieldToParticles(dm, sw, &user);
637:   TestFieldGradientProjection(dm, sw, &user);
638:   DMDestroy(&dm);
639:   DMDestroy(&sw);
640:   PetscFinalize();
641:   return 0;
642: }

644: /*TEST

646:   # Swarm does not handle complex or quad
647:   build:
648:     requires: !complex double

650:   test:
651:     suffix: proj_tri_0
652:     requires: triangle
653:     args: -dm_plex_box_faces 1,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
654:     filter: grep -v marker | grep -v atomic | grep -v usage

656:   test:
657:     suffix: proj_tri_2_faces
658:     requires: triangle
659:     args: -dm_plex_box_faces 2,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
660:     filter: grep -v marker | grep -v atomic | grep -v usage

662:   test:
663:     suffix: proj_quad_0
664:     requires: triangle
665:     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,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
666:     filter: grep -v marker | grep -v atomic | grep -v usage

668:   test:
669:     suffix: proj_quad_2_faces
670:     requires: triangle
671:     args: -dm_plex_simplex 0 -dm_plex_box_faces 2,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
672:     filter: grep -v marker | grep -v atomic | grep -v usage

674:   test:
675:     suffix: proj_tri_5P
676:     requires: triangle
677:     args: -dm_plex_box_faces 1,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
678:     filter: grep -v marker | grep -v atomic | grep -v usage

680:   test:
681:     suffix: proj_quad_5P
682:     requires: triangle
683:     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,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
684:     filter: grep -v marker | grep -v atomic | grep -v usage

686:   test:
687:     suffix: proj_tri_mdx
688:     requires: triangle
689:     args: -dm_plex_box_faces 1,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
690:     filter: grep -v marker | grep -v atomic | grep -v usage

692:   test:
693:     suffix: proj_tri_mdx_5P
694:     requires: triangle
695:     args: -dm_plex_box_faces 1,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
696:     filter: grep -v marker | grep -v atomic | grep -v usage

698:   test:
699:     suffix: proj_tri_3d
700:     requires: ctetgen
701:     args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,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
702:     filter: grep -v marker | grep -v atomic | grep -v usage

704:   test:
705:     suffix: proj_tri_3d_2_faces
706:     requires: ctetgen
707:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,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
708:     filter: grep -v marker | grep -v atomic | grep -v usage

710:   test:
711:     suffix: proj_tri_3d_5P
712:     requires: ctetgen
713:     args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,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
714:     filter: grep -v marker | grep -v atomic | grep -v usage

716:   test:
717:     suffix: proj_tri_3d_mdx
718:     requires: ctetgen
719:     args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,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
720:     filter: grep -v marker | grep -v atomic | grep -v usage

722:   test:
723:     suffix: proj_tri_3d_mdx_5P
724:     requires: ctetgen
725:     args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,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
726:     filter: grep -v marker | grep -v atomic | grep -v usage

728:   test:
729:     suffix: proj_tri_3d_mdx_2_faces
730:     requires: ctetgen
731:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,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
732:     filter: grep -v marker | grep -v atomic | grep -v usage

734:   test:
735:     suffix: proj_tri_3d_mdx_5P_2_faces
736:     requires: ctetgen
737:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,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
738:     filter: grep -v marker | grep -v atomic | grep -v usage

740:   test:
741:     suffix: proj_quad_lsqr_scale
742:     requires: !complex
743:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
744:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
745:       -particlesPerCell 200 \
746:       -ptof_pc_type lu  \
747:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none

749:   test:
750:     suffix: proj_quad_lsqr_prec_scale
751:     requires: !complex
752:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
753:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
754:       -particlesPerCell 200 \
755:       -ptof_pc_type lu  \
756:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type lu -ftop_pc_factor_shift_type nonzero -block_diag_prec

758: TEST*/