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