1 static char help[] = "Tests L2 projection with DMSwarm using delta function particles and deposition of linear shape.\n";
2
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;
19
20 /* const char *const ex2FunctionTypes[] = {"linear","x2_x4","sin","ex2FunctionTypes","EX2_FUNCTION_",0}; */
21
linear(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * a_ctx)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;
26
27 u[0] = 0.0;
28 for (d = 0; d < dim; ++d) u[0] += x[d] / (ctx->L[d]);
29 return PETSC_SUCCESS;
30 }
31
x2_x4(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * a_ctx)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;
36
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 }
41
sinx(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,void * a_ctx)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;
45
46 u[0] = PetscSinScalar(2 * PETSC_PI * ctx->k * x[0] / (ctx->L[0]));
47 return PETSC_SUCCESS;
48 }
49
ProcessOptions(MPI_Comm comm,AppCtx * options)50 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
51 {
52 char fstring[PETSC_MAX_PATH_LEN] = "linear";
53 PetscBool flag;
54
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;
62
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 }
87
PerturbVertices(DM dm,AppCtx * user)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;
96
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;
113
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 }
124
CreateMesh(MPI_Comm comm,DM * dm,AppCtx * user)125 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
126 {
127 PetscReal low[3], high[3];
128 PetscInt cdim, d;
129
130 PetscFunctionBeginUser;
131 PetscCall(DMCreate(comm, dm));
132 PetscCall(DMSetType(*dm, DMPLEX));
133 PetscCall(DMSetFromOptions(*dm));
134 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
135
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 }
142
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[])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 }
147
CreateFEM(DM dm,AppCtx * user)148 static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
149 {
150 PetscFE fe;
151 PetscDS ds;
152 DMPolytopeType ct;
153 PetscInt dim, cStart;
154
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 }
169
CreateParticles(DM dm,DM * sw,AppCtx * user)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;
182
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));
191
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));
198
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));
206
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);
211
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));
215
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];
228
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;
245
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 }
263
CreateParticles_Shape(DM dm,DM * sw,AppCtx * user)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;
276
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 }
321
computeParticleMoments(DM sw,PetscReal moments[3],AppCtx * user)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;
329
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;
340
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];
345
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 }
358
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[])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 }
363
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[])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 }
368
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[])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;
372
373 f0[0] = 0.0;
374 for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
375 }
376
computeFEMMoments(DM dm,Vec u,PetscReal moments[3],AppCtx * user)377 static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
378 {
379 PetscDS prob;
380 PetscScalar mom;
381
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 }
395
TestL2ProjectionParticlesToField(DM dm,DM sw,Vec fhat,AppCtx * user)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
402
403 PetscFunctionBeginUser;
404 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD));
405
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 }
416
TestL2ProjectionFieldToParticles(DM dm,DM sw,Vec fhat,AppCtx * user)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
423
424 PetscFunctionBeginUser;
425 PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE));
426
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 }
437
438 /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
InterpolateGradient(DM dm,Vec locX,Vec locC)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;
451
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;
464
465 PetscCall(DMGetField(dm, field, NULL, &obj));
466 PetscCall(PetscObjectGetClassId(obj, &id));
467 if (id == PETSCFE_CLASSID) {
468 PetscFE fe = (PetscFE)obj;
469
470 PetscCall(PetscFEGetQuadrature(fe, &quad));
471 PetscCall(PetscFEGetNumComponents(fe, &Nc));
472 } else if (id == PETSCFV_CLASSID) {
473 PetscFV fv = (PetscFV)obj;
474
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;
488
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;
495
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;
503
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];
520
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 }
547
TestFieldGradientProjection(DM dm,DM sw,AppCtx * user)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;
558
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));
566
567 PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
568 PetscCall(MatViewFromOptions(M_p, NULL, "-M_p_view"));
569
570 /* make particle weight vector */
571 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
572
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"));
578
579 PetscCall(DMCreateMatrix(dm, &M));
580 PetscCall(DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user));
581
582 PetscCall(InterpolateGradient(dm, fhat, grad));
583
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"));
590
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 }
598
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 }
607
TestL2Projection_Shape(DM dm,DM sw,AppCtx * user)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;
616
617 PetscFunctionBeginUser;
618 ctxs[0] = (void *)user;
619 funcs[0] = linear;
620
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 }
647
main(int argc,char * argv[])648 int main(int argc, char *argv[])
649 {
650 MPI_Comm comm;
651 DM dm, sw;
652 AppCtx user;
653
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;
662
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 }
678
679 /*TEST
680
681 # Swarm does not handle complex or quad
682 build:
683 requires: !complex double
684
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
690
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
696
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
702
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
725
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
731
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
737
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
743
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
749
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
755
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
761
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
767
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
773
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
779
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
785
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
794
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*/
833