Actual source code: ex2.c

  1: static char help[] = "Tests L2 projection with DMSwarm using delta function particles and deposition of linear shape.\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 shape;
 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->shape            = PETSC_FALSE;

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

156:   PetscFunctionBeginUser;
157:   PetscCall(DMGetDimension(dm, &dim));
158:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
159:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
160:   PetscCall(PetscFECreateByCell(PetscObjectComm((PetscObject)dm), dim, 1, ct, NULL, -1, &fe));
161:   PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
162:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
163:   PetscCall(DMCreateDS(dm));
164:   PetscCall(PetscFEDestroy(&fe));
165:   /* Setup to form mass matrix */
166:   PetscCall(DMGetDS(dm, &ds));
167:   PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
168:   PetscFunctionReturn(PETSC_SUCCESS);
169: }

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

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

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

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

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

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

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

256: static PetscErrorCode CreateParticles_Shape(DM dm, DM *sw, AppCtx *user)
257: {
258:   PetscDS          prob;
259:   PetscFE          fe;
260:   PetscQuadrature  quad;
261:   PetscScalar     *vals;
262:   PetscReal       *v0, *J, *invJ, detJ, *coords, *xi0;
263:   PetscInt        *cellid;
264:   const PetscReal *qpoints, *qweights;
265:   PetscInt         cStart, cEnd, c, Nq, q, dim;

267:   PetscFunctionBeginUser;
268:   PetscCall(DMGetDimension(dm, &dim));
269:   PetscCall(DMGetDS(dm, &prob));
270:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
271:   PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
272:   PetscCall(PetscFEGetQuadrature(fe, &quad));
273:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
274:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
275:   PetscCall(DMSetType(*sw, DMSWARM));
276:   PetscCall(DMSetDimension(*sw, dim));
277:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
278:   PetscCall(DMSwarmSetCellDM(*sw, dm));
279:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_REAL));
280:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
281:   PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
282:   PetscCall(DMSetFromOptions(*sw));
283:   PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
284:   for (c = 0; c < dim; c++) xi0[c] = -1.;
285:   PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
286:   PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
287:   PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
288:   for (c = cStart; c < cEnd; ++c) {
289:     for (q = 0; q < Nq; ++q) {
290:       PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
291:       cellid[c * Nq + q] = c;
292:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
293:       PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
294:       vals[c * Nq + q] *= qweights[q] * detJ;
295:     }
296:   }
297:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
298:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
299:   PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
300:   PetscCall(PetscFree4(xi0, v0, J, invJ));
301:   PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
302:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
303:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
308: {
309:   DM                 dm;
310:   const PetscReal   *coords;
311:   const PetscScalar *w;
312:   PetscReal          mom[3] = {0.0, 0.0, 0.0};
313:   PetscInt           cell, cStart, cEnd, dim;

315:   PetscFunctionBeginUser;
316:   PetscCall(DMGetDimension(sw, &dim));
317:   PetscCall(DMSwarmGetCellDM(sw, &dm));
318:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
319:   PetscCall(DMSwarmSortGetAccess(sw));
320:   PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
321:   PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
322:   for (cell = cStart; cell < cEnd; ++cell) {
323:     PetscInt *pidx;
324:     PetscInt  Np, p, d;

326:     PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
327:     for (p = 0; p < Np; ++p) {
328:       const PetscInt   idx = pidx[p];
329:       const PetscReal *c   = &coords[idx * dim];

331:       mom[0] += PetscRealPart(w[idx]);
332:       mom[1] += PetscRealPart(w[idx]) * c[0];
333:       for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
334:     }
335:     PetscCall(PetscFree(pidx));
336:   }
337:   PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
338:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
339:   PetscCall(DMSwarmSortRestoreAccess(sw));
340:   PetscCall(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: 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[])
345: {
346:   f0[0] = u[0];
347: }

349: 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[])
350: {
351:   f0[0] = x[0] * u[0];
352: }

354: 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[])
355: {
356:   PetscInt d;

358:   f0[0] = 0.0;
359:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
360: }

362: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
363: {
364:   PetscDS     prob;
365:   PetscScalar mom;

367:   PetscFunctionBeginUser;
368:   PetscCall(DMGetDS(dm, &prob));
369:   PetscCall(PetscDSSetObjective(prob, 0, &f0_1));
370:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
371:   moments[0] = PetscRealPart(mom);
372:   PetscCall(PetscDSSetObjective(prob, 0, &f0_x));
373:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
374:   moments[1] = PetscRealPart(mom);
375:   PetscCall(PetscDSSetObjective(prob, 0, &f0_r2));
376:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
377:   moments[2] = PetscRealPart(mom);
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, Vec fhat, AppCtx *user)
382: {
383:   const char *fieldnames[1] = {"w_q"};
384:   Vec         fields[1]     = {fhat};
385:   PetscReal   pmoments[3]; // \int f, \int x f, \int r^2 f
386:   PetscReal   fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f

388:   PetscFunctionBeginUser;
389:   PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_FORWARD));

391:   /* Check moments of field */
392:   PetscCall(computeParticleMoments(sw, pmoments, user));
393:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
394:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
395:   for (PetscInt m = 0; m < 3; ++m) {
396:     PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
397:                user->momentTol);
398:   }
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, Vec fhat, AppCtx *user)
403: {
404:   const char *fieldnames[1] = {"w_q"};
405:   Vec         fields[1]     = {fhat};
406:   PetscReal   pmoments[3]; // \int f, \int x f, \int r^2 f
407:   PetscReal   fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f

409:   PetscFunctionBeginUser;
410:   PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_REVERSE));

412:   /* Check moments */
413:   PetscCall(computeParticleMoments(sw, pmoments, user));
414:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
415:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
416:   for (PetscInt m = 0; m < 3; ++m) {
417:     PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
418:                user->momentTol);
419:   }
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
424: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
425: {
426:   DM_Plex         *mesh  = (DM_Plex *)dm->data;
427:   PetscInt         debug = mesh->printFEM;
428:   DM               dmC;
429:   PetscSection     section;
430:   PetscQuadrature  quad = NULL;
431:   PetscScalar     *interpolant, *gradsum;
432:   PetscFEGeom      fegeom;
433:   PetscReal       *coords;
434:   const PetscReal *quadPoints, *quadWeights;
435:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;

437:   PetscFunctionBegin;
438:   PetscCall(VecGetDM(locC, &dmC));
439:   PetscCall(VecSet(locC, 0.0));
440:   PetscCall(DMGetDimension(dm, &dim));
441:   PetscCall(DMGetCoordinateDim(dm, &coordDim));
442:   fegeom.dimEmbed = coordDim;
443:   PetscCall(DMGetLocalSection(dm, &section));
444:   PetscCall(PetscSectionGetNumFields(section, &numFields));
445:   for (field = 0; field < numFields; ++field) {
446:     PetscObject  obj;
447:     PetscClassId id;
448:     PetscInt     Nc;

450:     PetscCall(DMGetField(dm, field, NULL, &obj));
451:     PetscCall(PetscObjectGetClassId(obj, &id));
452:     if (id == PETSCFE_CLASSID) {
453:       PetscFE fe = (PetscFE)obj;

455:       PetscCall(PetscFEGetQuadrature(fe, &quad));
456:       PetscCall(PetscFEGetNumComponents(fe, &Nc));
457:     } else if (id == PETSCFV_CLASSID) {
458:       PetscFV fv = (PetscFV)obj;

460:       PetscCall(PetscFVGetQuadrature(fv, &quad));
461:       PetscCall(PetscFVGetNumComponents(fv, &Nc));
462:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
463:     numComponents += Nc;
464:   }
465:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
466:   PetscCheck(!(qNc != 1) || !(qNc != numComponents), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Quadrature components %" PetscInt_FMT " != %" PetscInt_FMT " field components", qNc, numComponents);
467:   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));
468:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
469:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
470:   for (v = vStart; v < vEnd; ++v) {
471:     PetscInt *star = NULL;
472:     PetscInt  starSize, st, d, fc;

474:     PetscCall(PetscArrayzero(gradsum, coordDim * numComponents));
475:     PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
476:     for (st = 0; st < starSize * 2; st += 2) {
477:       const PetscInt cell = star[st];
478:       PetscScalar   *grad = &gradsum[coordDim * numComponents];
479:       PetscScalar   *x    = NULL;

481:       if ((cell < cStart) || (cell >= cEnd)) continue;
482:       PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
483:       PetscCall(DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x));
484:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
485:         PetscObject  obj;
486:         PetscClassId id;
487:         PetscInt     Nb, Nc, q, qc = 0;

489:         PetscCall(PetscArrayzero(grad, coordDim * numComponents));
490:         PetscCall(DMGetField(dm, field, NULL, &obj));
491:         PetscCall(PetscObjectGetClassId(obj, &id));
492:         if (id == PETSCFE_CLASSID) {
493:           PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
494:           PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
495:         } else if (id == PETSCFV_CLASSID) {
496:           PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
497:           Nb = 1;
498:         } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
499:         for (q = 0; q < Nq; ++q) {
500:           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);
501:           if (id == PETSCFE_CLASSID) PetscCall(PetscFEInterpolateGradient_Static((PetscFE)obj, 1, &x[fieldOffset], &fegeom, q, interpolant));
502:           else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
503:           for (fc = 0; fc < Nc; ++fc) {
504:             const PetscReal wt = quadWeights[q * qNc + qc + fc];

506:             for (d = 0; d < coordDim; ++d) grad[fc * coordDim + d] += interpolant[fc * dim + d] * wt * fegeom.detJ[q];
507:           }
508:         }
509:         fieldOffset += Nb;
510:       }
511:       PetscCall(DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x));
512:       for (fc = 0; fc < numComponents; ++fc) {
513:         for (d = 0; d < coordDim; ++d) gradsum[fc * coordDim + d] += grad[fc * coordDim + d];
514:       }
515:       if (debug) {
516:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " gradient: [", cell));
517:         for (fc = 0; fc < numComponents; ++fc) {
518:           for (d = 0; d < coordDim; ++d) {
519:             if (fc || d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
520:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc * coordDim + d])));
521:           }
522:         }
523:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
524:       }
525:     }
526:     PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
527:     PetscCall(DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES));
528:   }
529:   PetscCall(PetscFree6(gradsum, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
534: {
535:   MPI_Comm  comm;
536:   KSP       ksp;
537:   Mat       M;                  /* FEM mass matrix */
538:   Mat       M_p;                /* Particle mass matrix */
539:   Vec       f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
540:   PetscReal pmoments[3];        /* \int f, \int x f, \int r^2 f */
541:   PetscReal fmoments[3];        /* \int \hat f, \int x \hat f, \int r^2 \hat f */
542:   PetscInt  m;

544:   PetscFunctionBeginUser;
545:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
546:   PetscCall(KSPCreate(comm, &ksp));
547:   PetscCall(KSPSetOptionsPrefix(ksp, "ptof_"));
548:   PetscCall(DMGetGlobalVector(dm, &fhat));
549:   PetscCall(DMGetGlobalVector(dm, &rhs));
550:   PetscCall(DMGetGlobalVector(dm, &grad));

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

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

558:   /* create matrix RHS vector */
559:   PetscCall(MatMultTranspose(M_p, f, rhs));
560:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
561:   PetscCall(PetscObjectSetName((PetscObject)rhs, "rhs"));
562:   PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));

564:   PetscCall(DMCreateMatrix(dm, &M));
565:   PetscCall(DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user));

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

569:   PetscCall(MatViewFromOptions(M, NULL, "-M_view"));
570:   PetscCall(KSPSetOperators(ksp, M, M));
571:   PetscCall(KSPSetFromOptions(ksp));
572:   PetscCall(KSPSolve(ksp, rhs, grad));
573:   PetscCall(PetscObjectSetName((PetscObject)fhat, "fhat"));
574:   PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view"));

576:   /* Check moments of field */
577:   PetscCall(computeParticleMoments(sw, pmoments, user));
578:   PetscCall(computeFEMMoments(dm, grad, fmoments, user));
579:   PetscCall(PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
580:   for (m = 0; m < 3; ++m) {
581:     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);
582:   }

584:   PetscCall(KSPDestroy(&ksp));
585:   PetscCall(MatDestroy(&M));
586:   PetscCall(MatDestroy(&M_p));
587:   PetscCall(DMRestoreGlobalVector(dm, &fhat));
588:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
589:   PetscCall(DMRestoreGlobalVector(dm, &grad));

591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }

594: static PetscErrorCode TestL2Projection_Shape(DM dm, DM sw, AppCtx *user)
595: {
596:   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
597:   KSP       ksp;
598:   Mat       mass;
599:   Vec       u, rhs, uproj;
600:   void     *ctxs[1];
601:   PetscReal error;

603:   PetscFunctionBeginUser;
604:   ctxs[0]  = (void *)user;
605:   funcs[0] = linear;

607:   PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &u));
608:   PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
609:   PetscCall(DMGetGlobalVector(dm, &rhs));
610:   PetscCall(DMCreateMassMatrix(sw, dm, &mass));
611:   PetscCall(MatMultTranspose(mass, u, rhs));
612:   PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
613:   PetscCall(MatDestroy(&mass));
614:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &u));
615:   PetscCall(DMGetGlobalVector(dm, &uproj));
616:   PetscCall(DMCreateMatrix(dm, &mass));
617:   PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
618:   PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
619:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
620:   PetscCall(KSPSetOperators(ksp, mass, mass));
621:   PetscCall(KSPSetFromOptions(ksp));
622:   PetscCall(KSPSolve(ksp, rhs, uproj));
623:   PetscCall(KSPDestroy(&ksp));
624:   PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
625:   PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
626:   PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
627:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
628:   PetscCall(MatDestroy(&mass));
629:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
630:   PetscCall(DMRestoreGlobalVector(dm, &uproj));
631:   PetscFunctionReturn(PETSC_SUCCESS);
632: }

634: int main(int argc, char *argv[])
635: {
636:   MPI_Comm comm;
637:   DM       dm, sw;
638:   AppCtx   user;

640:   PetscFunctionBeginUser;
641:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
642:   comm = PETSC_COMM_WORLD;
643:   PetscCall(ProcessOptions(comm, &user));
644:   PetscCall(CreateMesh(comm, &dm, &user));
645:   PetscCall(CreateFEM(dm, &user));
646:   if (!user.shape) {
647:     Vec fhat;

649:     PetscCall(CreateParticles(dm, &sw, &user));
650:     PetscCall(DMGetGlobalVector(dm, &fhat));
651:     PetscCall(TestL2ProjectionParticlesToField(dm, sw, fhat, &user));
652:     PetscCall(TestL2ProjectionFieldToParticles(dm, sw, fhat, &user));
653:     PetscCall(TestFieldGradientProjection(dm, sw, &user));
654:     PetscCall(DMRestoreGlobalVector(dm, &fhat));
655:   } else {
656:     PetscCall(CreateParticles_Shape(dm, &sw, &user));
657:     PetscCall(TestL2Projection_Shape(dm, sw, &user));
658:   }
659:   PetscCall(DMDestroy(&dm));
660:   PetscCall(DMDestroy(&sw));
661:   PetscCall(PetscFinalize());
662:   return 0;
663: }

665: /*TEST

667:   # Swarm does not handle complex or quad
668:   build:
669:     requires: !complex double

671:   test:
672:     suffix: proj_tri_0
673:     requires: triangle
674:     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
675:     filter: grep -v marker | grep -v atomic | grep -v usage

677:   test:
678:     suffix: proj_tri_2_faces
679:     requires: triangle
680:     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
681:     filter: grep -v marker | grep -v atomic | grep -v usage

683:   test:
684:     suffix: proj_quad_0
685:     requires: triangle
686:     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
687:     filter: grep -v marker | grep -v atomic | grep -v usage

689:   test:
690:     suffix: proj_quad_2_faces
691:     requires: triangle
692:     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
693:     filter: grep -v marker | grep -v atomic | grep -v usage

695:   test:
696:     suffix: proj_tri_5P
697:     requires: triangle
698:     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
699:     filter: grep -v marker | grep -v atomic | grep -v usage

701:   test:
702:     suffix: proj_quad_5P
703:     requires: triangle
704:     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
705:     filter: grep -v marker | grep -v atomic | grep -v usage

707:   test:
708:     suffix: proj_tri_mdx
709:     requires: triangle
710:     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
711:     filter: grep -v marker | grep -v atomic | grep -v usage

713:   test:
714:     suffix: proj_tri_mdx_5P
715:     requires: triangle
716:     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
717:     filter: grep -v marker | grep -v atomic | grep -v usage

719:   test:
720:     suffix: proj_tri_3d
721:     requires: ctetgen
722:     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
723:     filter: grep -v marker | grep -v atomic | grep -v usage

725:   test:
726:     suffix: proj_tri_3d_2_faces
727:     requires: ctetgen
728:     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
729:     filter: grep -v marker | grep -v atomic | grep -v usage

731:   test:
732:     suffix: proj_tri_3d_5P
733:     requires: ctetgen
734:     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
735:     filter: grep -v marker | grep -v atomic | grep -v usage

737:   test:
738:     suffix: proj_tri_3d_mdx
739:     requires: ctetgen
740:     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
741:     filter: grep -v marker | grep -v atomic | grep -v usage

743:   test:
744:     suffix: proj_tri_3d_mdx_5P
745:     requires: ctetgen
746:     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
747:     filter: grep -v marker | grep -v atomic | grep -v usage

749:   test:
750:     suffix: proj_tri_3d_mdx_2_faces
751:     requires: ctetgen
752:     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
753:     filter: grep -v marker | grep -v atomic | grep -v usage

755:   test:
756:     suffix: proj_tri_3d_mdx_5P_2_faces
757:     requires: ctetgen
758:     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
759:     filter: grep -v marker | grep -v atomic | grep -v usage

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

770:   test:
771:     suffix: proj_quad_lsqr_prec_scale
772:     requires: !complex
773:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
774:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
775:       -particlesPerCell 200 \
776:       -ptof_pc_type lu  \
777:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero
778:   test:
779:     suffix: proj_shape_linear_tri_2d
780:     requires: ctetgen triangle
781:     args: -shape -dm_plex_box_faces 2,2\
782:           -dm_view -sw_view\
783:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
784:           -pc_type lu
785:   test:
786:     suffix: proj_shape_linear_quad_2d
787:     requires: ctetgen triangle
788:     args: -shape -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
789:           -dm_view -sw_view\
790:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
791:           -pc_type lu
792:   test:
793:     suffix: proj_shape_linear_tri_3d
794:     requires: ctetgen triangle
795:     args: -shape -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
796:           -dm_view -sw_view\
797:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
798:           -pc_type lu
799:   test:
800:     suffix: proj_shape_linear_quad_3d
801:     requires: ctetgen triangle
802:     args: -shape -dm_plex_dim 3 -dm_plex_simplex 0\
803:           -dm_plex_box_faces 2,2,2\
804:           -dm_view -sw_view\
805:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
806:           -pc_type lu
807: TEST*/