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:   PetscRandom    rnd, rndp;
173:   PetscReal      interval = user->particleRelDx;
174:   DMPolytopeType ct;
175:   PetscBool      simplex;
176:   PetscScalar    value, *vals;
177:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
178:   PetscInt      *cellid;
179:   PetscInt       Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

663: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

768:   test:
769:     suffix: proj_quad_lsqr_prec_scale
770:     requires: !complex
771:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
772:       -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
773:       -particlesPerCell 200 \
774:       -ptof_pc_type lu \
775:       -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
776:   test:
777:     suffix: proj_shape_linear_tri_2d
778:     requires: ctetgen triangle
779:     args: -shape -dm_plex_box_faces 2,2\
780:           -dm_view -sw_view\
781:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
782:           -pc_type lu
783:   test:
784:     suffix: proj_shape_linear_quad_2d
785:     requires: ctetgen triangle
786:     args: -shape -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
787:           -dm_view -sw_view\
788:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
789:           -pc_type lu
790:   test:
791:     suffix: proj_shape_linear_tri_3d
792:     requires: ctetgen triangle
793:     args: -shape -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
794:           -dm_view -sw_view\
795:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
796:           -pc_type lu
797:   test:
798:     suffix: proj_shape_linear_quad_3d
799:     requires: ctetgen triangle
800:     args: -shape -dm_plex_dim 3 -dm_plex_simplex 0\
801:           -dm_plex_box_faces 2,2,2\
802:           -dm_view -sw_view\
803:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
804:           -pc_type lu
805: TEST*/