Actual source code: ex2.c
1: static char help[] = "Tests L2 projection with DMSwarm using delta function particles.\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 useBlockDiagPrec; /* Use the block diagonal of the normal equations as a preconditioner */
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 0;
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 0;
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 0;
48: }
50: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
51: {
52: char fstring[PETSC_MAX_PATH_LEN] = "linear";
53: PetscBool flag;
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->useBlockDiagPrec = PETSC_FALSE;
63: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
64: PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
65: PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
66: PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
67: PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
68: PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL);
69: PetscStrcmp(fstring, "linear", &flag);
70: if (flag) {
71: options->func = linear;
72: } else {
73: PetscStrcmp(fstring, "sin", &flag);
74: if (flag) {
75: options->func = sinx;
76: } else {
77: PetscStrcmp(fstring, "x2_x4", &flag);
78: options->func = x2_x4;
80: }
81: }
82: PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL);
83: PetscOptionsBool("-block_diag_prec", "Use the block diagonal of the normal equations to precondition the particle projection", "ex2.c", options->useBlockDiagPrec, &options->useBlockDiagPrec, NULL);
84: PetscOptionsEnd();
86: return 0;
87: }
89: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
90: {
91: PetscRandom rnd;
92: PetscReal interval = user->meshRelDx;
93: Vec coordinates;
94: PetscScalar *coords;
95: PetscReal *hh, low[3], high[3];
96: PetscInt d, cdim, cEnd, N, p, bs;
99: DMGetBoundingBox(dm, low, high);
100: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
101: PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd);
102: PetscRandomSetInterval(rnd, -interval, interval);
103: PetscRandomSetFromOptions(rnd);
104: DMGetCoordinatesLocal(dm, &coordinates);
105: DMGetCoordinateDim(dm, &cdim);
106: PetscCalloc1(cdim, &hh);
107: for (d = 0; d < cdim; ++d) hh[d] = (user->L[d]) / PetscPowReal(cEnd, 1. / cdim);
108: VecGetLocalSize(coordinates, &N);
109: VecGetBlockSize(coordinates, &bs);
111: VecGetArray(coordinates, &coords);
112: for (p = 0; p < N; p += cdim) {
113: PetscScalar *coord = &coords[p], value;
115: for (d = 0; d < cdim; ++d) {
116: PetscRandomGetValue(rnd, &value);
117: coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value * hh[d])));
118: }
119: }
120: VecRestoreArray(coordinates, &coords);
121: PetscRandomDestroy(&rnd);
122: PetscFree(hh);
123: return 0;
124: }
126: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
127: {
128: PetscReal low[3], high[3];
129: PetscInt cdim, d;
132: DMCreate(comm, dm);
133: DMSetType(*dm, DMPLEX);
134: DMSetFromOptions(*dm);
135: DMViewFromOptions(*dm, NULL, "-dm_view");
137: DMGetCoordinateDim(*dm, &cdim);
138: DMGetBoundingBox(*dm, low, high);
139: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
140: PerturbVertices(*dm, user);
141: return 0;
142: }
144: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
145: {
146: g0[0] = 1.0;
147: }
149: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
150: {
151: PetscFE fe;
152: PetscDS ds;
153: DMPolytopeType ct;
154: PetscBool simplex;
155: PetscInt dim, cStart;
158: DMGetDimension(dm, &dim);
159: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
160: DMPlexGetCellType(dm, cStart, &ct);
161: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
162: PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe);
163: PetscObjectSetName((PetscObject)fe, "fe");
164: DMSetField(dm, 0, NULL, (PetscObject)fe);
165: DMCreateDS(dm);
166: PetscFEDestroy(&fe);
167: /* Setup to form mass matrix */
168: DMGetDS(dm, &ds);
169: PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
170: return 0;
171: }
173: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
174: {
175: PetscRandom rnd, rndp;
176: PetscReal interval = user->particleRelDx;
177: DMPolytopeType ct;
178: PetscBool simplex;
179: PetscScalar value, *vals;
180: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
181: PetscInt *cellid;
182: PetscInt Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;
185: DMGetDimension(dm, &dim);
186: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
187: DMPlexGetCellType(dm, cStart, &ct);
188: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
189: DMCreate(PetscObjectComm((PetscObject)dm), sw);
190: DMSetType(*sw, DMSWARM);
191: DMSetDimension(*sw, dim);
193: PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd);
194: PetscRandomSetInterval(rnd, -1.0, 1.0);
195: PetscRandomSetFromOptions(rnd);
196: PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rndp);
197: PetscRandomSetInterval(rndp, -interval, interval);
198: PetscRandomSetFromOptions(rndp);
200: DMSwarmSetType(*sw, DMSWARM_PIC);
201: DMSwarmSetCellDM(*sw, dm);
202: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
203: DMSwarmFinalizeFieldRegister(*sw);
204: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
205: DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
206: DMSetFromOptions(*sw);
207: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
208: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
209: DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals);
211: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
212: for (c = 0; c < Ncell; ++c) {
213: if (Np == 1) {
214: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
215: cellid[c] = c;
216: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
217: } else {
218: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
219: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
220: for (p = 0; p < Np; ++p) {
221: const PetscInt n = c * Np + p;
222: PetscReal sum = 0.0, refcoords[3];
224: cellid[n] = c;
225: for (d = 0; d < dim; ++d) {
226: PetscRandomGetValue(rnd, &value);
227: refcoords[d] = PetscRealPart(value);
228: sum += refcoords[d];
229: }
230: if (simplex && sum > 0.0)
231: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
232: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
233: }
234: }
235: }
236: PetscFree5(centroid, xi0, v0, J, invJ);
237: for (c = 0; c < Ncell; ++c) {
238: for (p = 0; p < Np; ++p) {
239: const PetscInt n = c * Np + p;
241: for (d = 0; d < dim; ++d) {
242: PetscRandomGetValue(rndp, &value);
243: coords[n * dim + d] += PetscRealPart(value);
244: }
245: user->func(dim, 0.0, &coords[n * dim], 1, &vals[c], user);
246: }
247: }
248: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
249: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
250: DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals);
251: PetscRandomDestroy(&rnd);
252: PetscRandomDestroy(&rndp);
253: PetscObjectSetName((PetscObject)*sw, "Particles");
254: DMViewFromOptions(*sw, NULL, "-sw_view");
255: return 0;
256: }
258: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
259: {
260: DM dm;
261: const PetscReal *coords;
262: const PetscScalar *w;
263: PetscReal mom[3] = {0.0, 0.0, 0.0};
264: PetscInt cell, cStart, cEnd, dim;
267: DMGetDimension(sw, &dim);
268: DMSwarmGetCellDM(sw, &dm);
269: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
270: DMSwarmSortGetAccess(sw);
271: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
272: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w);
273: for (cell = cStart; cell < cEnd; ++cell) {
274: PetscInt *pidx;
275: PetscInt Np, p, d;
277: DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
278: for (p = 0; p < Np; ++p) {
279: const PetscInt idx = pidx[p];
280: const PetscReal *c = &coords[idx * dim];
282: mom[0] += PetscRealPart(w[idx]);
283: mom[1] += PetscRealPart(w[idx]) * c[0];
284: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
285: }
286: PetscFree(pidx);
287: }
288: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
289: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w);
290: DMSwarmSortRestoreAccess(sw);
291: MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw));
292: return 0;
293: }
295: 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[])
296: {
297: f0[0] = u[0];
298: }
300: 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[])
301: {
302: f0[0] = x[0] * u[0];
303: }
305: 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[])
306: {
307: PetscInt d;
309: f0[0] = 0.0;
310: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
311: }
313: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
314: {
315: PetscDS prob;
316: PetscScalar mom;
319: DMGetDS(dm, &prob);
320: PetscDSSetObjective(prob, 0, &f0_1);
321: DMPlexComputeIntegralFEM(dm, u, &mom, user);
322: moments[0] = PetscRealPart(mom);
323: PetscDSSetObjective(prob, 0, &f0_x);
324: DMPlexComputeIntegralFEM(dm, u, &mom, user);
325: moments[1] = PetscRealPart(mom);
326: PetscDSSetObjective(prob, 0, &f0_r2);
327: DMPlexComputeIntegralFEM(dm, u, &mom, user);
328: moments[2] = PetscRealPart(mom);
329: return 0;
330: }
332: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
333: {
334: MPI_Comm comm;
335: KSP ksp;
336: Mat M; /* FEM mass matrix */
337: Mat M_p; /* Particle mass matrix */
338: Vec f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
339: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
340: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
341: PetscInt m;
344: PetscObjectGetComm((PetscObject)dm, &comm);
345: KSPCreate(comm, &ksp);
346: KSPSetOptionsPrefix(ksp, "ptof_");
347: DMGetGlobalVector(dm, &fhat);
348: DMGetGlobalVector(dm, &rhs);
350: DMCreateMassMatrix(sw, dm, &M_p);
351: MatViewFromOptions(M_p, NULL, "-M_p_view");
353: /* make particle weight vector */
354: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
356: /* create matrix RHS vector */
357: MatMultTranspose(M_p, f, rhs);
358: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
359: PetscObjectSetName((PetscObject)rhs, "rhs");
360: VecViewFromOptions(rhs, NULL, "-rhs_view");
362: DMCreateMatrix(dm, &M);
363: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
364: MatViewFromOptions(M, NULL, "-M_view");
365: KSPSetOperators(ksp, M, M);
366: KSPSetFromOptions(ksp);
367: KSPSolve(ksp, rhs, fhat);
368: PetscObjectSetName((PetscObject)fhat, "fhat");
369: VecViewFromOptions(fhat, NULL, "-fhat_view");
371: /* Check moments of field */
372: computeParticleMoments(sw, pmoments, user);
373: computeFEMMoments(dm, fhat, fmoments, user);
374: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
375: for (m = 0; m < 3; ++m) {
377: }
379: KSPDestroy(&ksp);
380: MatDestroy(&M);
381: MatDestroy(&M_p);
382: DMRestoreGlobalVector(dm, &fhat);
383: DMRestoreGlobalVector(dm, &rhs);
385: return 0;
386: }
388: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
389: {
390: MPI_Comm comm;
391: KSP ksp;
392: Mat M; /* FEM mass matrix */
393: Mat M_p, PM_p; /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
394: Vec f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
395: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
396: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
397: PetscInt m;
400: PetscObjectGetComm((PetscObject)dm, &comm);
402: KSPCreate(comm, &ksp);
403: KSPSetOptionsPrefix(ksp, "ftop_");
404: KSPSetFromOptions(ksp);
406: DMGetGlobalVector(dm, &fhat);
407: DMGetGlobalVector(dm, &rhs);
409: DMCreateMassMatrix(sw, dm, &M_p);
410: MatViewFromOptions(M_p, NULL, "-M_p_view");
412: /* make particle weight vector */
413: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
415: /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
416: PetscObjectSetName((PetscObject)rhs, "rhs");
417: VecViewFromOptions(rhs, NULL, "-rhs_view");
418: DMCreateMatrix(dm, &M);
419: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
420: MatViewFromOptions(M, NULL, "-M_view");
421: MatMultTranspose(M, fhat, rhs);
422: if (user->useBlockDiagPrec) DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);
423: else {
424: PetscObjectReference((PetscObject)M_p);
425: PM_p = M_p;
426: }
428: KSPSetOperators(ksp, M_p, PM_p);
429: KSPSolveTranspose(ksp, rhs, f);
430: PetscObjectSetName((PetscObject)fhat, "fhat");
431: VecViewFromOptions(fhat, NULL, "-fhat_view");
433: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
435: /* Check moments */
436: computeParticleMoments(sw, pmoments, user);
437: computeFEMMoments(dm, fhat, fmoments, user);
438: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
439: for (m = 0; m < 3; ++m) {
441: }
442: KSPDestroy(&ksp);
443: MatDestroy(&M);
444: MatDestroy(&M_p);
445: MatDestroy(&PM_p);
446: DMRestoreGlobalVector(dm, &fhat);
447: DMRestoreGlobalVector(dm, &rhs);
448: return 0;
449: }
451: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
452: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
453: {
454: DM_Plex *mesh = (DM_Plex *)dm->data;
455: PetscInt debug = mesh->printFEM;
456: DM dmC;
457: PetscSection section;
458: PetscQuadrature quad = NULL;
459: PetscScalar *interpolant, *gradsum;
460: PetscFEGeom fegeom;
461: PetscReal *coords;
462: const PetscReal *quadPoints, *quadWeights;
463: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
465: VecGetDM(locC, &dmC);
466: VecSet(locC, 0.0);
467: DMGetDimension(dm, &dim);
468: DMGetCoordinateDim(dm, &coordDim);
469: fegeom.dimEmbed = coordDim;
470: DMGetLocalSection(dm, §ion);
471: PetscSectionGetNumFields(section, &numFields);
472: for (field = 0; field < numFields; ++field) {
473: PetscObject obj;
474: PetscClassId id;
475: PetscInt Nc;
477: DMGetField(dm, field, NULL, &obj);
478: PetscObjectGetClassId(obj, &id);
479: if (id == PETSCFE_CLASSID) {
480: PetscFE fe = (PetscFE)obj;
482: PetscFEGetQuadrature(fe, &quad);
483: PetscFEGetNumComponents(fe, &Nc);
484: } else if (id == PETSCFV_CLASSID) {
485: PetscFV fv = (PetscFV)obj;
487: PetscFVGetQuadrature(fv, &quad);
488: PetscFVGetNumComponents(fv, &Nc);
489: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
490: numComponents += Nc;
491: }
492: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
494: PetscMalloc6(coordDim * numComponents * 2, &gradsum, coordDim * numComponents, &interpolant, coordDim * Nq, &coords, Nq, &fegeom.detJ, coordDim * coordDim * Nq, &fegeom.J, coordDim * coordDim * Nq, &fegeom.invJ);
495: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
496: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
497: for (v = vStart; v < vEnd; ++v) {
498: PetscInt *star = NULL;
499: PetscInt starSize, st, d, fc;
501: PetscArrayzero(gradsum, coordDim * numComponents);
502: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
503: for (st = 0; st < starSize * 2; st += 2) {
504: const PetscInt cell = star[st];
505: PetscScalar *grad = &gradsum[coordDim * numComponents];
506: PetscScalar *x = NULL;
508: if ((cell < cStart) || (cell >= cEnd)) continue;
509: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
510: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
511: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
512: PetscObject obj;
513: PetscClassId id;
514: PetscInt Nb, Nc, q, qc = 0;
516: PetscArrayzero(grad, coordDim * numComponents);
517: DMGetField(dm, field, NULL, &obj);
518: PetscObjectGetClassId(obj, &id);
519: if (id == PETSCFE_CLASSID) {
520: PetscFEGetNumComponents((PetscFE)obj, &Nc);
521: PetscFEGetDimension((PetscFE)obj, &Nb);
522: } else if (id == PETSCFV_CLASSID) {
523: PetscFVGetNumComponents((PetscFV)obj, &Nc);
524: Nb = 1;
525: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
526: for (q = 0; q < Nq; ++q) {
528: if (id == PETSCFE_CLASSID) PetscFEInterpolateGradient_Static((PetscFE)obj, 1, &x[fieldOffset], &fegeom, q, interpolant);
529: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
530: for (fc = 0; fc < Nc; ++fc) {
531: const PetscReal wt = quadWeights[q * qNc + qc + fc];
533: for (d = 0; d < coordDim; ++d) grad[fc * coordDim + d] += interpolant[fc * dim + d] * wt * fegeom.detJ[q];
534: }
535: }
536: fieldOffset += Nb;
537: }
538: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
539: for (fc = 0; fc < numComponents; ++fc) {
540: for (d = 0; d < coordDim; ++d) gradsum[fc * coordDim + d] += grad[fc * coordDim + d];
541: }
542: if (debug) {
543: PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " gradient: [", cell);
544: for (fc = 0; fc < numComponents; ++fc) {
545: for (d = 0; d < coordDim; ++d) {
546: if (fc || d > 0) PetscPrintf(PETSC_COMM_SELF, ", ");
547: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc * coordDim + d]));
548: }
549: }
550: PetscPrintf(PETSC_COMM_SELF, "]\n");
551: }
552: }
553: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
554: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
555: }
556: PetscFree6(gradsum, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ);
557: return 0;
558: }
560: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
561: {
562: MPI_Comm comm;
563: KSP ksp;
564: Mat M; /* FEM mass matrix */
565: Mat M_p; /* Particle mass matrix */
566: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
567: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
568: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
569: PetscInt m;
572: PetscObjectGetComm((PetscObject)dm, &comm);
573: KSPCreate(comm, &ksp);
574: KSPSetOptionsPrefix(ksp, "ptof_");
575: DMGetGlobalVector(dm, &fhat);
576: DMGetGlobalVector(dm, &rhs);
577: DMGetGlobalVector(dm, &grad);
579: DMCreateMassMatrix(sw, dm, &M_p);
580: MatViewFromOptions(M_p, NULL, "-M_p_view");
582: /* make particle weight vector */
583: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
585: /* create matrix RHS vector */
586: MatMultTranspose(M_p, f, rhs);
587: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
588: PetscObjectSetName((PetscObject)rhs, "rhs");
589: VecViewFromOptions(rhs, NULL, "-rhs_view");
591: DMCreateMatrix(dm, &M);
592: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
594: InterpolateGradient(dm, fhat, grad);
596: MatViewFromOptions(M, NULL, "-M_view");
597: KSPSetOperators(ksp, M, M);
598: KSPSetFromOptions(ksp);
599: KSPSolve(ksp, rhs, grad);
600: PetscObjectSetName((PetscObject)fhat, "fhat");
601: VecViewFromOptions(fhat, NULL, "-fhat_view");
603: /* Check moments of field */
604: computeParticleMoments(sw, pmoments, user);
605: computeFEMMoments(dm, grad, fmoments, user);
606: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
607: for (m = 0; m < 3; ++m) {
609: }
611: KSPDestroy(&ksp);
612: MatDestroy(&M);
613: MatDestroy(&M_p);
614: DMRestoreGlobalVector(dm, &fhat);
615: DMRestoreGlobalVector(dm, &rhs);
616: DMRestoreGlobalVector(dm, &grad);
618: return 0;
619: }
621: int main(int argc, char *argv[])
622: {
623: MPI_Comm comm;
624: DM dm, sw;
625: AppCtx user;
628: PetscInitialize(&argc, &argv, NULL, help);
629: comm = PETSC_COMM_WORLD;
630: ProcessOptions(comm, &user);
631: CreateMesh(comm, &dm, &user);
632: CreateFEM(dm, &user);
633: CreateParticles(dm, &sw, &user);
634: TestL2ProjectionParticlesToField(dm, sw, &user);
635: TestL2ProjectionFieldToParticles(dm, sw, &user);
636: TestFieldGradientProjection(dm, sw, &user);
637: DMDestroy(&dm);
638: DMDestroy(&sw);
639: PetscFinalize();
640: return 0;
641: }
643: /*TEST
645: # Swarm does not handle complex or quad
646: build:
647: requires: !complex double
649: test:
650: suffix: proj_tri_0
651: requires: triangle
652: 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
653: filter: grep -v marker | grep -v atomic | grep -v usage
655: test:
656: suffix: proj_tri_2_faces
657: requires: triangle
658: 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
659: filter: grep -v marker | grep -v atomic | grep -v usage
661: test:
662: suffix: proj_quad_0
663: requires: triangle
664: 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
665: filter: grep -v marker | grep -v atomic | grep -v usage
667: test:
668: suffix: proj_quad_2_faces
669: requires: triangle
670: 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
671: filter: grep -v marker | grep -v atomic | grep -v usage
673: test:
674: suffix: proj_tri_5P
675: requires: triangle
676: 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
677: filter: grep -v marker | grep -v atomic | grep -v usage
679: test:
680: suffix: proj_quad_5P
681: requires: triangle
682: 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
683: filter: grep -v marker | grep -v atomic | grep -v usage
685: test:
686: suffix: proj_tri_mdx
687: requires: triangle
688: 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
689: filter: grep -v marker | grep -v atomic | grep -v usage
691: test:
692: suffix: proj_tri_mdx_5P
693: requires: triangle
694: 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
695: filter: grep -v marker | grep -v atomic | grep -v usage
697: test:
698: suffix: proj_tri_3d
699: requires: ctetgen
700: 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
701: filter: grep -v marker | grep -v atomic | grep -v usage
703: test:
704: suffix: proj_tri_3d_2_faces
705: requires: ctetgen
706: 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
707: filter: grep -v marker | grep -v atomic | grep -v usage
709: test:
710: suffix: proj_tri_3d_5P
711: requires: ctetgen
712: 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
713: filter: grep -v marker | grep -v atomic | grep -v usage
715: test:
716: suffix: proj_tri_3d_mdx
717: requires: ctetgen
718: 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
719: filter: grep -v marker | grep -v atomic | grep -v usage
721: test:
722: suffix: proj_tri_3d_mdx_5P
723: requires: ctetgen
724: 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
725: filter: grep -v marker | grep -v atomic | grep -v usage
727: test:
728: suffix: proj_tri_3d_mdx_2_faces
729: requires: ctetgen
730: 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
731: filter: grep -v marker | grep -v atomic | grep -v usage
733: test:
734: suffix: proj_tri_3d_mdx_5P_2_faces
735: requires: ctetgen
736: 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
737: filter: grep -v marker | grep -v atomic | grep -v usage
739: test:
740: suffix: proj_quad_lsqr_scale
741: requires: !complex
742: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
743: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
744: -particlesPerCell 200 \
745: -ptof_pc_type lu \
746: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none
748: test:
749: suffix: proj_quad_lsqr_prec_scale
750: requires: !complex
751: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
752: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
753: -particlesPerCell 200 \
754: -ptof_pc_type lu \
755: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type lu -ftop_pc_factor_shift_type nonzero -block_diag_prec
757: TEST*/