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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 48: }

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

 55:   PetscFunctionBeginUser;
 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:   PetscCall(PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL));
 65:   PetscCall(PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL));
 66:   PetscCall(PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL));
 67:   PetscCall(PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL));
 68:   PetscCall(PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL));
 69:   PetscCall(PetscStrcmp(fstring, "linear", &flag));
 70:   if (flag) {
 71:     options->func = linear;
 72:   } else {
 73:     PetscCall(PetscStrcmp(fstring, "sin", &flag));
 74:     if (flag) {
 75:       options->func = sinx;
 76:     } else {
 77:       PetscCall(PetscStrcmp(fstring, "x2_x4", &flag));
 78:       options->func = x2_x4;
 79:       PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown function %s", fstring);
 80:     }
 81:   }
 82:   PetscCall(PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL));
 83:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 98:   PetscFunctionBeginUser;
 99:   PetscCall(DMGetBoundingBox(dm, low, high));
100:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
101:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
102:   PetscCall(PetscRandomSetInterval(rnd, -interval, interval));
103:   PetscCall(PetscRandomSetFromOptions(rnd));
104:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
105:   PetscCall(DMGetCoordinateDim(dm, &cdim));
106:   PetscCall(PetscCalloc1(cdim, &hh));
107:   for (d = 0; d < cdim; ++d) hh[d] = (user->L[d]) / PetscPowReal(cEnd, 1. / cdim);
108:   PetscCall(VecGetLocalSize(coordinates, &N));
109:   PetscCall(VecGetBlockSize(coordinates, &bs));
110:   PetscCheck(bs == cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, cdim);
111:   PetscCall(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:       PetscCall(PetscRandomGetValue(rnd, &value));
117:       coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value * hh[d])));
118:     }
119:   }
120:   PetscCall(VecRestoreArray(coordinates, &coords));
121:   PetscCall(PetscRandomDestroy(&rnd));
122:   PetscCall(PetscFree(hh));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

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

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

137:   PetscCall(DMGetCoordinateDim(*dm, &cdim));
138:   PetscCall(DMGetBoundingBox(*dm, low, high));
139:   for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
140:   PetscCall(PerturbVertices(*dm, user));
141:   PetscFunctionReturn(PETSC_SUCCESS);
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;

157:   PetscFunctionBeginUser;
158:   PetscCall(DMGetDimension(dm, &dim));
159:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
160:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
161:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
162:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe));
163:   PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
164:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
165:   PetscCall(DMCreateDS(dm));
166:   PetscCall(PetscFEDestroy(&fe));
167:   /* Setup to form mass matrix */
168:   PetscCall(DMGetDS(dm, &ds));
169:   PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
170:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

211:   PetscCall(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:       PetscCall(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:       PetscCall(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:           PetscCall(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:   PetscCall(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:         PetscCall(PetscRandomGetValue(rndp, &value));
243:         coords[n * dim + d] += PetscRealPart(value);
244:       }
245:       PetscCall(user->func(dim, 0.0, &coords[n * dim], 1, &vals[c], user));
246:     }
247:   }
248:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
249:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
250:   PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
251:   PetscCall(PetscRandomDestroy(&rnd));
252:   PetscCall(PetscRandomDestroy(&rndp));
253:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
254:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
255:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

277:     PetscCall(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:     PetscCall(PetscFree(pidx));
287:   }
288:   PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
289:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
290:   PetscCall(DMSwarmSortRestoreAccess(sw));
291:   PetscCall(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
292:   PetscFunctionReturn(PETSC_SUCCESS);
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;

318:   PetscFunctionBeginUser;
319:   PetscCall(DMGetDS(dm, &prob));
320:   PetscCall(PetscDSSetObjective(prob, 0, &f0_1));
321:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
322:   moments[0] = PetscRealPart(mom);
323:   PetscCall(PetscDSSetObjective(prob, 0, &f0_x));
324:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
325:   moments[1] = PetscRealPart(mom);
326:   PetscCall(PetscDSSetObjective(prob, 0, &f0_r2));
327:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
328:   moments[2] = PetscRealPart(mom);
329:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

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

371:   /* Check moments of field */
372:   PetscCall(computeParticleMoments(sw, pmoments, user));
373:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
374:   PetscCall(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) {
376:     PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, comm, PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), user->momentTol);
377:   }

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

385:   PetscFunctionReturn(PETSC_SUCCESS);
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;

399:   PetscFunctionBeginUser;
400:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));

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

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

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

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

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

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

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

435:   /* Check moments */
436:   PetscCall(computeParticleMoments(sw, pmoments, user));
437:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
438:   PetscCall(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) {
440:     PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, comm, PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), user->momentTol);
441:   }
442:   PetscCall(KSPDestroy(&ksp));
443:   PetscCall(MatDestroy(&M));
444:   PetscCall(MatDestroy(&M_p));
445:   PetscCall(MatDestroy(&PM_p));
446:   PetscCall(DMRestoreGlobalVector(dm, &fhat));
447:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
448:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
466:   PetscCall(VecGetDM(locC, &dmC));
467:   PetscCall(VecSet(locC, 0.0));
468:   PetscCall(DMGetDimension(dm, &dim));
469:   PetscCall(DMGetCoordinateDim(dm, &coordDim));
470:   fegeom.dimEmbed = coordDim;
471:   PetscCall(DMGetLocalSection(dm, &section));
472:   PetscCall(PetscSectionGetNumFields(section, &numFields));
473:   for (field = 0; field < numFields; ++field) {
474:     PetscObject  obj;
475:     PetscClassId id;
476:     PetscInt     Nc;

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

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

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

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

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

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

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

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

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

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

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

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

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

595:   PetscCall(InterpolateGradient(dm, fhat, grad));

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

604:   /* Check moments of field */
605:   PetscCall(computeParticleMoments(sw, pmoments, user));
606:   PetscCall(computeFEMMoments(dm, grad, fmoments, user));
607:   PetscCall(PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
608:   for (m = 0; m < 3; ++m) {
609:     PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, comm, PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), user->momentTol);
610:   }

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

619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }

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

628:   PetscFunctionBeginUser;
629:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
630:   comm = PETSC_COMM_WORLD;
631:   PetscCall(ProcessOptions(comm, &user));
632:   PetscCall(CreateMesh(comm, &dm, &user));
633:   PetscCall(CreateFEM(dm, &user));
634:   PetscCall(CreateParticles(dm, &sw, &user));
635:   PetscCall(TestL2ProjectionParticlesToField(dm, sw, &user));
636:   PetscCall(TestL2ProjectionFieldToParticles(dm, sw, &user));
637:   PetscCall(TestFieldGradientProjection(dm, sw, &user));
638:   PetscCall(DMDestroy(&dm));
639:   PetscCall(DMDestroy(&sw));
640:   PetscCall(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*/