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();
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

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

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

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

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

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

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

143: 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[])
144: {
145:   g0[0] = 1.0;
146: }

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

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

170: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
171: {
172:   DMSwarmCellDM  celldm;
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      *swarm_cellid;
180:   PetscInt       Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d, Nfc;
181:   const char   **coordFields, *cellid;

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

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

199:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
200:   PetscCall(DMSwarmSetCellDM(*sw, dm));
201:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
202:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
203:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Ncell));
204:   PetscCall(DMSwarmSetLocalSizes(*sw, Ncell * Np, 0));
205:   PetscCall(DMSetFromOptions(*sw));

207:   PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
208:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
209:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
210:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);

212:   PetscCall(DMSwarmGetField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
213:   PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
214:   PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));

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

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

246:       for (d = 0; d < dim; ++d) {
247:         PetscCall(PetscRandomGetValue(rndp, &value));
248:         coords[n * dim + d] += PetscRealPart(value);
249:       }
250:       PetscCall(user->func(dim, 0.0, &coords[n * dim], 1, &vals[c], user));
251:     }
252:   }
253:   PetscCall(DMSwarmRestoreField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
254:   PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
255:   PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
256:   PetscCall(PetscRandomDestroy(&rnd));
257:   PetscCall(PetscRandomDestroy(&rndp));
258:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
259:   PetscCall(DMSwarmVectorDefineField(*sw, "w_q"));
260:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: static PetscErrorCode CreateParticles_Shape(DM dm, DM *sw, AppCtx *user)
265: {
266:   DMSwarmCellDM    celldm;
267:   PetscDS          prob;
268:   PetscFE          fe;
269:   PetscQuadrature  quad;
270:   PetscScalar     *vals;
271:   PetscReal       *v0, *J, *invJ, detJ, *coords, *xi0;
272:   PetscInt        *swarm_cellid;
273:   const PetscReal *qpoints, *qweights;
274:   PetscInt         cStart, cEnd, c, Nq, q, dim, Nfc;
275:   const char     **coordFields, *cellid;

277:   PetscFunctionBeginUser;
278:   PetscCall(DMGetDimension(dm, &dim));
279:   PetscCall(DMGetDS(dm, &prob));
280:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
281:   PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
282:   PetscCall(PetscFEGetQuadrature(fe, &quad));
283:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
284:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
285:   PetscCall(DMSetType(*sw, DMSWARM));
286:   PetscCall(DMSetDimension(*sw, dim));
287:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
288:   PetscCall(DMSwarmSetCellDM(*sw, dm));
289:   PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
290:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
291:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
292:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
293:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_REAL));
294:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
295:   PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
296:   PetscCall(DMSetFromOptions(*sw));
297:   PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
298:   for (c = 0; c < dim; c++) xi0[c] = -1.;
299:   PetscCall(DMSwarmGetField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
300:   PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
301:   PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
302:   for (c = cStart; c < cEnd; ++c) {
303:     for (q = 0; q < Nq; ++q) {
304:       PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
305:       swarm_cellid[c * Nq + q] = c;
306:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
307:       PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
308:       vals[c * Nq + q] *= qweights[q] * detJ;
309:     }
310:   }
311:   PetscCall(DMSwarmRestoreField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
312:   PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
313:   PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
314:   PetscCall(PetscFree4(xi0, v0, J, invJ));
315:   PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
316:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
317:   PetscCall(DMSwarmVectorDefineField(*sw, "w_q"));
318:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
323: {
324:   DM                 dm;
325:   const PetscReal   *coords;
326:   const PetscScalar *w;
327:   PetscReal          mom[3] = {0.0, 0.0, 0.0};
328:   PetscInt           cell, cStart, cEnd, dim;

330:   PetscFunctionBeginUser;
331:   PetscCall(DMGetDimension(sw, &dim));
332:   PetscCall(DMSwarmGetCellDM(sw, &dm));
333:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
334:   PetscCall(DMSwarmSortGetAccess(sw));
335:   PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
336:   PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
337:   for (cell = cStart; cell < cEnd; ++cell) {
338:     PetscInt *pidx;
339:     PetscInt  Np, p, d;

341:     PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
342:     for (p = 0; p < Np; ++p) {
343:       const PetscInt   idx = pidx[p];
344:       const PetscReal *c   = &coords[idx * dim];

346:       mom[0] += PetscRealPart(w[idx]);
347:       mom[1] += PetscRealPart(w[idx]) * c[0];
348:       for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
349:     }
350:     PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Np, &pidx));
351:   }
352:   PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
353:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
354:   PetscCall(DMSwarmSortRestoreAccess(sw));
355:   PetscCallMPI(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: 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[])
360: {
361:   f0[0] = u[0];
362: }

364: 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[])
365: {
366:   f0[0] = x[0] * u[0];
367: }

369: 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[])
370: {
371:   PetscInt d;

373:   f0[0] = 0.0;
374:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
375: }

377: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
378: {
379:   PetscDS     prob;
380:   PetscScalar mom;

382:   PetscFunctionBeginUser;
383:   PetscCall(DMGetDS(dm, &prob));
384:   PetscCall(PetscDSSetObjective(prob, 0, &f0_1));
385:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
386:   moments[0] = PetscRealPart(mom);
387:   PetscCall(PetscDSSetObjective(prob, 0, &f0_x));
388:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
389:   moments[1] = PetscRealPart(mom);
390:   PetscCall(PetscDSSetObjective(prob, 0, &f0_r2));
391:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
392:   moments[2] = PetscRealPart(mom);
393:   PetscFunctionReturn(PETSC_SUCCESS);
394: }

396: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, Vec fhat, AppCtx *user)
397: {
398:   const char *fieldnames[1] = {"w_q"};
399:   Vec         fields[1]     = {fhat};
400:   PetscReal   pmoments[3]; // \int f, \int x f, \int r^2 f
401:   PetscReal   fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f

403:   PetscFunctionBeginUser;
404:   PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD));

406:   /* Check moments of field */
407:   PetscCall(computeParticleMoments(sw, pmoments, user));
408:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
409:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
410:   for (PetscInt m = 0; m < 3; ++m) {
411:     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]),
412:                user->momentTol);
413:   }
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, Vec fhat, AppCtx *user)
418: {
419:   const char *fieldnames[1] = {"w_q"};
420:   Vec         fields[1]     = {fhat};
421:   PetscReal   pmoments[3]; // \int f, \int x f, \int r^2 f
422:   PetscReal   fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f

424:   PetscFunctionBeginUser;
425:   PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE));

427:   /* Check moments */
428:   PetscCall(computeParticleMoments(sw, pmoments, user));
429:   PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
430:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
431:   for (PetscInt m = 0; m < 3; ++m) {
432:     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]),
433:                user->momentTol);
434:   }
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
439: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
440: {
441:   DM_Plex         *mesh  = (DM_Plex *)dm->data;
442:   PetscInt         debug = mesh->printFEM;
443:   DM               dmC;
444:   PetscSection     section;
445:   PetscQuadrature  quad = NULL;
446:   PetscScalar     *interpolant, *gradsum;
447:   PetscFEGeom      fegeom;
448:   PetscReal       *coords;
449:   const PetscReal *quadPoints, *quadWeights;
450:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;

452:   PetscFunctionBegin;
453:   PetscCall(VecGetDM(locC, &dmC));
454:   PetscCall(VecSet(locC, 0.0));
455:   PetscCall(DMGetDimension(dm, &dim));
456:   PetscCall(DMGetCoordinateDim(dm, &coordDim));
457:   fegeom.dimEmbed = coordDim;
458:   PetscCall(DMGetLocalSection(dm, &section));
459:   PetscCall(PetscSectionGetNumFields(section, &numFields));
460:   for (field = 0; field < numFields; ++field) {
461:     PetscObject  obj;
462:     PetscClassId id;
463:     PetscInt     Nc;

465:     PetscCall(DMGetField(dm, field, NULL, &obj));
466:     PetscCall(PetscObjectGetClassId(obj, &id));
467:     if (id == PETSCFE_CLASSID) {
468:       PetscFE fe = (PetscFE)obj;

470:       PetscCall(PetscFEGetQuadrature(fe, &quad));
471:       PetscCall(PetscFEGetNumComponents(fe, &Nc));
472:     } else if (id == PETSCFV_CLASSID) {
473:       PetscFV fv = (PetscFV)obj;

475:       PetscCall(PetscFVGetQuadrature(fv, &quad));
476:       PetscCall(PetscFVGetNumComponents(fv, &Nc));
477:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
478:     numComponents += Nc;
479:   }
480:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
481:   PetscCheck(!(qNc != 1) || !(qNc != numComponents), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Quadrature components %" PetscInt_FMT " != %" PetscInt_FMT " field components", qNc, numComponents);
482:   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));
483:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
484:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
485:   for (v = vStart; v < vEnd; ++v) {
486:     PetscInt *star = NULL;
487:     PetscInt  starSize, st, d, fc;

489:     PetscCall(PetscArrayzero(gradsum, coordDim * numComponents));
490:     PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
491:     for (st = 0; st < starSize * 2; st += 2) {
492:       const PetscInt cell = star[st];
493:       PetscScalar   *grad = &gradsum[coordDim * numComponents];
494:       PetscScalar   *x    = NULL;

496:       if ((cell < cStart) || (cell >= cEnd)) continue;
497:       PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
498:       PetscCall(DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x));
499:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
500:         PetscObject  obj;
501:         PetscClassId id;
502:         PetscInt     Nb, Nc, q, qc = 0;

504:         PetscCall(PetscArrayzero(grad, coordDim * numComponents));
505:         PetscCall(DMGetField(dm, field, NULL, &obj));
506:         PetscCall(PetscObjectGetClassId(obj, &id));
507:         if (id == PETSCFE_CLASSID) {
508:           PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
509:           PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
510:         } else if (id == PETSCFV_CLASSID) {
511:           PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
512:           Nb = 1;
513:         } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
514:         for (q = 0; q < Nq; ++q) {
515:           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);
516:           if (id == PETSCFE_CLASSID) PetscCall(PetscFEInterpolateGradient_Static((PetscFE)obj, 1, &x[fieldOffset], &fegeom, q, interpolant));
517:           else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
518:           for (fc = 0; fc < Nc; ++fc) {
519:             const PetscReal wt = quadWeights[q * qNc + qc + fc];

521:             for (d = 0; d < coordDim; ++d) grad[fc * coordDim + d] += interpolant[fc * dim + d] * wt * fegeom.detJ[q];
522:           }
523:         }
524:         fieldOffset += Nb;
525:       }
526:       PetscCall(DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x));
527:       for (fc = 0; fc < numComponents; ++fc) {
528:         for (d = 0; d < coordDim; ++d) gradsum[fc * coordDim + d] += grad[fc * coordDim + d];
529:       }
530:       if (debug) {
531:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " gradient: [", cell));
532:         for (fc = 0; fc < numComponents; ++fc) {
533:           for (d = 0; d < coordDim; ++d) {
534:             if (fc || d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
535:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc * coordDim + d])));
536:           }
537:         }
538:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
539:       }
540:     }
541:     PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
542:     PetscCall(DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES));
543:   }
544:   PetscCall(PetscFree6(gradsum, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
549: {
550:   MPI_Comm  comm;
551:   KSP       ksp;
552:   Mat       M;                  /* FEM mass matrix */
553:   Mat       M_p;                /* Particle mass matrix */
554:   Vec       f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
555:   PetscReal pmoments[3];        /* \int f, \int x f, \int r^2 f */
556:   PetscReal fmoments[3];        /* \int \hat f, \int x \hat f, \int r^2 \hat f */
557:   PetscInt  m;

559:   PetscFunctionBeginUser;
560:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
561:   PetscCall(KSPCreate(comm, &ksp));
562:   PetscCall(KSPSetOptionsPrefix(ksp, "ptof_"));
563:   PetscCall(DMGetGlobalVector(dm, &fhat));
564:   PetscCall(DMGetGlobalVector(dm, &rhs));
565:   PetscCall(DMGetGlobalVector(dm, &grad));

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

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

573:   /* create matrix RHS vector */
574:   PetscCall(MatMultTranspose(M_p, f, rhs));
575:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
576:   PetscCall(PetscObjectSetName((PetscObject)rhs, "rhs"));
577:   PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));

579:   PetscCall(DMCreateMatrix(dm, &M));
580:   PetscCall(DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user));

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

584:   PetscCall(MatViewFromOptions(M, NULL, "-M_view"));
585:   PetscCall(KSPSetOperators(ksp, M, M));
586:   PetscCall(KSPSetFromOptions(ksp));
587:   PetscCall(KSPSolve(ksp, rhs, grad));
588:   PetscCall(PetscObjectSetName((PetscObject)fhat, "fhat"));
589:   PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view"));

591:   /* Check moments of field */
592:   PetscCall(computeParticleMoments(sw, pmoments, user));
593:   PetscCall(computeFEMMoments(dm, grad, fmoments, user));
594:   PetscCall(PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
595:   for (m = 0; m < 3; ++m) {
596:     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);
597:   }

599:   PetscCall(KSPDestroy(&ksp));
600:   PetscCall(MatDestroy(&M));
601:   PetscCall(MatDestroy(&M_p));
602:   PetscCall(DMRestoreGlobalVector(dm, &fhat));
603:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
604:   PetscCall(DMRestoreGlobalVector(dm, &grad));
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: static PetscErrorCode TestL2Projection_Shape(DM dm, DM sw, AppCtx *user)
609: {
610:   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
611:   KSP       ksp;
612:   Mat       mass;
613:   Vec       u, rhs, uproj;
614:   void     *ctxs[1];
615:   PetscReal error;

617:   PetscFunctionBeginUser;
618:   ctxs[0]  = (void *)user;
619:   funcs[0] = linear;

621:   PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &u));
622:   PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
623:   PetscCall(DMGetGlobalVector(dm, &rhs));
624:   PetscCall(DMCreateMassMatrix(sw, dm, &mass));
625:   PetscCall(MatMultTranspose(mass, u, rhs));
626:   PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
627:   PetscCall(MatDestroy(&mass));
628:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &u));
629:   PetscCall(DMGetGlobalVector(dm, &uproj));
630:   PetscCall(DMCreateMatrix(dm, &mass));
631:   PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
632:   PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
633:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
634:   PetscCall(KSPSetOperators(ksp, mass, mass));
635:   PetscCall(KSPSetFromOptions(ksp));
636:   PetscCall(KSPSolve(ksp, rhs, uproj));
637:   PetscCall(KSPDestroy(&ksp));
638:   PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
639:   PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
640:   PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
641:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
642:   PetscCall(MatDestroy(&mass));
643:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
644:   PetscCall(DMRestoreGlobalVector(dm, &uproj));
645:   PetscFunctionReturn(PETSC_SUCCESS);
646: }

648: int main(int argc, char *argv[])
649: {
650:   MPI_Comm comm;
651:   DM       dm, sw;
652:   AppCtx   user;

654:   PetscFunctionBeginUser;
655:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
656:   comm = PETSC_COMM_WORLD;
657:   PetscCall(ProcessOptions(comm, &user));
658:   PetscCall(CreateMesh(comm, &dm, &user));
659:   PetscCall(CreateFEM(dm, &user));
660:   if (!user.shape) {
661:     Vec fhat;

663:     PetscCall(CreateParticles(dm, &sw, &user));
664:     PetscCall(DMGetGlobalVector(dm, &fhat));
665:     PetscCall(TestL2ProjectionParticlesToField(dm, sw, fhat, &user));
666:     PetscCall(TestL2ProjectionFieldToParticles(dm, sw, fhat, &user));
667:     PetscCall(TestFieldGradientProjection(dm, sw, &user));
668:     PetscCall(DMRestoreGlobalVector(dm, &fhat));
669:   } else {
670:     PetscCall(CreateParticles_Shape(dm, &sw, &user));
671:     PetscCall(TestL2Projection_Shape(dm, sw, &user));
672:   }
673:   PetscCall(DMDestroy(&dm));
674:   PetscCall(DMDestroy(&sw));
675:   PetscCall(PetscFinalize());
676:   return 0;
677: }

679: /*TEST

681:   # Swarm does not handle complex or quad
682:   build:
683:     requires: !complex double

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

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

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

703:   test:
704:     suffix: proj_quad_2_faces
705:     requires: triangle
706:     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
707:     filter: grep -v marker | grep -v atomic | grep -v usage
708:   test:
709:     suffix: proj_quad_hip_serial
710:     requires: triangle hip
711:     nsize: 1
712:     args: -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -dm_vec_type hip -mat_type aijhipsparse -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1e-15 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
713:     filter: grep -v marker | grep -v atomic | grep -v usage
714:   test:
715:     suffix: proj_quad_hip_parallel
716:     requires: triangle hip
717:     nsize: 2
718:     args: -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -dm_vec_type hip -mat_type aijhipsparse -ptof_ksp_type cg -ptof_pc_type asm -ptof_sub_pc_type ilu -ptof_ksp_rtol 1e-15 -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
719:     filter: grep -v marker | grep -v atomic | grep -v usage
720:   test:
721:     suffix: proj_tri_5P
722:     requires: triangle
723:     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
724:     filter: grep -v marker | grep -v atomic | grep -v usage

726:   test:
727:     suffix: proj_quad_5P
728:     requires: triangle
729:     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
730:     filter: grep -v marker | grep -v atomic | grep -v usage

732:   test:
733:     suffix: proj_tri_mdx
734:     requires: triangle
735:     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
736:     filter: grep -v marker | grep -v atomic | grep -v usage

738:   test:
739:     suffix: proj_tri_mdx_5P
740:     requires: triangle
741:     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
742:     filter: grep -v marker | grep -v atomic | grep -v usage

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

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

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

762:   test:
763:     suffix: proj_tri_3d_mdx
764:     requires: ctetgen
765:     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
766:     filter: grep -v marker | grep -v atomic | grep -v usage

768:   test:
769:     suffix: proj_tri_3d_mdx_5P
770:     requires: ctetgen
771:     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
772:     filter: grep -v marker | grep -v atomic | grep -v usage

774:   test:
775:     suffix: proj_tri_3d_mdx_2_faces
776:     requires: ctetgen
777:     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
778:     filter: grep -v marker | grep -v atomic | grep -v usage

780:   test:
781:     suffix: proj_tri_3d_mdx_5P_2_faces
782:     requires: ctetgen
783:     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
784:     filter: grep -v marker | grep -v atomic | grep -v usage

786:   test:
787:     suffix: proj_quad_lsqr_scale
788:     requires: !complex
789:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
790:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
791:       -particlesPerCell 200 \
792:       -ptof_pc_type lu \
793:       -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none

795:   test:
796:     suffix: proj_quad_lsqr_prec_scale
797:     requires: !complex
798:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
799:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
800:       -particlesPerCell 200 \
801:       -ptof_pc_type lu \
802:       -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
803:   test:
804:     suffix: proj_shape_linear_tri_2d
805:     requires: ctetgen triangle
806:     args: -shape -dm_plex_box_faces 2,2\
807:           -dm_view -sw_view\
808:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
809:           -pc_type lu
810:   test:
811:     suffix: proj_shape_linear_quad_2d
812:     requires: ctetgen triangle
813:     args: -shape -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
814:           -dm_view -sw_view\
815:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
816:           -pc_type lu
817:   test:
818:     suffix: proj_shape_linear_tri_3d
819:     requires: ctetgen triangle
820:     args: -shape -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
821:           -dm_view -sw_view\
822:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
823:           -pc_type lu
824:   test:
825:     suffix: proj_shape_linear_quad_3d
826:     requires: ctetgen triangle
827:     args: -shape -dm_plex_dim 3 -dm_plex_simplex 0\
828:           -dm_plex_box_faces 2,2,2\
829:           -dm_view -sw_view\
830:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
831:           -pc_type lu
832: TEST*/