Actual source code: ex2.c
petsc-3.14.6 2021-03-30
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: char meshFilename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
11: PetscReal L[3]; /* Dimensions of the mesh bounding box */
12: PetscInt particlesPerCell; /* The number of partices per cell */
13: PetscReal particleRelDx; /* Relative particle position perturbation compared to average cell diameter h */
14: PetscReal meshRelDx; /* Relative vertex position perturbation compared to average cell diameter h */
15: PetscInt k; /* Mode number for test function */
16: PetscReal momentTol; /* Tolerance for checking moment conservation */
17: PetscBool useBlockDiagPrec; /* Use the block diagonal of the normal equations as a preconditioner */
18: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
19: } AppCtx;
21: /* const char *const ex2FunctionTypes[] = {"linear","x2_x4","sin","ex2FunctionTypes","EX2_FUNCTION_",0}; */
23: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
24: {
25: AppCtx *ctx = (AppCtx *) a_ctx;
26: PetscInt d;
28: u[0] = 0.0;
29: for (d = 0; d < dim; ++d) u[0] += x[d]/(ctx->L[d]);
30: return 0;
31: }
33: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
34: {
35: AppCtx *ctx = (AppCtx *) a_ctx;
36: PetscInt d;
38: u[0] = 1;
39: for (d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d])*PetscSqr(ctx->L[d]) - PetscPowRealInt(x[d], 4);
40: return 0;
41: }
43: static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
44: {
45: AppCtx *ctx = (AppCtx *) a_ctx;
47: u[0] = PetscSinScalar(2*PETSC_PI*ctx->k*x[0]/(ctx->L[0]));
48: return 0;
49: }
51: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
52: {
53: char fstring[PETSC_MAX_PATH_LEN] = "linear";
54: PetscBool flag;
58: options->particlesPerCell = 1;
59: options->k = 1;
60: options->particleRelDx = 1.e-20;
61: options->meshRelDx = 1.e-20;
62: options->momentTol = 100.*PETSC_MACHINE_EPSILON;
63: options->useBlockDiagPrec = PETSC_FALSE;
64: PetscStrcpy(options->meshFilename, "");
66: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
67: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex2.c", options->meshFilename, options->meshFilename, sizeof(options->meshFilename), NULL);
68: PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
69: PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
70: PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
71: PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
72: PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL);
73: PetscStrcmp(fstring, "linear", &flag);
74: if (flag) {
75: options->func = linear;
76: } else {
77: PetscStrcmp(fstring, "sin", &flag);
78: if (flag) {
79: options->func = sinx;
80: } else {
81: PetscStrcmp(fstring, "x2_x4", &flag);
82: options->func = x2_x4;
83: if (!flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown function %s",fstring);
84: }
85: }
86: PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL);
87: PetscOptionsBool("-block_diag_prec", "Use the block diagonal of the normal equations to precondition the particle projection", "ex2.c", options->useBlockDiagPrec, &options->useBlockDiagPrec, NULL);
88: PetscOptionsEnd();
90: return(0);
91: }
93: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
94: {
95: PetscRandom rnd;
96: PetscReal interval = user->meshRelDx;
97: Vec coordinates;
98: PetscScalar *coords;
99: PetscReal *hh, low[3], high[3];
100: PetscInt d, cdim, cEnd, N, p, bs;
104: DMGetBoundingBox(dm, low, high);
105: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
106: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
107: PetscRandomSetInterval(rnd, -interval, interval);
108: PetscRandomSetFromOptions(rnd);
109: DMGetCoordinatesLocal(dm, &coordinates);
110: DMGetCoordinateDim(dm, &cdim);
111: PetscCalloc1(cdim,&hh);
112: for (d = 0; d < cdim; ++d) hh[d] = (user->L[d])/PetscPowReal(cEnd, 1./cdim);
113: VecGetLocalSize(coordinates, &N);
114: VecGetBlockSize(coordinates, &bs);
115: if (bs != cdim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %D != %D", bs, cdim);
116: VecGetArray(coordinates, &coords);
117: for (p = 0; p < N; p += cdim) {
118: PetscScalar *coord = &coords[p], value;
120: for (d = 0; d < cdim; ++d) {
121: PetscRandomGetValue(rnd, &value);
122: coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value*hh[d])));
123: }
124: }
125: VecRestoreArray(coordinates, &coords);
126: PetscRandomDestroy(&rnd);
127: PetscFree(hh);
128: return(0);
129: }
131: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
132: {
133: PetscReal low[3], high[3];
134: PetscInt cdim, d;
135: PetscBool flg;
139: PetscStrcmp(user->meshFilename, "", &flg);
140: if (flg) {
141: DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
142: } else {
143: DMPlexCreateFromFile(comm, user->meshFilename, PETSC_TRUE, dm);
144: }
145: DMGetCoordinateDim(*dm, &cdim);
146: DMGetBoundingBox(*dm, low, high);
147: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
148: {
149: DM distributedMesh = NULL;
151: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
152: if (distributedMesh) {
153: DMDestroy(dm);
154: *dm = distributedMesh;
155: }
156: }
157: DMLocalizeCoordinates(*dm); /* needed for periodic */
158: DMSetFromOptions(*dm);
159: PerturbVertices(*dm, user);
160: PetscObjectSetName((PetscObject) *dm, "Mesh");
161: DMViewFromOptions(*dm, NULL, "-dm_view");
162: return(0);
163: }
165: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
166: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
167: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
168: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
169: {
170: g0[0] = 1.0;
171: }
173: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
174: {
175: PetscFE fe;
176: PetscDS ds;
177: DMPolytopeType ct;
178: PetscBool simplex;
179: PetscInt dim, cStart;
183: DMGetDimension(dm, &dim);
184: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
185: DMPlexGetCellType(dm, cStart, &ct);
186: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
187: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, NULL, -1, &fe);
188: PetscObjectSetName((PetscObject) fe, "fe");
189: DMSetField(dm, 0, NULL, (PetscObject) fe);
190: DMCreateDS(dm);
191: PetscFEDestroy(&fe);
192: /* Setup to form mass matrix */
193: DMGetDS(dm, &ds);
194: PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
195: return(0);
196: }
198: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
199: {
200: PetscRandom rnd, rndp;
201: PetscReal interval = user->particleRelDx;
202: DMPolytopeType ct;
203: PetscBool simplex;
204: PetscScalar value, *vals;
205: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
206: PetscInt *cellid;
207: PetscInt Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;
211: DMGetDimension(dm, &dim);
212: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
213: DMPlexGetCellType(dm, cStart, &ct);
214: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
215: DMCreate(PetscObjectComm((PetscObject) dm), sw);
216: DMSetType(*sw, DMSWARM);
217: DMSetDimension(*sw, dim);
219: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
220: PetscRandomSetInterval(rnd, -1.0, 1.0);
221: PetscRandomSetFromOptions(rnd);
222: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
223: PetscRandomSetInterval(rndp, -interval, interval);
224: PetscRandomSetFromOptions(rndp);
226: DMSwarmSetType(*sw, DMSWARM_PIC);
227: DMSwarmSetCellDM(*sw, dm);
228: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
229: DMSwarmFinalizeFieldRegister(*sw);
230: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
231: DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
232: DMSetFromOptions(*sw);
233: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
234: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
235: DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **) &vals);
237: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
238: for (c = 0; c < Ncell; ++c) {
239: if (Np == 1) {
240: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
241: cellid[c] = c;
242: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
243: } else {
244: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
245: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
246: for (p = 0; p < Np; ++p) {
247: const PetscInt n = c*Np + p;
248: PetscReal sum = 0.0, refcoords[3];
250: cellid[n] = c;
251: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rnd, &value); refcoords[d] = PetscRealPart(value); sum += refcoords[d];}
252: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
253: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
254: }
255: }
256: }
257: PetscFree5(centroid, xi0, v0, J, invJ);
258: for (c = 0; c < Ncell; ++c) {
259: for (p = 0; p < Np; ++p) {
260: const PetscInt n = c*Np + p;
262: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value); coords[n*dim+d] += PetscRealPart(value);}
263: user->func(dim, 0.0, &coords[n*dim], 1, &vals[c], user);
264: }
265: }
266: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
267: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
268: DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **) &vals);
269: PetscRandomDestroy(&rnd);
270: PetscRandomDestroy(&rndp);
271: PetscObjectSetName((PetscObject) *sw, "Particles");
272: DMViewFromOptions(*sw, NULL, "-sw_view");
273: return(0);
274: }
276: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
277: {
278: DM dm;
279: const PetscReal *coords;
280: const PetscScalar *w;
281: PetscReal mom[3] = {0.0, 0.0, 0.0};
282: PetscInt cell, cStart, cEnd, dim;
283: PetscErrorCode ierr;
286: DMGetDimension(sw, &dim);
287: DMSwarmGetCellDM(sw, &dm);
288: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
289: DMSwarmSortGetAccess(sw);
290: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
291: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **) &w);
292: for (cell = cStart; cell < cEnd; ++cell) {
293: PetscInt *pidx;
294: PetscInt Np, p, d;
296: DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
297: for (p = 0; p < Np; ++p) {
298: const PetscInt idx = pidx[p];
299: const PetscReal *c = &coords[idx*dim];
301: mom[0] += PetscRealPart(w[idx]);
302: mom[1] += PetscRealPart(w[idx]) * c[0];
303: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d]*c[d];
304: }
305: PetscFree(pidx);
306: }
307: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
308: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &w);
309: DMSwarmSortRestoreAccess(sw);
310: MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) sw));
311: return(0);
312: }
314: static void f0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux,
315: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
316: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
317: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
318: {
319: f0[0] = u[0];
320: }
322: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
323: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
324: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
325: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
326: {
327: f0[0] = x[0]*u[0];
328: }
330: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
331: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
332: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
333: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
334: {
335: PetscInt d;
337: f0[0] = 0.0;
338: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
339: }
341: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
342: {
343: PetscDS prob;
345: PetscScalar mom;
348: DMGetDS(dm, &prob);
349: PetscDSSetObjective(prob, 0, &f0_1);
350: DMPlexComputeIntegralFEM(dm, u, &mom, user);
351: moments[0] = PetscRealPart(mom);
352: PetscDSSetObjective(prob, 0, &f0_x);
353: DMPlexComputeIntegralFEM(dm, u, &mom, user);
354: moments[1] = PetscRealPart(mom);
355: PetscDSSetObjective(prob, 0, &f0_r2);
356: DMPlexComputeIntegralFEM(dm, u, &mom, user);
357: moments[2] = PetscRealPart(mom);
358: return(0);
359: }
361: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
362: {
363: MPI_Comm comm;
364: KSP ksp;
365: Mat M; /* FEM mass matrix */
366: Mat M_p; /* Particle mass matrix */
367: Vec f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
368: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
369: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
370: PetscInt m;
374: PetscObjectGetComm((PetscObject) dm, &comm);
375: KSPCreate(comm, &ksp);
376: KSPSetOptionsPrefix(ksp, "ptof_");
377: DMGetGlobalVector(dm, &fhat);
378: DMGetGlobalVector(dm, &rhs);
380: DMCreateMassMatrix(sw, dm, &M_p);
381: MatViewFromOptions(M_p, NULL, "-M_p_view");
383: /* make particle weight vector */
384: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
386: /* create matrix RHS vector */
387: MatMultTranspose(M_p, f, rhs);
388: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
389: PetscObjectSetName((PetscObject) rhs,"rhs");
390: VecViewFromOptions(rhs, NULL, "-rhs_view");
392: DMCreateMatrix(dm, &M);
393: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
394: MatViewFromOptions(M, NULL, "-M_view");
395: KSPSetOperators(ksp, M, M);
396: KSPSetFromOptions(ksp);
397: KSPSolve(ksp, rhs, fhat);
398: PetscObjectSetName((PetscObject) fhat,"fhat");
399: VecViewFromOptions(fhat, NULL, "-fhat_view");
401: /* Check moments of field */
402: computeParticleMoments(sw, pmoments, user);
403: computeFEMMoments(dm, fhat, fmoments, user);
404: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
405: for (m = 0; m < 3; ++m) {
406: 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);
407: }
409: KSPDestroy(&ksp);
410: MatDestroy(&M);
411: MatDestroy(&M_p);
412: DMRestoreGlobalVector(dm, &fhat);
413: DMRestoreGlobalVector(dm, &rhs);
415: return(0);
416: }
419: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
420: {
422: MPI_Comm comm;
423: KSP ksp;
424: Mat M; /* FEM mass matrix */
425: Mat M_p, PM_p; /* Particle mass matrix M_p, and the preconditioning matrix, e.g. M_p M^T_p */
426: Vec f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
427: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
428: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
429: PetscInt m;
433: PetscObjectGetComm((PetscObject) dm, &comm);
435: KSPCreate(comm, &ksp);
436: KSPSetOptionsPrefix(ksp, "ftop_");
437: KSPSetFromOptions(ksp);
439: DMGetGlobalVector(dm, &fhat);
440: DMGetGlobalVector(dm, &rhs);
442: DMCreateMassMatrix(sw, dm, &M_p);
443: MatViewFromOptions(M_p, NULL, "-M_p_view");
445: /* make particle weight vector */
446: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
448: /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
449: PetscObjectSetName((PetscObject) rhs,"rhs");
450: VecViewFromOptions(rhs, NULL, "-rhs_view");
451: DMCreateMatrix(dm, &M);
452: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
453: MatViewFromOptions(M, NULL, "-M_view");
454: MatMultTranspose(M, fhat, rhs);
455: if (user->useBlockDiagPrec) {DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);}
456: else {PetscObjectReference((PetscObject) M_p); PM_p = M_p;}
458: KSPSetOperators(ksp, M_p, PM_p);
459: KSPSolveTranspose(ksp, rhs, f);
460: PetscObjectSetName((PetscObject) fhat,"fhat");
461: VecViewFromOptions(fhat, NULL, "-fhat_view");
463: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
465: /* Check moments */
466: computeParticleMoments(sw, pmoments, user);
467: computeFEMMoments(dm, fhat, fmoments, user);
468: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
469: for (m = 0; m < 3; ++m) {
470: 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);
471: }
472: KSPDestroy(&ksp);
473: MatDestroy(&M);
474: MatDestroy(&M_p);
475: MatDestroy(&PM_p);
476: DMRestoreGlobalVector(dm, &fhat);
477: DMRestoreGlobalVector(dm, &rhs);
478: return(0);
479: }
481: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
482: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC){
484: DM_Plex *mesh = (DM_Plex *) dm->data;
485: PetscInt debug = mesh->printFEM;
486: DM dmC;
487: PetscSection section;
488: PetscQuadrature quad = NULL;
489: PetscScalar *interpolant, *gradsum;
490: PetscFEGeom fegeom;
491: PetscReal *coords;
492: const PetscReal *quadPoints, *quadWeights;
493: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
494: PetscErrorCode ierr;
497: VecGetDM(locC, &dmC);
498: VecSet(locC, 0.0);
499: DMGetDimension(dm, &dim);
500: DMGetCoordinateDim(dm, &coordDim);
501: fegeom.dimEmbed = coordDim;
502: DMGetLocalSection(dm, §ion);
503: PetscSectionGetNumFields(section, &numFields);
504: for (field = 0; field < numFields; ++field) {
505: PetscObject obj;
506: PetscClassId id;
507: PetscInt Nc;
509: DMGetField(dm, field, NULL, &obj);
510: PetscObjectGetClassId(obj, &id);
511: if (id == PETSCFE_CLASSID) {
512: PetscFE fe = (PetscFE) obj;
514: PetscFEGetQuadrature(fe, &quad);
515: PetscFEGetNumComponents(fe, &Nc);
516: } else if (id == PETSCFV_CLASSID) {
517: PetscFV fv = (PetscFV) obj;
519: PetscFVGetQuadrature(fv, &quad);
520: PetscFVGetNumComponents(fv, &Nc);
521: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
522: numComponents += Nc;
523: }
524: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
525: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
526: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
527: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
528: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
529: for (v = vStart; v < vEnd; ++v) {
530: PetscInt *star = NULL;
531: PetscInt starSize, st, d, fc;
533: PetscArrayzero(gradsum, coordDim*numComponents);
534: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
535: for (st = 0; st < starSize*2; st += 2) {
536: const PetscInt cell = star[st];
537: PetscScalar *grad = &gradsum[coordDim*numComponents];
538: PetscScalar *x = NULL;
540: if ((cell < cStart) || (cell >= cEnd)) continue;
541: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
542: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
543: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
544: PetscObject obj;
545: PetscClassId id;
546: PetscInt Nb, Nc, q, qc = 0;
548: PetscArrayzero(grad, coordDim*numComponents);
549: DMGetField(dm, field, NULL, &obj);
550: PetscObjectGetClassId(obj, &id);
551: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
552: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
553: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
554: for (q = 0; q < Nq; ++q) {
555: 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);
556: if (ierr) {
557: PetscErrorCode ierr2;
558: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
559: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
560: ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
561:
562: }
563: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &fegeom, q, interpolant);}
564: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
565: for (fc = 0; fc < Nc; ++fc) {
566: const PetscReal wt = quadWeights[q*qNc+qc+fc];
568: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
569: }
570: }
571: fieldOffset += Nb;
572: }
573: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
574: for (fc = 0; fc < numComponents; ++fc) {
575: for (d = 0; d < coordDim; ++d) {
576: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
577: }
578: }
579: if (debug) {
580: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
581: for (fc = 0; fc < numComponents; ++fc) {
582: for (d = 0; d < coordDim; ++d) {
583: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
584: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
585: }
586: }
587: PetscPrintf(PETSC_COMM_SELF, "]\n");
588: }
589: }
590: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
591: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
592: }
593: PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
594: return(0);
595: }
597: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
598: {
600: MPI_Comm comm;
601: KSP ksp;
602: Mat M; /* FEM mass matrix */
603: Mat M_p; /* Particle mass matrix */
604: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
605: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
606: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
607: PetscInt m;
611: PetscObjectGetComm((PetscObject) dm, &comm);
612: KSPCreate(comm, &ksp);
613: KSPSetOptionsPrefix(ksp, "ptof_");
614: DMGetGlobalVector(dm, &fhat);
615: DMGetGlobalVector(dm, &rhs);
616: DMGetGlobalVector(dm, &grad);
618: DMCreateMassMatrix(sw, dm, &M_p);
619: MatViewFromOptions(M_p, NULL, "-M_p_view");
621: /* make particle weight vector */
622: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
624: /* create matrix RHS vector */
625: MatMultTranspose(M_p, f, rhs);
626: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
627: PetscObjectSetName((PetscObject) rhs,"rhs");
628: VecViewFromOptions(rhs, NULL, "-rhs_view");
630: DMCreateMatrix(dm, &M);
631: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
633: InterpolateGradient(dm, fhat, grad);
635: MatViewFromOptions(M, NULL, "-M_view");
636: KSPSetOperators(ksp, M, M);
637: KSPSetFromOptions(ksp);
638: KSPSolve(ksp, rhs, grad);
639: PetscObjectSetName((PetscObject) fhat,"fhat");
640: VecViewFromOptions(fhat, NULL, "-fhat_view");
642: /* Check moments of field */
643: computeParticleMoments(sw, pmoments, user);
644: computeFEMMoments(dm, grad, fmoments, user);
645: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
646: for (m = 0; m < 3; ++m) {
647: 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);
648: }
650: KSPDestroy(&ksp);
651: MatDestroy(&M);
652: MatDestroy(&M_p);
653: DMRestoreGlobalVector(dm, &fhat);
654: DMRestoreGlobalVector(dm, &rhs);
655: DMRestoreGlobalVector(dm, &grad);
657: return(0);
658: }
660: int main (int argc, char * argv[]) {
661: MPI_Comm comm;
662: DM dm, sw;
663: AppCtx user;
666: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
667: comm = PETSC_COMM_WORLD;
668: ProcessOptions(comm, &user);
669: CreateMesh(comm, &dm, &user);
670: CreateFEM(dm, &user);
671: CreateParticles(dm, &sw, &user);
672: TestL2ProjectionParticlesToField(dm, sw, &user);
673: TestL2ProjectionFieldToParticles(dm, sw, &user);
674: TestFieldGradientProjection(dm, sw, &user);
675: DMDestroy(&dm);
676: DMDestroy(&sw);
677: PetscFinalize();
678: return ierr;
679: }
682: /*TEST
684: # Swarm does not handle complex or quad
685: build:
686: requires: !complex double
688: test:
689: suffix: proj_tri_0
690: requires: triangle
691: 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
692: filter: grep -v marker | grep -v atomic | grep -v usage
694: test:
695: suffix: proj_tri_2_faces
696: requires: triangle
697: 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
698: filter: grep -v marker | grep -v atomic | grep -v usage
700: test:
701: suffix: proj_quad_0
702: requires: triangle
703: args: -dm_plex_box_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
704: filter: grep -v marker | grep -v atomic | grep -v usage
706: test:
707: suffix: proj_quad_2_faces
708: requires: triangle
709: args: -dm_plex_box_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
710: filter: grep -v marker | grep -v atomic | grep -v usage
712: test:
713: suffix: proj_tri_5P
714: requires: triangle
715: 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
716: filter: grep -v marker | grep -v atomic | grep -v usage
718: test:
719: suffix: proj_quad_5P
720: requires: triangle
721: args: -dm_plex_box_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
722: filter: grep -v marker | grep -v atomic | grep -v usage
724: test:
725: suffix: proj_tri_mdx
726: requires: triangle
727: 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
728: filter: grep -v marker | grep -v atomic | grep -v usage
730: test:
731: suffix: proj_tri_mdx_5P
732: requires: triangle
733: 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
734: filter: grep -v marker | grep -v atomic | grep -v usage
736: test:
737: suffix: proj_tri_3d
738: requires: ctetgen
739: args: -dm_plex_box_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
740: filter: grep -v marker | grep -v atomic | grep -v usage
742: test:
743: suffix: proj_tri_3d_2_faces
744: requires: ctetgen
745: args: -dm_plex_box_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
746: filter: grep -v marker | grep -v atomic | grep -v usage
748: test:
749: suffix: proj_tri_3d_5P
750: requires: ctetgen
751: args: -dm_plex_box_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
752: filter: grep -v marker | grep -v atomic | grep -v usage
754: test:
755: suffix: proj_tri_3d_mdx
756: requires: ctetgen
757: args: -dm_plex_box_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
758: filter: grep -v marker | grep -v atomic | grep -v usage
760: test:
761: suffix: proj_tri_3d_mdx_5P
762: requires: ctetgen
763: args: -dm_plex_box_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
764: filter: grep -v marker | grep -v atomic | grep -v usage
766: test:
767: suffix: proj_tri_3d_mdx_2_faces
768: requires: ctetgen
769: args: -dm_plex_box_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
770: filter: grep -v marker | grep -v atomic | grep -v usage
772: test:
773: suffix: proj_tri_3d_mdx_5P_2_faces
774: requires: ctetgen
775: args: -dm_plex_box_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
776: filter: grep -v marker | grep -v atomic | grep -v usage
778: test:
779: suffix: proj_quad_lsqr_scale
780: requires: !complex
781: args: -dm_plex_box_simplex 0 -dm_plex_box_faces 4,4 \
782: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
783: -particlesPerCell 200 \
784: -ptof_pc_type lu \
785: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none
787: test:
788: suffix: proj_quad_lsqr_prec_scale
789: requires: !complex
790: args: -dm_plex_box_simplex 0 -dm_plex_box_faces 4,4 \
791: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
792: -particlesPerCell 200 \
793: -ptof_pc_type lu \
794: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type lu -ftop_pc_factor_shift_type nonzero -block_diag_prec
796: TEST*/