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;
57: options->particlesPerCell = 1;
58: options->k = 1;
59: options->particleRelDx = 1.e-20;
60: options->meshRelDx = 1.e-20;
61: options->momentTol = 100.*PETSC_MACHINE_EPSILON;
62: options->useBlockDiagPrec = PETSC_FALSE;
64: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
65: PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
66: PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
67: PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
68: PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
69: PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL);
70: PetscStrcmp(fstring, "linear", &flag);
71: if (flag) {
72: options->func = linear;
73: } else {
74: PetscStrcmp(fstring, "sin", &flag);
75: if (flag) {
76: options->func = sinx;
77: } else {
78: PetscStrcmp(fstring, "x2_x4", &flag);
79: options->func = x2_x4;
80: if (!flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown function %s",fstring);
81: }
82: }
83: PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL);
84: PetscOptionsBool("-block_diag_prec", "Use the block diagonal of the normal equations to precondition the particle projection", "ex2.c", options->useBlockDiagPrec, &options->useBlockDiagPrec, NULL);
85: PetscOptionsEnd();
87: return(0);
88: }
90: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
91: {
92: PetscRandom rnd;
93: PetscReal interval = user->meshRelDx;
94: Vec coordinates;
95: PetscScalar *coords;
96: PetscReal *hh, low[3], high[3];
97: PetscInt d, cdim, cEnd, N, p, bs;
101: DMGetBoundingBox(dm, low, high);
102: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
103: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
104: PetscRandomSetInterval(rnd, -interval, interval);
105: PetscRandomSetFromOptions(rnd);
106: DMGetCoordinatesLocal(dm, &coordinates);
107: DMGetCoordinateDim(dm, &cdim);
108: PetscCalloc1(cdim,&hh);
109: for (d = 0; d < cdim; ++d) hh[d] = (user->L[d])/PetscPowReal(cEnd, 1./cdim);
110: VecGetLocalSize(coordinates, &N);
111: VecGetBlockSize(coordinates, &bs);
112: if (bs != cdim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %D != %D", bs, cdim);
113: VecGetArray(coordinates, &coords);
114: for (p = 0; p < N; p += cdim) {
115: PetscScalar *coord = &coords[p], value;
117: for (d = 0; d < cdim; ++d) {
118: PetscRandomGetValue(rnd, &value);
119: coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value*hh[d])));
120: }
121: }
122: VecRestoreArray(coordinates, &coords);
123: PetscRandomDestroy(&rnd);
124: PetscFree(hh);
125: return(0);
126: }
128: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
129: {
130: PetscReal low[3], high[3];
131: PetscInt cdim, d;
135: DMCreate(comm, dm);
136: DMSetType(*dm, DMPLEX);
137: DMSetFromOptions(*dm);
138: DMViewFromOptions(*dm, NULL, "-dm_view");
140: DMGetCoordinateDim(*dm, &cdim);
141: DMGetBoundingBox(*dm, low, high);
142: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
143: PerturbVertices(*dm, user);
144: return(0);
145: }
147: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
148: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
149: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
150: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
151: {
152: g0[0] = 1.0;
153: }
155: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
156: {
157: PetscFE fe;
158: PetscDS ds;
159: DMPolytopeType ct;
160: PetscBool simplex;
161: PetscInt dim, cStart;
165: DMGetDimension(dm, &dim);
166: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
167: DMPlexGetCellType(dm, cStart, &ct);
168: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
169: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, NULL, -1, &fe);
170: PetscObjectSetName((PetscObject) fe, "fe");
171: DMSetField(dm, 0, NULL, (PetscObject) fe);
172: DMCreateDS(dm);
173: PetscFEDestroy(&fe);
174: /* Setup to form mass matrix */
175: DMGetDS(dm, &ds);
176: PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
177: return(0);
178: }
180: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
181: {
182: PetscRandom rnd, rndp;
183: PetscReal interval = user->particleRelDx;
184: DMPolytopeType ct;
185: PetscBool simplex;
186: PetscScalar value, *vals;
187: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
188: PetscInt *cellid;
189: PetscInt Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;
193: DMGetDimension(dm, &dim);
194: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
195: DMPlexGetCellType(dm, cStart, &ct);
196: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
197: DMCreate(PetscObjectComm((PetscObject) dm), sw);
198: DMSetType(*sw, DMSWARM);
199: DMSetDimension(*sw, dim);
201: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
202: PetscRandomSetInterval(rnd, -1.0, 1.0);
203: PetscRandomSetFromOptions(rnd);
204: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
205: PetscRandomSetInterval(rndp, -interval, interval);
206: PetscRandomSetFromOptions(rndp);
208: DMSwarmSetType(*sw, DMSWARM_PIC);
209: DMSwarmSetCellDM(*sw, dm);
210: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
211: DMSwarmFinalizeFieldRegister(*sw);
212: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
213: DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
214: DMSetFromOptions(*sw);
215: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
216: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
217: DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **) &vals);
219: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
220: for (c = 0; c < Ncell; ++c) {
221: if (Np == 1) {
222: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
223: cellid[c] = c;
224: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
225: } else {
226: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
227: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
228: for (p = 0; p < Np; ++p) {
229: const PetscInt n = c*Np + p;
230: PetscReal sum = 0.0, refcoords[3];
232: cellid[n] = c;
233: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rnd, &value); refcoords[d] = PetscRealPart(value); sum += refcoords[d];}
234: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
235: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
236: }
237: }
238: }
239: PetscFree5(centroid, xi0, v0, J, invJ);
240: for (c = 0; c < Ncell; ++c) {
241: for (p = 0; p < Np; ++p) {
242: const PetscInt n = c*Np + p;
244: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value); coords[n*dim+d] += PetscRealPart(value);}
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;
265: PetscErrorCode ierr;
268: DMGetDimension(sw, &dim);
269: DMSwarmGetCellDM(sw, &dm);
270: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
271: DMSwarmSortGetAccess(sw);
272: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
273: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **) &w);
274: for (cell = cStart; cell < cEnd; ++cell) {
275: PetscInt *pidx;
276: PetscInt Np, p, d;
278: DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
279: for (p = 0; p < Np; ++p) {
280: const PetscInt idx = pidx[p];
281: const PetscReal *c = &coords[idx*dim];
283: mom[0] += PetscRealPart(w[idx]);
284: mom[1] += PetscRealPart(w[idx]) * c[0];
285: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d]*c[d];
286: }
287: PetscFree(pidx);
288: }
289: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
290: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &w);
291: DMSwarmSortRestoreAccess(sw);
292: MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) sw));
293: return(0);
294: }
296: static void f0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux,
297: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
298: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
299: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
300: {
301: f0[0] = u[0];
302: }
304: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
305: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
306: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
307: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
308: {
309: f0[0] = x[0]*u[0];
310: }
312: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
313: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
314: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
315: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
316: {
317: PetscInt d;
319: f0[0] = 0.0;
320: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
321: }
323: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
324: {
325: PetscDS prob;
327: PetscScalar mom;
330: DMGetDS(dm, &prob);
331: PetscDSSetObjective(prob, 0, &f0_1);
332: DMPlexComputeIntegralFEM(dm, u, &mom, user);
333: moments[0] = PetscRealPart(mom);
334: PetscDSSetObjective(prob, 0, &f0_x);
335: DMPlexComputeIntegralFEM(dm, u, &mom, user);
336: moments[1] = PetscRealPart(mom);
337: PetscDSSetObjective(prob, 0, &f0_r2);
338: DMPlexComputeIntegralFEM(dm, u, &mom, user);
339: moments[2] = PetscRealPart(mom);
340: return(0);
341: }
343: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
344: {
345: MPI_Comm comm;
346: KSP ksp;
347: Mat M; /* FEM mass matrix */
348: Mat M_p; /* Particle mass matrix */
349: Vec f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
350: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
351: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
352: PetscInt m;
356: PetscObjectGetComm((PetscObject) dm, &comm);
357: KSPCreate(comm, &ksp);
358: KSPSetOptionsPrefix(ksp, "ptof_");
359: DMGetGlobalVector(dm, &fhat);
360: DMGetGlobalVector(dm, &rhs);
362: DMCreateMassMatrix(sw, dm, &M_p);
363: MatViewFromOptions(M_p, NULL, "-M_p_view");
365: /* make particle weight vector */
366: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
368: /* create matrix RHS vector */
369: MatMultTranspose(M_p, f, rhs);
370: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
371: PetscObjectSetName((PetscObject) rhs,"rhs");
372: VecViewFromOptions(rhs, NULL, "-rhs_view");
374: DMCreateMatrix(dm, &M);
375: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
376: MatViewFromOptions(M, NULL, "-M_view");
377: KSPSetOperators(ksp, M, M);
378: KSPSetFromOptions(ksp);
379: KSPSolve(ksp, rhs, fhat);
380: PetscObjectSetName((PetscObject) fhat,"fhat");
381: VecViewFromOptions(fhat, NULL, "-fhat_view");
383: /* Check moments of field */
384: computeParticleMoments(sw, pmoments, user);
385: computeFEMMoments(dm, fhat, fmoments, user);
386: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
387: for (m = 0; m < 3; ++m) {
388: if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
389: }
391: KSPDestroy(&ksp);
392: MatDestroy(&M);
393: MatDestroy(&M_p);
394: DMRestoreGlobalVector(dm, &fhat);
395: DMRestoreGlobalVector(dm, &rhs);
397: return(0);
398: }
400: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
401: {
403: MPI_Comm comm;
404: KSP ksp;
405: Mat M; /* FEM mass matrix */
406: Mat M_p, PM_p; /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
407: Vec f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
408: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
409: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
410: PetscInt m;
414: PetscObjectGetComm((PetscObject) dm, &comm);
416: KSPCreate(comm, &ksp);
417: KSPSetOptionsPrefix(ksp, "ftop_");
418: KSPSetFromOptions(ksp);
420: DMGetGlobalVector(dm, &fhat);
421: DMGetGlobalVector(dm, &rhs);
423: DMCreateMassMatrix(sw, dm, &M_p);
424: MatViewFromOptions(M_p, NULL, "-M_p_view");
426: /* make particle weight vector */
427: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
429: /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
430: PetscObjectSetName((PetscObject) rhs,"rhs");
431: VecViewFromOptions(rhs, NULL, "-rhs_view");
432: DMCreateMatrix(dm, &M);
433: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
434: MatViewFromOptions(M, NULL, "-M_view");
435: MatMultTranspose(M, fhat, rhs);
436: if (user->useBlockDiagPrec) {DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);}
437: else {PetscObjectReference((PetscObject) M_p); PM_p = M_p;}
439: KSPSetOperators(ksp, M_p, PM_p);
440: KSPSolveTranspose(ksp, rhs, f);
441: PetscObjectSetName((PetscObject) fhat,"fhat");
442: VecViewFromOptions(fhat, NULL, "-fhat_view");
444: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
446: /* Check moments */
447: computeParticleMoments(sw, pmoments, user);
448: computeFEMMoments(dm, fhat, fmoments, user);
449: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
450: for (m = 0; m < 3; ++m) {
451: if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
452: }
453: KSPDestroy(&ksp);
454: MatDestroy(&M);
455: MatDestroy(&M_p);
456: MatDestroy(&PM_p);
457: DMRestoreGlobalVector(dm, &fhat);
458: DMRestoreGlobalVector(dm, &rhs);
459: return(0);
460: }
462: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
463: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
464: {
465: DM_Plex *mesh = (DM_Plex *) dm->data;
466: PetscInt debug = mesh->printFEM;
467: DM dmC;
468: PetscSection section;
469: PetscQuadrature quad = NULL;
470: PetscScalar *interpolant, *gradsum;
471: PetscFEGeom fegeom;
472: PetscReal *coords;
473: const PetscReal *quadPoints, *quadWeights;
474: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
475: PetscErrorCode ierr;
478: VecGetDM(locC, &dmC);
479: VecSet(locC, 0.0);
480: DMGetDimension(dm, &dim);
481: DMGetCoordinateDim(dm, &coordDim);
482: fegeom.dimEmbed = coordDim;
483: DMGetLocalSection(dm, §ion);
484: PetscSectionGetNumFields(section, &numFields);
485: for (field = 0; field < numFields; ++field) {
486: PetscObject obj;
487: PetscClassId id;
488: PetscInt Nc;
490: DMGetField(dm, field, NULL, &obj);
491: PetscObjectGetClassId(obj, &id);
492: if (id == PETSCFE_CLASSID) {
493: PetscFE fe = (PetscFE) obj;
495: PetscFEGetQuadrature(fe, &quad);
496: PetscFEGetNumComponents(fe, &Nc);
497: } else if (id == PETSCFV_CLASSID) {
498: PetscFV fv = (PetscFV) obj;
500: PetscFVGetQuadrature(fv, &quad);
501: PetscFVGetNumComponents(fv, &Nc);
502: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
503: numComponents += Nc;
504: }
505: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
506: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
507: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
508: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
509: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
510: for (v = vStart; v < vEnd; ++v) {
511: PetscInt *star = NULL;
512: PetscInt starSize, st, d, fc;
514: PetscArrayzero(gradsum, coordDim*numComponents);
515: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
516: for (st = 0; st < starSize*2; st += 2) {
517: const PetscInt cell = star[st];
518: PetscScalar *grad = &gradsum[coordDim*numComponents];
519: PetscScalar *x = NULL;
521: if ((cell < cStart) || (cell >= cEnd)) continue;
522: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
523: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
524: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
525: PetscObject obj;
526: PetscClassId id;
527: PetscInt Nb, Nc, q, qc = 0;
529: PetscArrayzero(grad, coordDim*numComponents);
530: DMGetField(dm, field, NULL, &obj);
531: PetscObjectGetClassId(obj, &id);
532: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
533: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
534: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
535: for (q = 0; q < Nq; ++q) {
536: if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], cell, q);
537: if (ierr) {
538: PetscErrorCode ierr2;
539: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
540: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
541: ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
542:
543: }
544: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, 1, &x[fieldOffset], &fegeom, q, interpolant);}
545: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
546: for (fc = 0; fc < Nc; ++fc) {
547: const PetscReal wt = quadWeights[q*qNc+qc+fc];
549: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
550: }
551: }
552: fieldOffset += Nb;
553: }
554: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
555: for (fc = 0; fc < numComponents; ++fc) {
556: for (d = 0; d < coordDim; ++d) {
557: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
558: }
559: }
560: if (debug) {
561: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
562: for (fc = 0; fc < numComponents; ++fc) {
563: for (d = 0; d < coordDim; ++d) {
564: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
565: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
566: }
567: }
568: PetscPrintf(PETSC_COMM_SELF, "]\n");
569: }
570: }
571: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
572: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
573: }
574: PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
575: return(0);
576: }
578: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
579: {
581: MPI_Comm comm;
582: KSP ksp;
583: Mat M; /* FEM mass matrix */
584: Mat M_p; /* Particle mass matrix */
585: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
586: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
587: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
588: PetscInt m;
592: PetscObjectGetComm((PetscObject) dm, &comm);
593: KSPCreate(comm, &ksp);
594: KSPSetOptionsPrefix(ksp, "ptof_");
595: DMGetGlobalVector(dm, &fhat);
596: DMGetGlobalVector(dm, &rhs);
597: DMGetGlobalVector(dm, &grad);
599: DMCreateMassMatrix(sw, dm, &M_p);
600: MatViewFromOptions(M_p, NULL, "-M_p_view");
602: /* make particle weight vector */
603: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
605: /* create matrix RHS vector */
606: MatMultTranspose(M_p, f, rhs);
607: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
608: PetscObjectSetName((PetscObject) rhs,"rhs");
609: VecViewFromOptions(rhs, NULL, "-rhs_view");
611: DMCreateMatrix(dm, &M);
612: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
614: InterpolateGradient(dm, fhat, grad);
616: MatViewFromOptions(M, NULL, "-M_view");
617: KSPSetOperators(ksp, M, M);
618: KSPSetFromOptions(ksp);
619: KSPSolve(ksp, rhs, grad);
620: PetscObjectSetName((PetscObject) fhat,"fhat");
621: VecViewFromOptions(fhat, NULL, "-fhat_view");
623: /* Check moments of field */
624: computeParticleMoments(sw, pmoments, user);
625: computeFEMMoments(dm, grad, fmoments, user);
626: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
627: for (m = 0; m < 3; ++m) {
628: if (PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]) > user->momentTol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Moment %D error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m])/fmoments[m]), user->momentTol);
629: }
631: KSPDestroy(&ksp);
632: MatDestroy(&M);
633: MatDestroy(&M_p);
634: DMRestoreGlobalVector(dm, &fhat);
635: DMRestoreGlobalVector(dm, &rhs);
636: DMRestoreGlobalVector(dm, &grad);
638: return(0);
639: }
641: int main (int argc, char * argv[]) {
642: MPI_Comm comm;
643: DM dm, sw;
644: AppCtx user;
647: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
648: comm = PETSC_COMM_WORLD;
649: ProcessOptions(comm, &user);
650: CreateMesh(comm, &dm, &user);
651: CreateFEM(dm, &user);
652: CreateParticles(dm, &sw, &user);
653: TestL2ProjectionParticlesToField(dm, sw, &user);
654: TestL2ProjectionFieldToParticles(dm, sw, &user);
655: TestFieldGradientProjection(dm, sw, &user);
656: DMDestroy(&dm);
657: DMDestroy(&sw);
658: PetscFinalize();
659: return ierr;
660: }
662: /*TEST
664: # Swarm does not handle complex or quad
665: build:
666: requires: !complex double
668: test:
669: suffix: proj_tri_0
670: requires: triangle
671: 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
672: filter: grep -v marker | grep -v atomic | grep -v usage
674: test:
675: suffix: proj_tri_2_faces
676: requires: triangle
677: 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
678: filter: grep -v marker | grep -v atomic | grep -v usage
680: test:
681: suffix: proj_quad_0
682: requires: triangle
683: 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
684: filter: grep -v marker | grep -v atomic | grep -v usage
686: test:
687: suffix: proj_quad_2_faces
688: requires: triangle
689: 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
690: filter: grep -v marker | grep -v atomic | grep -v usage
692: test:
693: suffix: proj_tri_5P
694: requires: triangle
695: 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
696: filter: grep -v marker | grep -v atomic | grep -v usage
698: test:
699: suffix: proj_quad_5P
700: requires: triangle
701: 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
702: filter: grep -v marker | grep -v atomic | grep -v usage
704: test:
705: suffix: proj_tri_mdx
706: requires: triangle
707: 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
708: filter: grep -v marker | grep -v atomic | grep -v usage
710: test:
711: suffix: proj_tri_mdx_5P
712: requires: triangle
713: 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
714: filter: grep -v marker | grep -v atomic | grep -v usage
716: test:
717: suffix: proj_tri_3d
718: requires: ctetgen
719: 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
720: filter: grep -v marker | grep -v atomic | grep -v usage
722: test:
723: suffix: proj_tri_3d_2_faces
724: requires: ctetgen
725: 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
726: filter: grep -v marker | grep -v atomic | grep -v usage
728: test:
729: suffix: proj_tri_3d_5P
730: requires: ctetgen
731: 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
732: filter: grep -v marker | grep -v atomic | grep -v usage
734: test:
735: suffix: proj_tri_3d_mdx
736: requires: ctetgen
737: 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
738: filter: grep -v marker | grep -v atomic | grep -v usage
740: test:
741: suffix: proj_tri_3d_mdx_5P
742: requires: ctetgen
743: 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
744: filter: grep -v marker | grep -v atomic | grep -v usage
746: test:
747: suffix: proj_tri_3d_mdx_2_faces
748: requires: ctetgen
749: 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
750: filter: grep -v marker | grep -v atomic | grep -v usage
752: test:
753: suffix: proj_tri_3d_mdx_5P_2_faces
754: requires: ctetgen
755: 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
756: filter: grep -v marker | grep -v atomic | grep -v usage
758: test:
759: suffix: proj_quad_lsqr_scale
760: requires: !complex
761: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
762: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
763: -particlesPerCell 200 \
764: -ptof_pc_type lu \
765: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none
767: test:
768: suffix: proj_quad_lsqr_prec_scale
769: requires: !complex
770: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
771: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
772: -particlesPerCell 200 \
773: -ptof_pc_type lu \
774: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type lu -ftop_pc_factor_shift_type nonzero -block_diag_prec
776: TEST*/