Actual source code: ex2.c
petsc-3.13.4 2020-08-01
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: PetscInt dim; /* The topological mesh dimension */
11: PetscBool simplex; /* Flag for simplices or tensor cells */
12: char meshFilename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
13: PetscInt faces; /* Number of faces per edge if unit square/cube generated */
14: PetscReal domain_lo[3], domain_hi[3]; /* Lower left and upper right mesh corners */
15: DMBoundaryType boundary[3]; /* The domain boundary type, e.g. periodic */
16: PetscInt particlesPerCell; /* The number of partices per cell */
17: PetscReal particleRelDx; /* Relative particle position perturbation compared to average cell diameter h */
18: PetscReal meshRelDx; /* Relative vertex position perturbation compared to average cell diameter h */
19: PetscInt k; /* Mode number for test function */
20: PetscReal momentTol; /* Tolerance for checking moment conservation */
21: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
22: } AppCtx;
24: /* const char *const ex2FunctionTypes[] = {"linear","x2_x4","sin","ex2FunctionTypes","EX2_FUNCTION_",0}; */
26: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
27: {
28: AppCtx *ctx = (AppCtx *) a_ctx;
29: PetscInt d;
31: u[0] = 0.0;
32: for (d = 0; d < dim; ++d) u[0] += x[d]/(ctx->domain_hi[d] - ctx->domain_lo[d]);
33: return 0;
34: }
36: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
37: {
38: AppCtx *ctx = (AppCtx *) a_ctx;
39: PetscInt d;
41: u[0] = 1;
42: for (d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d])*PetscSqr(ctx->domain_hi[d]) - PetscPowRealInt(x[d], 4);
43: return 0;
44: }
46: static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
47: {
48: AppCtx *ctx = (AppCtx *) a_ctx;
50: u[0] = PetscSinScalar(2*PETSC_PI*ctx->k*x[0]/(ctx->domain_hi[0] - ctx->domain_lo[0]));
51: return 0;
52: }
54: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
55: {
56: PetscInt ii, bd;
57: char fstring[PETSC_MAX_PATH_LEN] = "linear";
58: PetscBool flag;
62: options->dim = 2;
63: options->simplex = PETSC_TRUE;
64: options->faces = 1;
65: options->domain_lo[0] = 0.0;
66: options->domain_lo[1] = 0.0;
67: options->domain_lo[2] = 0.0;
68: options->domain_hi[0] = 1.0;
69: options->domain_hi[1] = 1.0;
70: options->domain_hi[2] = 1.0;
71: options->boundary[0] = DM_BOUNDARY_NONE; /* PERIODIC (plotting does not work in parallel, moments not conserved) */
72: options->boundary[1] = DM_BOUNDARY_NONE;
73: options->boundary[2] = DM_BOUNDARY_NONE;
74: options->particlesPerCell = 1;
75: options->k = 1;
76: options->particleRelDx = 1.e-20;
77: options->meshRelDx = 1.e-20;
78: options->momentTol = 100.*PETSC_MACHINE_EPSILON;
79: PetscStrcpy(options->meshFilename, "");
81: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
82: PetscOptionsInt("-dim", "The topological mesh dimension", "ex2.c", options->dim, &options->dim, NULL);
83: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex2.c", options->simplex, &options->simplex, NULL);
84: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex2.c", options->meshFilename, options->meshFilename, PETSC_MAX_PATH_LEN, NULL);
85: PetscOptionsInt("-faces", "Number of faces per edge if unit square/cube generated", "ex2.c", options->faces, &options->faces, NULL);
86: PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL);
87: PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL);
88: PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL);
89: PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL);
90: ii = options->dim;
91: PetscOptionsRealArray("-domain_hi", "Domain size", "ex2.c", options->domain_hi, &ii, NULL);
92: ii = options->dim;
93: PetscOptionsRealArray("-domain_lo", "Domain size", "ex2.c", options->domain_lo, &ii, NULL);
94: bd = options->boundary[0];
95: PetscOptionsEList("-x_boundary", "The x-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[0]], &bd, NULL);
96: options->boundary[0] = (DMBoundaryType) bd;
97: bd = options->boundary[1];
98: PetscOptionsEList("-y_boundary", "The y-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[1]], &bd, NULL);
99: options->boundary[1] = (DMBoundaryType) bd;
100: bd = options->boundary[2];
101: PetscOptionsEList("-z_boundary", "The z-boundary", "ex2.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->boundary[2]], &bd, NULL);
102: options->boundary[2] = (DMBoundaryType) bd;
103: PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, PETSC_MAX_PATH_LEN, NULL);
104: PetscStrcmp(fstring, "linear", &flag);
105: if (flag) {
106: options->func = linear;
107: } else {
108: PetscStrcmp(fstring, "sin", &flag);
109: if (flag) {
110: options->func = sinx;
111: } else {
112: PetscStrcmp(fstring, "x2_x4", &flag);
113: options->func = x2_x4;
114: if (!flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown function %s",fstring);
115: }
116: }
117: PetscOptionsEnd();
119: return(0);
120: }
122: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
123: {
124: PetscRandom rnd;
125: PetscReal interval = user->meshRelDx;
126: Vec coordinates;
127: PetscScalar *coords;
128: PetscReal *hh;
129: PetscInt d, cdim, N, p, bs;
133: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
134: PetscRandomSetInterval(rnd, -interval, interval);
135: PetscRandomSetFromOptions(rnd);
136: DMGetCoordinatesLocal(dm, &coordinates);
137: DMGetCoordinateDim(dm, &cdim);
138: PetscCalloc1(PetscMax(user->dim,cdim),&hh);
139: for (d = 0; d < user->dim; ++d) hh[d] = (user->domain_hi[d] - user->domain_lo[d])/user->faces;
140: VecGetLocalSize(coordinates, &N);
141: VecGetBlockSize(coordinates, &bs);
142: if (bs != cdim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %D != %D", bs, cdim);
143: VecGetArray(coordinates, &coords);
144: for (p = 0; p < N; p += cdim) {
145: PetscScalar *coord = &coords[p], value;
147: for (d = 0; d < cdim; ++d) {
148: PetscRandomGetValue(rnd, &value);
149: coord[d] = PetscMax(user->domain_lo[d], PetscMin(user->domain_hi[d], PetscRealPart(coord[d] + value*hh[d])));
150: }
151: }
152: VecRestoreArray(coordinates, &coords);
153: PetscRandomDestroy(&rnd);
154: PetscFree(hh);
155: return(0);
156: }
158: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
159: {
160: PetscBool flg;
164: PetscStrcmp(user->meshFilename, "", &flg);
165: if (flg) {
166: PetscInt faces[3];
168: faces[0] = user->faces; faces[1] = user->faces; faces[2] = user->faces;
169: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, faces, user->domain_lo, user->domain_hi, user->boundary, PETSC_TRUE, dm);
170: } else {
171: DMPlexCreateFromFile(comm, user->meshFilename, PETSC_TRUE, dm);
172: DMGetDimension(*dm, &user->dim);
173: }
174: {
175: DM distributedMesh = NULL;
177: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
178: if (distributedMesh) {
179: DMDestroy(dm);
180: *dm = distributedMesh;
181: }
182: }
183: DMLocalizeCoordinates(*dm); /* needed for periodic */
184: DMSetFromOptions(*dm);
185: PerturbVertices(*dm, user);
186: PetscObjectSetName((PetscObject) *dm, "Mesh");
187: DMViewFromOptions(*dm, NULL, "-dm_view");
188: return(0);
189: }
191: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
192: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
193: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
194: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
195: {
196: g0[0] = 1.0;
197: }
199: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
200: {
201: PetscFE fe;
202: PetscDS ds;
203: PetscInt dim;
207: DMGetDimension(dm, &dim);
208: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, user->simplex, NULL, -1, &fe);
209: PetscObjectSetName((PetscObject) fe, "fe");
210: DMSetField(dm, 0, NULL, (PetscObject) fe);
211: DMCreateDS(dm);
212: PetscFEDestroy(&fe);
213: /* Setup to form mass matrix */
214: DMGetDS(dm, &ds);
215: PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL);
216: return(0);
217: }
219: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
220: {
221: PetscRandom rnd, rndp;
222: PetscReal interval = user->particleRelDx;
223: PetscScalar value, *vals;
224: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
225: PetscInt *cellid;
226: PetscInt Ncell, Np = user->particlesPerCell, p, c, dim, d;
230: DMGetDimension(dm, &dim);
231: DMCreate(PetscObjectComm((PetscObject) dm), sw);
232: DMSetType(*sw, DMSWARM);
233: DMSetDimension(*sw, dim);
235: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
236: PetscRandomSetInterval(rnd, -1.0, 1.0);
237: PetscRandomSetFromOptions(rnd);
238: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rndp);
239: PetscRandomSetInterval(rndp, -interval, interval);
240: PetscRandomSetFromOptions(rndp);
242: DMSwarmSetType(*sw, DMSWARM_PIC);
243: DMSwarmSetCellDM(*sw, dm);
244: DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
245: DMSwarmFinalizeFieldRegister(*sw);
246: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
247: DMSwarmSetLocalSizes(*sw, Ncell * Np, 0);
248: DMSetFromOptions(*sw);
249: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
250: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
251: DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **) &vals);
253: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
254: for (c = 0; c < Ncell; ++c) {
255: if (Np == 1) {
256: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
257: cellid[c] = c;
258: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
259: } else {
260: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
261: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
262: for (p = 0; p < Np; ++p) {
263: const PetscInt n = c*Np + p;
264: PetscReal sum = 0.0, refcoords[3];
266: cellid[n] = c;
267: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rnd, &value); refcoords[d] = PetscRealPart(value); sum += refcoords[d];}
268: if (user->simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
269: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
270: }
271: }
272: }
273: PetscFree5(centroid, xi0, v0, J, invJ);
274: for (c = 0; c < Ncell; ++c) {
275: for (p = 0; p < Np; ++p) {
276: const PetscInt n = c*Np + p;
278: for (d = 0; d < dim; ++d) {PetscRandomGetValue(rndp, &value); coords[n*dim+d] += PetscRealPart(value);}
279: user->func(dim, 0.0, &coords[n*dim], 1, &vals[c], user);
280: }
281: }
282: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
283: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
284: DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **) &vals);
285: PetscRandomDestroy(&rnd);
286: PetscRandomDestroy(&rndp);
287: PetscObjectSetName((PetscObject) *sw, "Particles");
288: DMViewFromOptions(*sw, NULL, "-sw_view");
289: return(0);
290: }
292: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
293: {
294: DM dm;
295: const PetscReal *coords;
296: const PetscScalar *w;
297: PetscReal mom[3] = {0.0, 0.0, 0.0};
298: PetscInt cell, cStart, cEnd, dim;
299: PetscErrorCode ierr;
302: DMGetDimension(sw, &dim);
303: DMSwarmGetCellDM(sw, &dm);
304: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
305: DMSwarmSortGetAccess(sw);
306: DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
307: DMSwarmGetField(sw, "w_q", NULL, NULL, (void **) &w);
308: for (cell = cStart; cell < cEnd; ++cell) {
309: PetscInt *pidx;
310: PetscInt Np, p, d;
312: DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx);
313: for (p = 0; p < Np; ++p) {
314: const PetscInt idx = pidx[p];
315: const PetscReal *c = &coords[idx*dim];
317: mom[0] += PetscRealPart(w[idx]);
318: mom[1] += PetscRealPart(w[idx]) * c[0];
319: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d]*c[d];
320: }
321: PetscFree(pidx);
322: }
323: DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
324: DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &w);
325: DMSwarmSortRestoreAccess(sw);
326: MPI_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) sw));
327: return(0);
328: }
330: static void f0_1(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: f0[0] = u[0];
336: }
338: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux,
339: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
340: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
341: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
342: {
343: f0[0] = x[0]*u[0];
344: }
346: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
347: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
348: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
349: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
350: {
351: PetscInt d;
353: f0[0] = 0.0;
354: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d])*u[0];
355: }
357: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
358: {
359: PetscDS prob;
361: PetscScalar mom;
364: DMGetDS(dm, &prob);
365: PetscDSSetObjective(prob, 0, &f0_1);
366: DMPlexComputeIntegralFEM(dm, u, &mom, user);
367: moments[0] = PetscRealPart(mom);
368: PetscDSSetObjective(prob, 0, &f0_x);
369: DMPlexComputeIntegralFEM(dm, u, &mom, user);
370: moments[1] = PetscRealPart(mom);
371: PetscDSSetObjective(prob, 0, &f0_r2);
372: DMPlexComputeIntegralFEM(dm, u, &mom, user);
373: moments[2] = PetscRealPart(mom);
374: return(0);
375: }
377: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, AppCtx *user)
378: {
379: MPI_Comm comm;
380: KSP ksp;
381: Mat M; /* FEM mass matrix */
382: Mat M_p; /* Particle mass matrix */
383: Vec f, rhs, fhat; /* Particle field f, \int phi_i f, FEM field */
384: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
385: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
386: PetscInt m;
390: PetscObjectGetComm((PetscObject) dm, &comm);
391: KSPCreate(comm, &ksp);
392: KSPSetOptionsPrefix(ksp, "ptof_");
393: DMGetGlobalVector(dm, &fhat);
394: DMGetGlobalVector(dm, &rhs);
396: DMCreateMassMatrix(sw, dm, &M_p);
397: MatViewFromOptions(M_p, NULL, "-M_p_view");
399: /* make particle weight vector */
400: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
402: /* create matrix RHS vector */
403: MatMultTranspose(M_p, f, rhs);
404: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
405: PetscObjectSetName((PetscObject) rhs,"rhs");
406: VecViewFromOptions(rhs, NULL, "-rhs_view");
408: DMCreateMatrix(dm, &M);
409: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
410: MatViewFromOptions(M, NULL, "-M_view");
411: KSPSetOperators(ksp, M, M);
412: KSPSetFromOptions(ksp);
413: KSPSolve(ksp, rhs, fhat);
414: PetscObjectSetName((PetscObject) fhat,"fhat");
415: VecViewFromOptions(fhat, NULL, "-fhat_view");
417: /* Check moments of field */
418: computeParticleMoments(sw, pmoments, user);
419: computeFEMMoments(dm, fhat, fmoments, user);
420: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
421: for (m = 0; m < 3; ++m) {
422: 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);
423: }
425: KSPDestroy(&ksp);
426: MatDestroy(&M);
427: MatDestroy(&M_p);
428: DMRestoreGlobalVector(dm, &fhat);
429: DMRestoreGlobalVector(dm, &rhs);
431: return(0);
432: }
435: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, AppCtx *user)
436: {
438: MPI_Comm comm;
439: KSP ksp;
440: Mat M; /* FEM mass matrix */
441: Mat M_p; /* Particle mass matrix */
442: Vec f, rhs, fhat; /* Particle field f, \int phi_i fhat, FEM field */
443: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
444: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
445: PetscInt m;
449: PetscObjectGetComm((PetscObject) dm, &comm);
451: KSPCreate(comm, &ksp);
452: KSPSetOptionsPrefix(ksp, "ftop_");
453: KSPSetFromOptions(ksp);
455: DMGetGlobalVector(dm, &fhat);
456: DMGetGlobalVector(dm, &rhs);
458: DMCreateMassMatrix(sw, dm, &M_p);
459: MatViewFromOptions(M_p, NULL, "-M_p_view");
461: /* make particle weight vector */
462: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
464: /* create matrix RHS vector, in this case the FEM field fhat with the coefficients vector #alpha */
465: PetscObjectSetName((PetscObject) rhs,"rhs");
466: VecViewFromOptions(rhs, NULL, "-rhs_view");
467: DMCreateMatrix(dm, &M);
468: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
469: MatViewFromOptions(M, NULL, "-M_view");
470: MatMultTranspose(M, fhat, rhs);
472: KSPSetOperators(ksp, M_p, M_p);
473: KSPSolveTranspose(ksp, rhs, f);
474: PetscObjectSetName((PetscObject) fhat,"fhat");
475: VecViewFromOptions(fhat, NULL, "-fhat_view");
477: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
479: /* Check moments */
480: computeParticleMoments(sw, pmoments, user);
481: computeFEMMoments(dm, fhat, fmoments, user);
482: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
483: for (m = 0; m < 3; ++m) {
484: 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);
485: }
486: KSPDestroy(&ksp);
487: MatDestroy(&M);
488: MatDestroy(&M_p);
489: DMRestoreGlobalVector(dm, &fhat);
490: DMRestoreGlobalVector(dm, &rhs);
491: return(0);
492: }
494: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
495: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC){
497: DM_Plex *mesh = (DM_Plex *) dm->data;
498: PetscInt debug = mesh->printFEM;
499: DM dmC;
500: PetscSection section;
501: PetscQuadrature quad = NULL;
502: PetscScalar *interpolant, *gradsum;
503: PetscFEGeom fegeom;
504: PetscReal *coords;
505: const PetscReal *quadPoints, *quadWeights;
506: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
507: PetscErrorCode ierr;
510: VecGetDM(locC, &dmC);
511: VecSet(locC, 0.0);
512: DMGetDimension(dm, &dim);
513: DMGetCoordinateDim(dm, &coordDim);
514: fegeom.dimEmbed = coordDim;
515: DMGetLocalSection(dm, §ion);
516: PetscSectionGetNumFields(section, &numFields);
517: for (field = 0; field < numFields; ++field) {
518: PetscObject obj;
519: PetscClassId id;
520: PetscInt Nc;
522: DMGetField(dm, field, NULL, &obj);
523: PetscObjectGetClassId(obj, &id);
524: if (id == PETSCFE_CLASSID) {
525: PetscFE fe = (PetscFE) obj;
527: PetscFEGetQuadrature(fe, &quad);
528: PetscFEGetNumComponents(fe, &Nc);
529: } else if (id == PETSCFV_CLASSID) {
530: PetscFV fv = (PetscFV) obj;
532: PetscFVGetQuadrature(fv, &quad);
533: PetscFVGetNumComponents(fv, &Nc);
534: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
535: numComponents += Nc;
536: }
537: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
538: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
539: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
540: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
541: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
542: for (v = vStart; v < vEnd; ++v) {
543: PetscInt *star = NULL;
544: PetscInt starSize, st, d, fc;
546: PetscArrayzero(gradsum, coordDim*numComponents);
547: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
548: for (st = 0; st < starSize*2; st += 2) {
549: const PetscInt cell = star[st];
550: PetscScalar *grad = &gradsum[coordDim*numComponents];
551: PetscScalar *x = NULL;
553: if ((cell < cStart) || (cell >= cEnd)) continue;
554: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
555: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
556: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
557: PetscObject obj;
558: PetscClassId id;
559: PetscInt Nb, Nc, q, qc = 0;
561: PetscArrayzero(grad, coordDim*numComponents);
562: DMGetField(dm, field, NULL, &obj);
563: PetscObjectGetClassId(obj, &id);
564: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
565: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
566: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
567: for (q = 0; q < Nq; ++q) {
568: 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);
569: if (ierr) {
570: PetscErrorCode ierr2;
571: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
572: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
573: ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
574:
575: }
576: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &fegeom, q, interpolant);}
577: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
578: for (fc = 0; fc < Nc; ++fc) {
579: const PetscReal wt = quadWeights[q*qNc+qc+fc];
581: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
582: }
583: }
584: fieldOffset += Nb;
585: }
586: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
587: for (fc = 0; fc < numComponents; ++fc) {
588: for (d = 0; d < coordDim; ++d) {
589: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
590: }
591: }
592: if (debug) {
593: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
594: for (fc = 0; fc < numComponents; ++fc) {
595: for (d = 0; d < coordDim; ++d) {
596: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
597: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
598: }
599: }
600: PetscPrintf(PETSC_COMM_SELF, "]\n");
601: }
602: }
603: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
604: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
605: }
606: PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
607: return(0);
608: }
610: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
611: {
613: MPI_Comm comm;
614: KSP ksp;
615: Mat M; /* FEM mass matrix */
616: Mat M_p; /* Particle mass matrix */
617: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
618: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
619: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
620: PetscInt m;
624: PetscObjectGetComm((PetscObject) dm, &comm);
625: KSPCreate(comm, &ksp);
626: KSPSetOptionsPrefix(ksp, "ptof_");
627: DMGetGlobalVector(dm, &fhat);
628: DMGetGlobalVector(dm, &rhs);
629: DMGetGlobalVector(dm, &grad);
631: DMCreateMassMatrix(sw, dm, &M_p);
632: MatViewFromOptions(M_p, NULL, "-M_p_view");
634: /* make particle weight vector */
635: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
637: /* create matrix RHS vector */
638: MatMultTranspose(M_p, f, rhs);
639: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
640: PetscObjectSetName((PetscObject) rhs,"rhs");
641: VecViewFromOptions(rhs, NULL, "-rhs_view");
643: DMCreateMatrix(dm, &M);
644: DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user);
646: InterpolateGradient(dm, fhat, grad);
648: MatViewFromOptions(M, NULL, "-M_view");
649: KSPSetOperators(ksp, M, M);
650: KSPSetFromOptions(ksp);
651: KSPSolve(ksp, rhs, grad);
652: PetscObjectSetName((PetscObject) fhat,"fhat");
653: VecViewFromOptions(fhat, NULL, "-fhat_view");
655: /* Check moments of field */
656: computeParticleMoments(sw, pmoments, user);
657: computeFEMMoments(dm, grad, fmoments, user);
658: PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]);
659: for (m = 0; m < 3; ++m) {
660: 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);
661: }
663: KSPDestroy(&ksp);
664: MatDestroy(&M);
665: MatDestroy(&M_p);
666: DMRestoreGlobalVector(dm, &fhat);
667: DMRestoreGlobalVector(dm, &rhs);
668: DMRestoreGlobalVector(dm, &grad);
670: return(0);
671: }
673: int main (int argc, char * argv[]) {
674: MPI_Comm comm;
675: DM dm, sw;
676: AppCtx user;
679: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
680: comm = PETSC_COMM_WORLD;
681: ProcessOptions(comm, &user);
682: CreateMesh(comm, &dm, &user);
683: CreateFEM(dm, &user);
684: CreateParticles(dm, &sw, &user);
685: TestL2ProjectionParticlesToField(dm, sw, &user);
686: TestL2ProjectionFieldToParticles(dm, sw, &user);
687: TestFieldGradientProjection(dm, sw, &user);
688: DMDestroy(&dm);
689: DMDestroy(&sw);
690: PetscFinalize();
691: return ierr;
692: }
695: /*TEST
697: test:
698: suffix: proj_tri_0
699: requires: triangle !complex
700: args: -dim 2 -faces 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_2_faces
705: requires: triangle !complex
706: args: -dim 2 -faces 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_quad_0
711: requires: triangle !complex
712: args: -dim 2 -simplex 0 -faces 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
713: filter: grep -v marker | grep -v atomic | grep -v usage
715: test:
716: suffix: proj_quad_2_faces
717: requires: triangle !complex
718: args: -dim 2 -simplex 0 -faces 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
719: filter: grep -v marker | grep -v atomic | grep -v usage
721: test:
722: suffix: proj_tri_5P
723: requires: triangle !complex
724: args: -dim 2 -faces 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
725: filter: grep -v marker | grep -v atomic | grep -v usage
727: test:
728: suffix: proj_quad_5P
729: requires: triangle !complex
730: args: -dim 2 -simplex 0 -faces 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
731: filter: grep -v marker | grep -v atomic | grep -v usage
733: test:
734: suffix: proj_tri_mdx
735: requires: triangle !complex
736: args: -dim 2 -faces 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
737: filter: grep -v marker | grep -v atomic | grep -v usage
739: test:
740: suffix: proj_tri_mdx_5P
741: requires: triangle !complex
742: args: -dim 2 -faces 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
743: filter: grep -v marker | grep -v atomic | grep -v usage
745: test:
746: suffix: proj_tri_3d
747: requires: ctetgen !complex
748: args: -dim 3 -faces 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
749: filter: grep -v marker | grep -v atomic | grep -v usage
751: test:
752: suffix: proj_tri_3d_2_faces
753: requires: ctetgen !complex
754: args: -dim 3 -faces 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
755: filter: grep -v marker | grep -v atomic | grep -v usage
757: test:
758: suffix: proj_tri_3d_5P
759: requires: ctetgen !complex
760: args: -dim 3 -faces 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
761: filter: grep -v marker | grep -v atomic | grep -v usage
763: test:
764: suffix: proj_tri_3d_mdx
765: requires: ctetgen !complex
766: args: -dim 3 -faces 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
767: filter: grep -v marker | grep -v atomic | grep -v usage
769: test:
770: suffix: proj_tri_3d_mdx_5P
771: requires: ctetgen !complex
772: args: -dim 3 -faces 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
773: filter: grep -v marker | grep -v atomic | grep -v usage
775: test:
776: suffix: proj_tri_3d_mdx_2_faces
777: requires: ctetgen !complex
778: args: -dim 3 -faces 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
779: filter: grep -v marker | grep -v atomic | grep -v usage
781: test:
782: suffix: proj_tri_3d_mdx_5P_2_faces
783: requires: ctetgen !complex
784: args: -dim 3 -faces 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
785: filter: grep -v marker | grep -v atomic | grep -v usage
787: TEST*/