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;

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

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

 86:   return 0;
 87: }

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

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

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

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

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

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

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

149: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
150: {
151:   PetscFE        fe;
152:   PetscDS        ds;
153:   DMPolytopeType ct;
154:   PetscBool      simplex;
155:   PetscInt       dim, cStart;

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

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

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

193:   PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd);
194:   PetscRandomSetInterval(rnd, -1.0, 1.0);
195:   PetscRandomSetFromOptions(rnd);
196:   PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rndp);
197:   PetscRandomSetInterval(rndp, -interval, interval);
198:   PetscRandomSetFromOptions(rndp);

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

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

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

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

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

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

277:     DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
278:     for (p = 0; p < Np; ++p) {
279:       const PetscInt   idx = pidx[p];
280:       const PetscReal *c   = &coords[idx * dim];

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

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

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

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

309:   f0[0] = 0.0;
310:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
311: }

313: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
314: {
315:   PetscDS     prob;
316:   PetscScalar mom;

319:   DMGetDS(dm, &prob);
320:   PetscDSSetObjective(prob, 0, &f0_1);
321:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
322:   moments[0] = PetscRealPart(mom);
323:   PetscDSSetObjective(prob, 0, &f0_x);
324:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
325:   moments[1] = PetscRealPart(mom);
326:   PetscDSSetObjective(prob, 0, &f0_r2);
327:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
328:   moments[2] = PetscRealPart(mom);
329:   return 0;
330: }

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

344:   PetscObjectGetComm((PetscObject)dm, &comm);
345:   KSPCreate(comm, &ksp);
346:   KSPSetOptionsPrefix(ksp, "ptof_");
347:   DMGetGlobalVector(dm, &fhat);
348:   DMGetGlobalVector(dm, &rhs);

350:   DMCreateMassMatrix(sw, dm, &M_p);
351:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

356:   /* create matrix RHS vector */
357:   MatMultTranspose(M_p, f, rhs);
358:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
359:   PetscObjectSetName((PetscObject)rhs, "rhs");
360:   VecViewFromOptions(rhs, NULL, "-rhs_view");

362:   DMCreateMatrix(dm, &M);
363:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
364:   MatViewFromOptions(M, NULL, "-M_view");
365:   KSPSetOperators(ksp, M, M);
366:   KSPSetFromOptions(ksp);
367:   KSPSolve(ksp, rhs, fhat);
368:   PetscObjectSetName((PetscObject)fhat, "fhat");
369:   VecViewFromOptions(fhat, NULL, "-fhat_view");

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

379:   KSPDestroy(&ksp);
380:   MatDestroy(&M);
381:   MatDestroy(&M_p);
382:   DMRestoreGlobalVector(dm, &fhat);
383:   DMRestoreGlobalVector(dm, &rhs);

385:   return 0;
386: }

388: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
389: {
390:   MPI_Comm  comm;
391:   KSP       ksp;
392:   Mat       M;            /* FEM mass matrix */
393:   Mat       M_p, PM_p;    /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
394:   Vec       f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
395:   PetscReal pmoments[3];  /* \int f, \int x f, \int r^2 f */
396:   PetscReal fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
397:   PetscInt  m;

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

402:   KSPCreate(comm, &ksp);
403:   KSPSetOptionsPrefix(ksp, "ftop_");
404:   KSPSetFromOptions(ksp);

406:   DMGetGlobalVector(dm, &fhat);
407:   DMGetGlobalVector(dm, &rhs);

409:   DMCreateMassMatrix(sw, dm, &M_p);
410:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

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

428:   KSPSetOperators(ksp, M_p, PM_p);
429:   KSPSolveTranspose(ksp, rhs, f);
430:   PetscObjectSetName((PetscObject)fhat, "fhat");
431:   VecViewFromOptions(fhat, NULL, "-fhat_view");

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

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

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

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

477:     DMGetField(dm, field, NULL, &obj);
478:     PetscObjectGetClassId(obj, &id);
479:     if (id == PETSCFE_CLASSID) {
480:       PetscFE fe = (PetscFE)obj;

482:       PetscFEGetQuadrature(fe, &quad);
483:       PetscFEGetNumComponents(fe, &Nc);
484:     } else if (id == PETSCFV_CLASSID) {
485:       PetscFV fv = (PetscFV)obj;

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

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

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

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

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

560: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
561: {
562:   MPI_Comm  comm;
563:   KSP       ksp;
564:   Mat       M;                  /* FEM mass matrix */
565:   Mat       M_p;                /* Particle mass matrix */
566:   Vec       f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
567:   PetscReal pmoments[3];        /* \int f, \int x f, \int r^2 f */
568:   PetscReal fmoments[3];        /* \int \hat f, \int x \hat f, \int r^2 \hat f */
569:   PetscInt  m;

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

579:   DMCreateMassMatrix(sw, dm, &M_p);
580:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

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

591:   DMCreateMatrix(dm, &M);
592:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);

594:   InterpolateGradient(dm, fhat, grad);

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

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

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

618:   return 0;
619: }

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

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

643: /*TEST

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

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

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

661:   test:
662:     suffix: proj_quad_0
663:     requires: triangle
664:     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
665:     filter: grep -v marker | grep -v atomic | grep -v usage

667:   test:
668:     suffix: proj_quad_2_faces
669:     requires: triangle
670:     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
671:     filter: grep -v marker | grep -v atomic | grep -v usage

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

679:   test:
680:     suffix: proj_quad_5P
681:     requires: triangle
682:     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
683:     filter: grep -v marker | grep -v atomic | grep -v usage

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

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

697:   test:
698:     suffix: proj_tri_3d
699:     requires: ctetgen
700:     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
701:     filter: grep -v marker | grep -v atomic | grep -v usage

703:   test:
704:     suffix: proj_tri_3d_2_faces
705:     requires: ctetgen
706:     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
707:     filter: grep -v marker | grep -v atomic | grep -v usage

709:   test:
710:     suffix: proj_tri_3d_5P
711:     requires: ctetgen
712:     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
713:     filter: grep -v marker | grep -v atomic | grep -v usage

715:   test:
716:     suffix: proj_tri_3d_mdx
717:     requires: ctetgen
718:     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
719:     filter: grep -v marker | grep -v atomic | grep -v usage

721:   test:
722:     suffix: proj_tri_3d_mdx_5P
723:     requires: ctetgen
724:     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
725:     filter: grep -v marker | grep -v atomic | grep -v usage

727:   test:
728:     suffix: proj_tri_3d_mdx_2_faces
729:     requires: ctetgen
730:     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
731:     filter: grep -v marker | grep -v atomic | grep -v usage

733:   test:
734:     suffix: proj_tri_3d_mdx_5P_2_faces
735:     requires: ctetgen
736:     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
737:     filter: grep -v marker | grep -v atomic | grep -v usage

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

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

757: TEST*/