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;
 80:       if (!flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown function %s",fstring);
 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;

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

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

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

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

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

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

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

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

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

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

201:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
202:   PetscRandomSetInterval(rnd, -1.0, 1.0);
203:   PetscRandomSetFromOptions(rnd);
204:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
205:   PetscRandomSetInterval(rndp, -interval, interval);
206:   PetscRandomSetFromOptions(rndp);

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

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

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

244:       for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value); coords[n*dim+d] += PetscRealPart(value);}
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;
265:   PetscErrorCode     ierr;

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

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

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

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

304: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
305:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
306:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
307:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
308: {
309:   f0[0] = x[0]*u[0];
310: }

312: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
313:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
314:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
315:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
316: {
317:   PetscInt d;

319:   f0[0] = 0.0;
320:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
321: }

323: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
324: {
325:   PetscDS        prob;
327:   PetscScalar    mom;

330:   DMGetDS(dm, &prob);
331:   PetscDSSetObjective(prob, 0, &f0_1);
332:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
333:   moments[0] = PetscRealPart(mom);
334:   PetscDSSetObjective(prob, 0, &f0_x);
335:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
336:   moments[1] = PetscRealPart(mom);
337:   PetscDSSetObjective(prob, 0, &f0_r2);
338:   DMPlexComputeIntegralFEM(dm, u, &mom, user);
339:   moments[2] = PetscRealPart(mom);
340:   return(0);
341: }

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

356:   PetscObjectGetComm((PetscObject) dm, &comm);
357:   KSPCreate(comm, &ksp);
358:   KSPSetOptionsPrefix(ksp, "ptof_");
359:   DMGetGlobalVector(dm, &fhat);
360:   DMGetGlobalVector(dm, &rhs);

362:   DMCreateMassMatrix(sw, dm, &M_p);
363:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

368:   /* create matrix RHS vector */
369:   MatMultTranspose(M_p, f, rhs);
370:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
371:   PetscObjectSetName((PetscObject) rhs,"rhs");
372:   VecViewFromOptions(rhs, NULL, "-rhs_view");

374:   DMCreateMatrix(dm, &M);
375:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
376:   MatViewFromOptions(M, NULL, "-M_view");
377:   KSPSetOperators(ksp, M, M);
378:   KSPSetFromOptions(ksp);
379:   KSPSolve(ksp, rhs, fhat);
380:   PetscObjectSetName((PetscObject) fhat,"fhat");
381:   VecViewFromOptions(fhat, NULL, "-fhat_view");

383:   /* Check moments of field */
384:   computeParticleMoments(sw, pmoments, user);
385:   computeFEMMoments(dm, fhat, fmoments, user);
386:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
387:   for (m = 0; m < 3; ++m) {
388:     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);
389:   }

391:   KSPDestroy(&ksp);
392:   MatDestroy(&M);
393:   MatDestroy(&M_p);
394:   DMRestoreGlobalVector(dm, &fhat);
395:   DMRestoreGlobalVector(dm, &rhs);

397:   return(0);
398: }

400: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
401: {

403:   MPI_Comm       comm;
404:   KSP            ksp;
405:   Mat            M;            /* FEM mass matrix */
406:   Mat            M_p, PM_p;    /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
407:   Vec            f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
408:   PetscReal      pmoments[3];  /* \int f, \int x f, \int r^2 f */
409:   PetscReal      fmoments[3];  /* \int \hat f, \int x \hat f, \int r^2 \hat f */
410:   PetscInt       m;

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

416:   KSPCreate(comm, &ksp);
417:   KSPSetOptionsPrefix(ksp, "ftop_");
418:   KSPSetFromOptions(ksp);

420:   DMGetGlobalVector(dm, &fhat);
421:   DMGetGlobalVector(dm, &rhs);

423:   DMCreateMassMatrix(sw, dm, &M_p);
424:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

429:   /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
430:   PetscObjectSetName((PetscObject) rhs,"rhs");
431:   VecViewFromOptions(rhs, NULL, "-rhs_view");
432:   DMCreateMatrix(dm, &M);
433:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
434:   MatViewFromOptions(M, NULL, "-M_view");
435:   MatMultTranspose(M, fhat, rhs);
436:   if (user->useBlockDiagPrec) {DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);}
437:   else                        {PetscObjectReference((PetscObject) M_p); PM_p = M_p;}

439:   KSPSetOperators(ksp, M_p, PM_p);
440:   KSPSolveTranspose(ksp, rhs, f);
441:   PetscObjectSetName((PetscObject) fhat,"fhat");
442:   VecViewFromOptions(fhat, NULL, "-fhat_view");

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

446:   /* Check moments */
447:   computeParticleMoments(sw, pmoments, user);
448:   computeFEMMoments(dm, fhat, fmoments, user);
449:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
450:   for (m = 0; m < 3; ++m) {
451:     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);
452:   }
453:   KSPDestroy(&ksp);
454:   MatDestroy(&M);
455:   MatDestroy(&M_p);
456:   MatDestroy(&PM_p);
457:   DMRestoreGlobalVector(dm, &fhat);
458:   DMRestoreGlobalVector(dm, &rhs);
459:   return(0);
460: }

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

478:   VecGetDM(locC, &dmC);
479:   VecSet(locC, 0.0);
480:   DMGetDimension(dm, &dim);
481:   DMGetCoordinateDim(dm, &coordDim);
482:   fegeom.dimEmbed = coordDim;
483:   DMGetLocalSection(dm, &section);
484:   PetscSectionGetNumFields(section, &numFields);
485:   for (field = 0; field < numFields; ++field) {
486:     PetscObject  obj;
487:     PetscClassId id;
488:     PetscInt     Nc;

490:     DMGetField(dm, field, NULL, &obj);
491:     PetscObjectGetClassId(obj, &id);
492:     if (id == PETSCFE_CLASSID) {
493:       PetscFE fe = (PetscFE) obj;

495:       PetscFEGetQuadrature(fe, &quad);
496:       PetscFEGetNumComponents(fe, &Nc);
497:     } else if (id == PETSCFV_CLASSID) {
498:       PetscFV fv = (PetscFV) obj;

500:       PetscFVGetQuadrature(fv, &quad);
501:       PetscFVGetNumComponents(fv, &Nc);
502:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
503:     numComponents += Nc;
504:   }
505:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
506:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
507:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
508:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
509:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
510:   for (v = vStart; v < vEnd; ++v) {
511:     PetscInt   *star = NULL;
512:     PetscInt    starSize, st, d, fc;

514:     PetscArrayzero(gradsum, coordDim*numComponents);
515:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
516:     for (st = 0; st < starSize*2; st += 2) {
517:       const PetscInt cell = star[st];
518:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
519:       PetscScalar   *x    = NULL;

521:       if ((cell < cStart) || (cell >= cEnd)) continue;
522:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
523:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
524:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
525:         PetscObject  obj;
526:         PetscClassId id;
527:         PetscInt     Nb, Nc, q, qc = 0;

529:         PetscArrayzero(grad, coordDim*numComponents);
530:         DMGetField(dm, field, NULL, &obj);
531:         PetscObjectGetClassId(obj, &id);
532:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
533:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
534:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
535:         for (q = 0; q < Nq; ++q) {
536:           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);
537:           if (ierr) {
538:             PetscErrorCode ierr2;
539:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
540:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
541:             ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
542:             
543:           }
544:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, 1, &x[fieldOffset], &fegeom, q, interpolant);}
545:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
546:           for (fc = 0; fc < Nc; ++fc) {
547:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

549:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
550:           }
551:         }
552:         fieldOffset += Nb;
553:       }
554:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
555:       for (fc = 0; fc < numComponents; ++fc) {
556:         for (d = 0; d < coordDim; ++d) {
557:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
558:         }
559:       }
560:       if (debug) {
561:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
562:         for (fc = 0; fc < numComponents; ++fc) {
563:           for (d = 0; d < coordDim; ++d) {
564:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
565:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
566:           }
567:         }
568:         PetscPrintf(PETSC_COMM_SELF, "]\n");
569:       }
570:     }
571:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
572:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
573:   }
574:   PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
575:   return(0);
576: }

578: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
579: {

581:   MPI_Comm       comm;
582:   KSP            ksp;
583:   Mat            M;                   /* FEM mass matrix */
584:   Mat            M_p;                 /* Particle mass matrix */
585:   Vec            f, rhs, fhat, grad;  /* Particle field f, \int phi_i f, FEM field */
586:   PetscReal      pmoments[3];         /* \int f, \int x f, \int r^2 f */
587:   PetscReal      fmoments[3];         /* \int \hat f, \int x \hat f, \int r^2 \hat f */
588:   PetscInt       m;

592:   PetscObjectGetComm((PetscObject) dm, &comm);
593:   KSPCreate(comm, &ksp);
594:   KSPSetOptionsPrefix(ksp, "ptof_");
595:   DMGetGlobalVector(dm, &fhat);
596:   DMGetGlobalVector(dm, &rhs);
597:   DMGetGlobalVector(dm, &grad);

599:   DMCreateMassMatrix(sw, dm, &M_p);
600:   MatViewFromOptions(M_p, NULL, "-M_p_view");

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

605:   /* create matrix RHS vector */
606:   MatMultTranspose(M_p, f, rhs);
607:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
608:   PetscObjectSetName((PetscObject) rhs,"rhs");
609:   VecViewFromOptions(rhs, NULL, "-rhs_view");

611:   DMCreateMatrix(dm, &M);
612:   DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);

614:   InterpolateGradient(dm, fhat, grad);

616:   MatViewFromOptions(M, NULL, "-M_view");
617:   KSPSetOperators(ksp, M, M);
618:   KSPSetFromOptions(ksp);
619:   KSPSolve(ksp, rhs, grad);
620:   PetscObjectSetName((PetscObject) fhat,"fhat");
621:   VecViewFromOptions(fhat, NULL, "-fhat_view");

623:   /* Check moments of field */
624:   computeParticleMoments(sw, pmoments, user);
625:   computeFEMMoments(dm, grad, fmoments, user);
626:   PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
627:   for (m = 0; m < 3; ++m) {
628:     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);
629:   }

631:   KSPDestroy(&ksp);
632:   MatDestroy(&M);
633:   MatDestroy(&M_p);
634:   DMRestoreGlobalVector(dm, &fhat);
635:   DMRestoreGlobalVector(dm, &rhs);
636:   DMRestoreGlobalVector(dm, &grad);

638:   return(0);
639: }

641: int main (int argc, char * argv[]) {
642:   MPI_Comm       comm;
643:   DM             dm, sw;
644:   AppCtx         user;

647:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
648:   comm = PETSC_COMM_WORLD;
649:   ProcessOptions(comm, &user);
650:   CreateMesh(comm, &dm, &user);
651:   CreateFEM(dm, &user);
652:   CreateParticles(dm, &sw, &user);
653:   TestL2ProjectionParticlesToField(dm, sw, &user);
654:   TestL2ProjectionFieldToParticles(dm, sw, &user);
655:   TestFieldGradientProjection(dm, sw, &user);
656:   DMDestroy(&dm);
657:   DMDestroy(&sw);
658:   PetscFinalize();
659:   return ierr;
660: }

662: /*TEST

664:   # Swarm does not handle complex or quad
665:   build:
666:     requires: !complex double

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

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

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

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

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

698:   test:
699:     suffix: proj_quad_5P
700:     requires: triangle
701:     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
702:     filter: grep -v marker | grep -v atomic | grep -v usage

704:   test:
705:     suffix: proj_tri_mdx
706:     requires: triangle
707:     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
708:     filter: grep -v marker | grep -v atomic | grep -v usage

710:   test:
711:     suffix: proj_tri_mdx_5P
712:     requires: triangle
713:     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
714:     filter: grep -v marker | grep -v atomic | grep -v usage

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

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

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

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

740:   test:
741:     suffix: proj_tri_3d_mdx_5P
742:     requires: ctetgen
743:     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
744:     filter: grep -v marker | grep -v atomic | grep -v usage

746:   test:
747:     suffix: proj_tri_3d_mdx_2_faces
748:     requires: ctetgen
749:     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
750:     filter: grep -v marker | grep -v atomic | grep -v usage

752:   test:
753:     suffix: proj_tri_3d_mdx_5P_2_faces
754:     requires: ctetgen
755:     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
756:     filter: grep -v marker | grep -v atomic | grep -v usage

758:   test:
759:     suffix: proj_quad_lsqr_scale
760:     requires: !complex
761:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
762:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
763:       -particlesPerCell 200 \
764:       -ptof_pc_type lu  \
765:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none

767:   test:
768:     suffix: proj_quad_lsqr_prec_scale
769:     requires: !complex
770:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
771:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
772:       -particlesPerCell 200 \
773:       -ptof_pc_type lu  \
774:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type lu -ftop_pc_factor_shift_type nonzero -block_diag_prec

776: TEST*/