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, ¢roid, 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, §ion));
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*/