Actual source code: dmproject.c
1: #include <petsc/private/dmimpl.h>
2: #include <petscdm.h>
3: #include <petscdmda.h>
4: #include <petscdmplex.h>
5: #include <petscdmswarm.h>
6: #include <petscksp.h>
7: #include <petscblaslapack.h>
9: #include <petsc/private/dmswarmimpl.h>
10: #include "../src/dm/impls/swarm/data_bucket.h" // For DataBucket internals
11: #include "petscmath.h"
13: typedef struct _projectConstraintsCtx {
14: DM dm;
15: Vec mask;
16: } projectConstraintsCtx;
18: static PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y)
19: {
20: DM dm;
21: Vec local, mask;
22: projectConstraintsCtx *ctx;
24: PetscFunctionBegin;
25: PetscCall(MatShellGetContext(CtC, &ctx));
26: dm = ctx->dm;
27: mask = ctx->mask;
28: PetscCall(DMGetLocalVector(dm, &local));
29: PetscCall(DMGlobalToLocalBegin(dm, x, INSERT_VALUES, local));
30: PetscCall(DMGlobalToLocalEnd(dm, x, INSERT_VALUES, local));
31: if (mask) PetscCall(VecPointwiseMult(local, mask, local));
32: PetscCall(VecSet(y, 0.));
33: PetscCall(DMLocalToGlobalBegin(dm, local, ADD_VALUES, y));
34: PetscCall(DMLocalToGlobalEnd(dm, local, ADD_VALUES, y));
35: PetscCall(DMRestoreLocalVector(dm, &local));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode DMGlobalToLocalSolve_project1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
40: {
41: PetscFunctionBegin;
42: for (PetscInt f = 0; f < Nf; f++) u[f] = 1.;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by `DMGlobalToLocalBegin()`/`DMGlobalToLocalEnd()` with mode
48: `INSERT_VALUES`.
50: Collective
52: Input Parameters:
53: + dm - The `DM` object
54: . x - The local vector
55: - y - The global vector: the input value of this variable is used as an initial guess
57: Output Parameter:
58: . y - The least-squares solution
60: Level: advanced
62: Note:
63: It is assumed that the sum of all the local vector sizes is greater than or equal to the global vector size, so the solution is
64: a least-squares solution. It is also assumed that `DMLocalToGlobalBegin()`/`DMLocalToGlobalEnd()` with mode `ADD_VALUES` is the adjoint of the
65: global-to-local map, so that the least-squares solution may be found by the normal equations.
67: If the `DM` is of type `DMPLEX`, then `y` is the solution of $ L^T * D * L * y = L^T * D * x $, where $D$ is a diagonal mask that is 1 for every point in
68: the union of the closures of the local cells and 0 otherwise. This difference is only relevant if there are anchor points that are not in the
69: closure of any local cell (see `DMPlexGetAnchors()`/`DMPlexSetAnchors()`).
71: What is L?
73: If this solves for a global vector from a local vector why is not called `DMLocalToGlobalSolve()`?
75: .seealso: [](ch_dmbase), `DM`, `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToGlobalEnd()`, `DMPlexGetAnchors()`, `DMPlexSetAnchors()`
76: @*/
77: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
78: {
79: Mat CtC;
80: PetscInt n, N, cStart, cEnd, c;
81: PetscBool isPlex;
82: KSP ksp;
83: PC pc;
84: Vec global, mask = NULL;
85: projectConstraintsCtx ctx;
87: PetscFunctionBegin;
88: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
89: if (isPlex) {
90: /* mark points in the closure */
91: PetscCall(DMCreateLocalVector(dm, &mask));
92: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
93: if (cEnd > cStart) {
94: PetscScalar *ones;
95: PetscInt numValues;
97: PetscCall(DMPlexVecGetClosure(dm, NULL, mask, cStart, &numValues, NULL));
98: PetscCall(PetscMalloc1(numValues, &ones));
99: for (PetscInt i = 0; i < numValues; i++) ones[i] = 1.;
100: for (c = cStart; c < cEnd; c++) PetscCall(DMPlexVecSetClosure(dm, NULL, mask, c, ones, INSERT_VALUES));
101: PetscCall(PetscFree(ones));
102: }
103: } else {
104: PetscBool hasMask;
106: PetscCall(DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask));
107: if (!hasMask) {
108: PetscErrorCode (**func)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
109: void **ctx;
110: PetscInt Nf;
112: PetscCall(DMGetNumFields(dm, &Nf));
113: PetscCall(PetscMalloc2(Nf, &func, Nf, &ctx));
114: for (PetscInt f = 0; f < Nf; ++f) {
115: func[f] = DMGlobalToLocalSolve_project1;
116: ctx[f] = NULL;
117: }
118: PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
119: PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctx, INSERT_ALL_VALUES, mask));
120: PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
121: PetscCall(PetscFree2(func, ctx));
122: }
123: PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
124: }
125: ctx.dm = dm;
126: ctx.mask = mask;
127: PetscCall(VecGetSize(y, &N));
128: PetscCall(VecGetLocalSize(y, &n));
129: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &CtC));
130: PetscCall(MatSetSizes(CtC, n, n, N, N));
131: PetscCall(MatSetType(CtC, MATSHELL));
132: PetscCall(MatSetUp(CtC));
133: PetscCall(MatShellSetContext(CtC, &ctx));
134: PetscCall(MatShellSetOperation(CtC, MATOP_MULT, (PetscErrorCodeFn *)MatMult_GlobalToLocalNormal));
135: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
136: PetscCall(KSPSetOperators(ksp, CtC, CtC));
137: PetscCall(KSPSetType(ksp, KSPCG));
138: PetscCall(KSPGetPC(ksp, &pc));
139: PetscCall(PCSetType(pc, PCNONE));
140: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
141: PetscCall(KSPSetUp(ksp));
142: PetscCall(DMGetGlobalVector(dm, &global));
143: PetscCall(VecSet(global, 0.));
144: if (mask) PetscCall(VecPointwiseMult(x, mask, x));
145: PetscCall(DMLocalToGlobalBegin(dm, x, ADD_VALUES, global));
146: PetscCall(DMLocalToGlobalEnd(dm, x, ADD_VALUES, global));
147: PetscCall(KSPSolve(ksp, global, y));
148: PetscCall(DMRestoreGlobalVector(dm, &global));
149: /* clean up */
150: PetscCall(KSPDestroy(&ksp));
151: PetscCall(MatDestroy(&CtC));
152: if (isPlex) {
153: PetscCall(VecDestroy(&mask));
154: } else {
155: PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
156: }
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: /*@
161: DMProjectField - This projects a given function of the input fields into the function space provided by a `DM`, putting the coefficients in a global vector.
163: Collective
165: Input Parameters:
166: + dm - The `DM`
167: . time - The time
168: . U - The input field vector
169: . funcs - The functions to evaluate, one per field, see `PetscPointFn`
170: - mode - The insertion mode for values
172: Output Parameter:
173: . X - The output vector
175: Level: advanced
177: Note:
178: There are three different `DM`s that potentially interact in this function. The output `dm`, specifies the layout of the values calculates by the function.
179: The input `DM`, attached to `U`, may be different. For example, you can input the solution over the full domain, but output over a piece of the boundary, or
180: a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary `DM`, attached to the
181: auxiliary field vector, which is attached to `dm`, can also be different. It can have a different topology, number of fields, and discretizations.
183: .seealso: [](ch_dmbase), `DM`, `PetscPointFn`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`, `DMProjectFunction()`, `DMComputeL2Diff()`
184: @*/
185: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U, PetscPointFn **funcs, InsertMode mode, Vec X)
186: {
187: Vec localX, localU;
188: DM dmIn;
190: PetscFunctionBegin;
192: PetscCall(DMGetLocalVector(dm, &localX));
193: /* We currently check whether locU == locX to see if we need to apply BC */
194: if (U != X) {
195: PetscCall(VecGetDM(U, &dmIn));
196: PetscCall(DMGetLocalVector(dmIn, &localU));
197: } else {
198: dmIn = dm;
199: localU = localX;
200: }
201: PetscCall(DMGlobalToLocalBegin(dmIn, U, INSERT_VALUES, localU));
202: PetscCall(DMGlobalToLocalEnd(dmIn, U, INSERT_VALUES, localU));
203: PetscCall(DMProjectFieldLocal(dm, time, localU, funcs, mode, localX));
204: PetscCall(DMLocalToGlobalBegin(dm, localX, mode, X));
205: PetscCall(DMLocalToGlobalEnd(dm, localX, mode, X));
206: if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
207: Mat cMat;
209: PetscCall(DMGetDefaultConstraints(dm, NULL, &cMat, NULL));
210: if (cMat) PetscCall(DMGlobalToLocalSolve(dm, localX, X));
211: }
212: PetscCall(DMRestoreLocalVector(dm, &localX));
213: if (U != X) PetscCall(DMRestoreLocalVector(dmIn, &localU));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /********************* Adaptive Interpolation **************************/
219: /* See the discussion of Adaptive Interpolation in manual/high_level_mg.rst */
220: /*@
221: DMAdaptInterpolator - Adapts a grid interpolator so that it accurately reproduces a set of sample fine-grid vectors
223: Collective
225: Input Parameters:
226: + dmc - the coarse `DM`
227: . dmf - the fine `DM`
228: . In - the input (unadapted) interpolation matrix from `dmc` to `dmf`
229: . smoother - a `KSP` whose operator provides the fine-grid matrix used to weight modes by their Rayleigh quotient
230: . MF - a dense matrix whose columns are fine-grid sample vectors
231: . MC - a dense matrix whose columns are the corresponding coarse-grid sample vectors (may be `NULL`, in which case $I_n^T M_F$ is used)
232: - user - unused application context
234: Output Parameter:
235: . InAdapt - the adapted interpolation matrix (created inside the routine)
237: Options Database Key:
238: . -dm_interpolator_adapt_debug flag - print diagnostic information about the least-squares systems solved for each row
240: Level: developer
242: Note:
243: For each row of `In` a small weighted least-squares problem is solved (using LAPACK GELSS) so that the adapted
244: interpolation reproduces the fine-grid samples as accurately as possible; see the discussion of adaptive interpolation
245: in `manual/high_level_mg.rst`.
247: .seealso: [](ch_ksp), `DM`, `Mat`, `KSP`, `DMCheckInterpolator()`, `DMCreateInterpolation()`, `PCMG`
248: @*/
249: PetscErrorCode DMAdaptInterpolator(DM dmc, DM dmf, Mat In, KSP smoother, Mat MF, Mat MC, Mat *InAdapt, void *user)
250: {
251: Mat globalA, AF;
252: Vec tmp;
253: const PetscScalar *af, *ac;
254: PetscScalar *A, *b, *x, *workscalar;
255: PetscReal *w, *sing, *workreal, rcond = PETSC_SMALL;
256: PetscBLASInt M, N, one = 1, irank, lwrk;
257: PetscInt debug = 0, rStart, rEnd, r, maxcols = 0, k, Nc, ldac, ldaf;
258: PetscBool allocVc = PETSC_FALSE;
260: PetscFunctionBegin;
261: PetscCall(PetscLogEventBegin(DM_AdaptInterpolator, dmc, dmf, 0, 0));
262: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dm_interpolator_adapt_debug", &debug, NULL));
263: PetscCall(MatGetSize(MF, NULL, &Nc));
264: PetscCall(MatDuplicate(In, MAT_SHARE_NONZERO_PATTERN, InAdapt));
265: PetscCall(MatGetOwnershipRange(In, &rStart, &rEnd));
266: #if 0
267: PetscCall(MatGetMaxRowLen(In, &maxcols));
268: #else
269: for (r = rStart; r < rEnd; ++r) {
270: PetscInt ncols;
272: PetscCall(MatGetRow(In, r, &ncols, NULL, NULL));
273: maxcols = PetscMax(maxcols, ncols);
274: PetscCall(MatRestoreRow(In, r, &ncols, NULL, NULL));
275: }
276: #endif
277: if (Nc < maxcols) PetscCall(PetscPrintf(PETSC_COMM_SELF, "The number of input vectors %" PetscInt_FMT " < %" PetscInt_FMT " the maximum number of column entries\n", Nc, maxcols));
278: for (k = 0; k < Nc && debug; ++k) {
279: char name[PETSC_MAX_PATH_LEN];
280: const char *prefix;
281: Vec vc, vf;
283: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)smoother, &prefix));
285: if (MC) {
286: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sCoarse Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
287: PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
288: PetscCall(PetscObjectSetName((PetscObject)vc, name));
289: PetscCall(VecViewFromOptions(vc, NULL, "-dm_adapt_interp_view_coarse"));
290: PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
291: }
292: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sFine Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
293: PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
294: PetscCall(PetscObjectSetName((PetscObject)vf, name));
295: PetscCall(VecViewFromOptions(vf, NULL, "-dm_adapt_interp_view_fine"));
296: PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
297: }
298: PetscCall(PetscBLASIntCast(3 * PetscMin(Nc, maxcols) + PetscMax(2 * PetscMin(Nc, maxcols), PetscMax(Nc, maxcols)), &lwrk));
299: PetscCall(PetscMalloc7(Nc * maxcols, &A, PetscMax(Nc, maxcols), &b, Nc, &w, maxcols, &x, maxcols, &sing, lwrk, &workscalar, 5 * PetscMin(Nc, maxcols), &workreal));
300: /* w_k = \frac{\HC{v_k} B_l v_k}{\HC{v_k} A_l v_k} or the inverse Rayleigh quotient, which we calculate using \frac{\HC{v_k} v_k}{\HC{v_k} B^{-1}_l A_l v_k} */
301: PetscCall(KSPGetOperators(smoother, &globalA, NULL));
303: PetscCall(MatMatMult(globalA, MF, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AF));
304: for (k = 0; k < Nc; ++k) {
305: PetscScalar vnorm, vAnorm;
306: Vec vf;
308: w[k] = 1.0;
309: PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
310: PetscCall(MatDenseGetColumnVecRead(AF, k, &tmp));
311: PetscCall(VecDot(vf, vf, &vnorm));
312: #if 0
313: PetscCall(DMGetGlobalVector(dmf, &tmp2));
314: PetscCall(KSPSolve(smoother, tmp, tmp2));
315: PetscCall(VecDot(vf, tmp2, &vAnorm));
316: PetscCall(DMRestoreGlobalVector(dmf, &tmp2));
317: #else
318: PetscCall(VecDot(vf, tmp, &vAnorm));
319: #endif
320: w[k] = PetscRealPart(vnorm) / PetscRealPart(vAnorm);
321: PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
322: PetscCall(MatDenseRestoreColumnVecRead(AF, k, &tmp));
323: }
324: PetscCall(MatDestroy(&AF));
325: if (!MC) {
326: allocVc = PETSC_TRUE;
327: PetscCall(MatTransposeMatMult(In, MF, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &MC));
328: }
329: /* Solve a LS system for each fine row
330: MATT: Can we generalize to the case where Nc for the fine space
331: is different for Nc for the coarse? */
332: PetscCall(MatDenseGetArrayRead(MF, &af));
333: PetscCall(MatDenseGetLDA(MF, &ldaf));
334: PetscCall(MatDenseGetArrayRead(MC, &ac));
335: PetscCall(MatDenseGetLDA(MC, &ldac));
336: for (r = rStart; r < rEnd; ++r) {
337: PetscInt ncols;
338: const PetscInt *cols;
339: const PetscScalar *vals;
341: PetscCall(MatGetRow(In, r, &ncols, &cols, &vals));
342: for (k = 0; k < Nc; ++k) {
343: /* Need to fit lowest mode exactly */
344: const PetscReal wk = ((ncols == 1) && (k > 0)) ? 0.0 : PetscSqrtReal(w[k]);
346: /* b_k = \sqrt{w_k} f^{F,k}_r */
347: b[k] = wk * af[r - rStart + k * ldaf];
348: /* A_{kc} = \sqrt{w_k} f^{C,k}_c */
349: /* TODO Must pull out VecScatter from In, scatter in vc[k] values up front, and access them indirectly just as in MatMult() */
350: for (PetscInt c = 0; c < ncols; ++c) {
351: /* This is element (k, c) of A */
352: A[c * Nc + k] = wk * ac[cols[c] - rStart + k * ldac];
353: }
354: }
355: PetscCall(PetscBLASIntCast(Nc, &M));
356: PetscCall(PetscBLASIntCast(ncols, &N));
357: if (debug) {
358: #if PetscDefined(USE_COMPLEX)
359: PetscScalar *tmp;
361: PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
362: for (PetscInt j = 0; j < Nc; ++j) tmp[j] = w[j];
363: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, tmp));
364: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
365: for (PetscInt j = 0; j < Nc; ++j) tmp[j] = b[j];
366: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, tmp));
367: PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
368: #else
369: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, w));
370: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
371: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, b));
372: #endif
373: }
374: #if PetscDefined(USE_COMPLEX)
375: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, workreal, &info));
376: #else
377: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, &info));
378: #endif
379: if (debug) {
380: PetscCall(PetscPrintf(PETSC_COMM_SELF, "rank %" PetscBLASInt_FMT " rcond %g\n", irank, (double)rcond));
381: #if PetscDefined(USE_COMPLEX)
382: {
383: PetscScalar *tmp;
385: PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
386: for (PetscInt j = 0; j < PetscMin(Nc, ncols); ++j) tmp[j] = sing[j];
387: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, tmp));
388: PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
389: }
390: #else
391: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, sing));
392: #endif
393: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS old P", ncols, 1, vals));
394: PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS sol", ncols, 1, b));
395: }
396: PetscCall(MatSetValues(*InAdapt, 1, &r, ncols, cols, b, INSERT_VALUES));
397: PetscCall(MatRestoreRow(In, r, &ncols, &cols, &vals));
398: }
399: PetscCall(MatDenseRestoreArrayRead(MF, &af));
400: PetscCall(MatDenseRestoreArrayRead(MC, &ac));
401: PetscCall(PetscFree7(A, b, w, x, sing, workscalar, workreal));
402: if (allocVc) PetscCall(MatDestroy(&MC));
403: PetscCall(MatAssemblyBegin(*InAdapt, MAT_FINAL_ASSEMBLY));
404: PetscCall(MatAssemblyEnd(*InAdapt, MAT_FINAL_ASSEMBLY));
405: PetscCall(PetscLogEventEnd(DM_AdaptInterpolator, dmc, dmf, 0, 0));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*@
410: DMCheckInterpolator - Check that an interpolation matrix accurately reproduces a set of sample fine-grid vectors
412: Collective
414: Input Parameters:
415: + dmf - the fine `DM`
416: . In - the interpolation matrix from a coarse `DM` to `dmf`
417: . MC - a dense matrix whose columns are coarse-grid sample vectors
418: . MF - a dense matrix whose columns are the corresponding fine-grid sample vectors
419: - tol - tolerance on the maximum 2-norm of $v_f - I v_c$ across all sample vectors
421: Options Database Key:
422: . -dm_interpolator_adapt_error view - view the coarse, fine, and error vectors for each sample
424: Level: developer
426: Note:
427: For each column `k`, the residual $v_f^k - I v_c^k$ is computed and its infinity and 2 norms are printed. An error
428: is raised if the maximum 2-norm exceeds `tol`. Typically used with `DMAdaptInterpolator()` to validate the adapted operator.
430: .seealso: [](ch_ksp), `DM`, `Mat`, `DMAdaptInterpolator()`, `DMCreateInterpolation()`, `PCMG`
431: @*/
432: PetscErrorCode DMCheckInterpolator(DM dmf, Mat In, Mat MC, Mat MF, PetscReal tol)
433: {
434: Vec tmp;
435: PetscReal norminf, norm2, maxnorminf = 0.0, maxnorm2 = 0.0;
436: PetscInt k, Nc;
438: PetscFunctionBegin;
439: PetscCall(DMGetGlobalVector(dmf, &tmp));
440: PetscCall(MatViewFromOptions(In, NULL, "-dm_interpolator_adapt_error"));
441: PetscCall(MatGetSize(MF, NULL, &Nc));
442: for (k = 0; k < Nc; ++k) {
443: Vec vc, vf;
445: PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
446: PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
447: PetscCall(MatMult(In, vc, tmp));
448: PetscCall(VecAXPY(tmp, -1.0, vf));
449: PetscCall(VecViewFromOptions(vc, NULL, "-dm_interpolator_adapt_error"));
450: PetscCall(VecViewFromOptions(vf, NULL, "-dm_interpolator_adapt_error"));
451: PetscCall(VecViewFromOptions(tmp, NULL, "-dm_interpolator_adapt_error"));
452: PetscCall(VecNorm(tmp, NORM_INFINITY, &norminf));
453: PetscCall(VecNorm(tmp, NORM_2, &norm2));
454: maxnorminf = PetscMax(maxnorminf, norminf);
455: maxnorm2 = PetscMax(maxnorm2, norm2);
456: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmf), "Coarse vec %" PetscInt_FMT " ||vf - P vc||_\\infty %g, ||vf - P vc||_2 %g\n", k, (double)norminf, (double)norm2));
457: PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
458: PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
459: }
460: PetscCall(DMRestoreGlobalVector(dmf, &tmp));
461: PetscCheck(maxnorm2 <= tol, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_WRONG, "max_k ||vf_k - P vc_k||_2 %g > tol %g", (double)maxnorm2, (double)tol);
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: // Project particles to field
466: // M_f u_f = M_p u_p
467: // u_f = M^{-1}_f M_p u_p
468: static PetscErrorCode DMSwarmProjectField_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
469: {
470: KSP ksp;
471: Mat M_f, M_p; // TODO Should cache these
472: Vec rhs;
473: const char *prefix;
475: PetscFunctionBegin;
476: PetscCall(DMCreateMassMatrix(dm, dm, &M_f));
477: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
478: PetscCall(DMGetGlobalVector(dm, &rhs));
479: PetscCall(MatMultTranspose(M_p, u_p, rhs));
481: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
482: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
483: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
484: PetscCall(KSPAppendOptionsPrefix(ksp, "ptof_"));
485: PetscCall(KSPSetFromOptions(ksp));
487: PetscCall(KSPSetOperators(ksp, M_f, M_f));
488: PetscCall(KSPSolve(ksp, rhs, u_f));
490: PetscCall(DMRestoreGlobalVector(dm, &rhs));
491: PetscCall(KSPDestroy(&ksp));
492: PetscCall(MatDestroy(&M_f));
493: PetscCall(MatDestroy(&M_p));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: // Project field to particles
498: // M_p u_p = M_f u_f
499: // u_p = M^+_p M_f u_f
500: static PetscErrorCode DMSwarmProjectParticles_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
501: {
502: KSP ksp;
503: PC pc;
504: Mat M_f, M_p, PM_p;
505: Vec rhs;
506: PetscBool isBjacobi;
507: const char *prefix;
509: PetscFunctionBegin;
510: PetscCall(DMCreateMassMatrix(dm, dm, &M_f));
511: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
512: PetscCall(DMGetGlobalVector(dm, &rhs));
513: PetscCall(MatMult(M_f, u_f, rhs));
515: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
516: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
517: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
518: PetscCall(KSPAppendOptionsPrefix(ksp, "ftop_"));
519: PetscCall(KSPSetFromOptions(ksp));
521: PetscCall(KSPGetPC(ksp, &pc));
522: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &isBjacobi));
523: if (isBjacobi) {
524: PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p));
525: } else {
526: PM_p = M_p;
527: PetscCall(PetscObjectReference((PetscObject)PM_p));
528: }
529: PetscCall(KSPSetOperators(ksp, M_p, PM_p));
530: PetscCall(KSPSolveTranspose(ksp, rhs, u_p));
532: PetscCall(DMRestoreGlobalVector(dm, &rhs));
533: PetscCall(KSPDestroy(&ksp));
534: PetscCall(MatDestroy(&M_f));
535: PetscCall(MatDestroy(&M_p));
536: PetscCall(MatDestroy(&PM_p));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode DMSwarmProjectFields_Plex_Internal(DM sw, DM dm, PetscInt Nf, const char *fieldnames[], Vec vec, ScatterMode mode)
541: {
542: PetscDS ds;
543: Vec u;
544: PetscInt f = 0, bs, *Nc;
546: PetscFunctionBegin;
547: PetscCall(DMGetDS(dm, &ds));
548: PetscCall(PetscDSGetComponents(ds, &Nc));
549: PetscCall(PetscCitationsRegister(SwarmProjCitation, &SwarmProjcite));
550: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Currently supported only for a single field");
551: PetscCall(DMSwarmVectorDefineFields(sw, Nf, fieldnames));
552: PetscCall(DMSwarmCreateGlobalVectorFromFields(sw, Nf, fieldnames, &u));
553: PetscCall(VecGetBlockSize(u, &bs));
554: PetscCheck(Nc[f] == bs, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %" PetscInt_FMT " components %" PetscInt_FMT " != %" PetscInt_FMT " blocksize for swarm field %s", f, Nc[f], bs, fieldnames[f]);
555: if (mode == SCATTER_FORWARD) {
556: PetscCall(DMSwarmProjectField_Conservative_PLEX(sw, dm, u, vec));
557: } else {
558: PetscCall(DMSwarmProjectParticles_Conservative_PLEX(sw, dm, u, vec));
559: }
560: PetscCall(DMSwarmDestroyGlobalVectorFromFields(sw, Nf, fieldnames, &u));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
565: {
566: DMSwarmCellDM celldm;
567: Vec v_field_l, denom_l, coor_l, denom;
568: PetscScalar *_field_l, *_denom_l;
569: PetscInt k, p, e, npoints, nel, npe, Nfc;
570: PetscInt *mpfield_cell;
571: PetscReal *mpfield_coor;
572: const PetscInt *element_list;
573: const PetscInt *element;
574: PetscScalar xi_p[2], Ni[4];
575: const PetscScalar *_coor;
576: const char **coordFields, *cellid;
578: PetscFunctionBegin;
579: PetscCall(VecZeroEntries(v_field));
581: PetscCall(DMGetLocalVector(dm, &v_field_l));
582: PetscCall(DMGetGlobalVector(dm, &denom));
583: PetscCall(DMGetLocalVector(dm, &denom_l));
584: PetscCall(VecZeroEntries(v_field_l));
585: PetscCall(VecZeroEntries(denom));
586: PetscCall(VecZeroEntries(denom_l));
588: PetscCall(VecGetArray(v_field_l, &_field_l));
589: PetscCall(VecGetArray(denom_l, &_denom_l));
591: PetscCall(DMGetCoordinatesLocal(dm, &coor_l));
592: PetscCall(VecGetArrayRead(coor_l, &_coor));
594: PetscCall(DMSwarmGetCellDMActive(swarm, &celldm));
595: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
596: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
597: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
599: PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
600: PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
601: PetscCall(DMSwarmGetField(swarm, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
602: PetscCall(DMSwarmGetField(swarm, cellid, NULL, NULL, (void **)&mpfield_cell));
604: for (p = 0; p < npoints; p++) {
605: PetscReal *coor_p;
606: const PetscScalar *x0;
607: const PetscScalar *x2;
608: PetscScalar dx[2];
610: e = mpfield_cell[p];
611: coor_p = &mpfield_coor[2 * p];
612: element = &element_list[npe * e];
614: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
615: x0 = &_coor[2 * element[0]];
616: x2 = &_coor[2 * element[2]];
618: dx[0] = x2[0] - x0[0];
619: dx[1] = x2[1] - x0[1];
621: xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
622: xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
624: /* evaluate basis functions */
625: Ni[0] = 0.25 * (1.0 - xi_p[0]) * (1.0 - xi_p[1]);
626: Ni[1] = 0.25 * (1.0 + xi_p[0]) * (1.0 - xi_p[1]);
627: Ni[2] = 0.25 * (1.0 + xi_p[0]) * (1.0 + xi_p[1]);
628: Ni[3] = 0.25 * (1.0 - xi_p[0]) * (1.0 + xi_p[1]);
630: for (k = 0; k < npe; k++) {
631: _field_l[element[k]] += Ni[k] * swarm_field[p];
632: _denom_l[element[k]] += Ni[k];
633: }
634: }
636: PetscCall(DMSwarmRestoreField(swarm, cellid, NULL, NULL, (void **)&mpfield_cell));
637: PetscCall(DMSwarmRestoreField(swarm, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
638: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
639: PetscCall(VecRestoreArrayRead(coor_l, &_coor));
640: PetscCall(VecRestoreArray(v_field_l, &_field_l));
641: PetscCall(VecRestoreArray(denom_l, &_denom_l));
643: PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field));
644: PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field));
645: PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom));
646: PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom));
648: PetscCall(VecPointwiseDivide(v_field, v_field, denom));
650: PetscCall(DMRestoreLocalVector(dm, &v_field_l));
651: PetscCall(DMRestoreLocalVector(dm, &denom_l));
652: PetscCall(DMRestoreGlobalVector(dm, &denom));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode DMSwarmProjectFields_DA_Internal(DM swarm, DM celldm, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[], ScatterMode mode)
657: {
658: PetscInt dim;
659: DMDAElementType etype;
661: PetscFunctionBegin;
662: PetscCall(DMDAGetElementType(celldm, &etype));
663: PetscCheck(etype != DMDA_ELEMENT_P1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Only Q1 DMDA supported");
664: PetscCheck(mode == SCATTER_FORWARD, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Mapping the continuum to particles is not currently supported for DMDA");
666: PetscCall(DMGetDimension(swarm, &dim));
667: switch (dim) {
668: case 2:
669: for (PetscInt f = 0; f < nfields; f++) {
670: PetscReal *swarm_field;
672: PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field));
673: PetscCall(DMSwarmProjectField_ApproxQ1_DA_2D(swarm, swarm_field, celldm, vecs[f]));
674: }
675: break;
676: case 3:
677: SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
678: default:
679: break;
680: }
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: /*@
685: DMSwarmProjectFields - Project a set of swarm fields onto another `DM`
687: Collective
689: Input Parameters:
690: + sw - the `DMSWARM`
691: . dm - the `DM`, or `NULL` to use the cell `DM`
692: . nfields - the number of swarm fields to project
693: . fieldnames - the textual names of the swarm fields to project
694: . fields - an array of `Vec`'s of length nfields
695: - mode - if `SCATTER_FORWARD` then map particles to the continuum, and if `SCATTER_REVERSE` map the continuum to particles
697: Level: beginner
699: Notes:
700: Currently, there are two available projection methods. The first is conservative projection, used for a `DMPLEX` cell `DM`.
701: The second is the averaging which is used for a `DMDA` cell `DM`
703: $$
704: \phi_i = \sum_{p=0}^{np} N_i(x_p) \phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
705: $$
707: where $\phi_p $ is the swarm field at point $p$, $N_i()$ is the cell `DM` basis function at vertex $i$, $dJ$ is the determinant of the cell Jacobian and
708: $\phi_i$ is the projected vertex value of the field $\phi$.
710: The user is responsible for destroying both the array and the individual `Vec` objects.
712: For the `DMPLEX` case, there is only a single vector, so the field layout in the `DMPLEX` must match the requested fields from the `DMSwarm`.
714: For averaging projection, nly swarm fields registered with data type of `PETSC_REAL` can be projected onto the cell `DM`, and only swarm fields of block size = 1 can currently be projected.
716: .seealso: [](ch_dmbase), `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
717: @*/
718: PetscErrorCode DMSwarmProjectFields(DM sw, DM dm, PetscInt nfields, const char *fieldnames[], Vec fields[], ScatterMode mode)
719: {
720: DM_Swarm *swarm = (DM_Swarm *)sw->data;
721: DMSwarmDataField *gfield;
722: PetscBool isDA, isPlex;
723: MPI_Comm comm;
725: PetscFunctionBegin;
726: DMSWARMPICVALID(sw);
727: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
728: if (!dm) PetscCall(DMSwarmGetCellDM(sw, &dm));
729: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isDA));
730: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
731: PetscCall(PetscMalloc1(nfields, &gfield));
732: for (PetscInt f = 0; f < nfields; ++f) PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldnames[f], &gfield[f]));
734: if (isDA) {
735: for (PetscInt f = 0; f < nfields; f++) {
736: PetscCheck(gfield[f]->petsc_type == PETSC_REAL, comm, PETSC_ERR_SUP, "Projection only valid for fields using a data type = PETSC_REAL");
737: PetscCheck(gfield[f]->bs == 1, comm, PETSC_ERR_SUP, "Projection only valid for fields with block size = 1");
738: }
739: PetscCall(DMSwarmProjectFields_DA_Internal(sw, dm, nfields, gfield, fields, mode));
740: } else if (isPlex) {
741: PetscInt Nf;
743: PetscCall(DMGetNumFields(dm, &Nf));
744: PetscCheck(Nf == nfields, comm, PETSC_ERR_ARG_WRONG, "Number of DM fields %" PetscInt_FMT " != %" PetscInt_FMT " number of requested Swarm fields", Nf, nfields);
745: PetscCall(DMSwarmProjectFields_Plex_Internal(sw, dm, nfields, fieldnames, fields[0], mode));
746: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
748: PetscCall(PetscFree(gfield));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: // Project weak divergence of particles to field
753: // \int_X psi_i div u_f = \int_X psi_i div u_p
754: // \int_X grad psi_i . \sum_j u_f \psi_j = \int_X grad psi_i . \sum_p u_p \delta(x - x_p)
755: // D_f u_f = D_p u_p
756: // u_f = D^+_f D_p u_p
757: static PetscErrorCode DMSwarmProjectGradientField_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
758: {
759: DM gdm;
760: KSP ksp;
761: Mat D_f, D_p; // TODO Should cache these
762: Vec rhs;
763: const char *prefix;
765: PetscFunctionBegin;
766: PetscCall(VecGetDM(u_f, &gdm));
767: PetscCall(DMCreateGradientMatrix(dm, gdm, &D_f));
768: PetscCall(DMCreateGradientMatrix(sw, dm, &D_p));
769: PetscCall(DMGetGlobalVector(dm, &rhs));
770: PetscCall(PetscObjectSetName((PetscObject)rhs, "D u"));
771: PetscCall(MatMultTranspose(D_p, u_p, rhs));
772: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
774: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
775: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
776: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
777: PetscCall(KSPAppendOptionsPrefix(ksp, "gptof_"));
778: PetscCall(KSPSetFromOptions(ksp));
780: PetscCall(KSPSetOperators(ksp, D_f, D_f));
781: PetscCall(KSPSolveTranspose(ksp, rhs, u_f));
783: PetscCall(MatMultTranspose(D_f, u_f, rhs));
784: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
786: PetscCall(DMRestoreGlobalVector(dm, &rhs));
787: PetscCall(KSPDestroy(&ksp));
788: PetscCall(MatDestroy(&D_f));
789: PetscCall(MatDestroy(&D_p));
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: // Project weak divergence of field to particles
794: // D_p u_p = D_f u_f
795: // u_p = D^+_p D_f u_f
796: static PetscErrorCode DMSwarmProjectGradientParticles_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
797: {
798: KSP ksp;
799: PC pc;
800: Mat D_f, D_p, PD_p;
801: Vec rhs;
802: PetscBool isBjacobi;
803: const char *prefix;
805: PetscFunctionBegin;
806: PetscCall(DMCreateGradientMatrix(dm, dm, &D_f));
807: PetscCall(DMCreateGradientMatrix(sw, dm, &D_p));
808: PetscCall(DMGetGlobalVector(dm, &rhs));
809: PetscCall(MatMult(D_f, u_f, rhs));
811: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
812: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
813: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
814: PetscCall(KSPAppendOptionsPrefix(ksp, "gftop_"));
815: PetscCall(KSPSetFromOptions(ksp));
817: PetscCall(KSPGetPC(ksp, &pc));
818: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &isBjacobi));
819: if (isBjacobi) {
820: PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PD_p));
821: } else {
822: PD_p = D_p;
823: PetscCall(PetscObjectReference((PetscObject)PD_p));
824: }
825: PetscCall(KSPSetOperators(ksp, D_p, PD_p));
826: PetscCall(KSPSolveTranspose(ksp, rhs, u_p));
828: PetscCall(DMRestoreGlobalVector(dm, &rhs));
829: PetscCall(KSPDestroy(&ksp));
830: PetscCall(MatDestroy(&D_f));
831: PetscCall(MatDestroy(&D_p));
832: PetscCall(MatDestroy(&PD_p));
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: static PetscErrorCode DMSwarmProjectGradientFields_Plex_Internal(DM sw, DM dm, PetscInt Nf, const char *fieldnames[], Vec vec, ScatterMode mode)
837: {
838: PetscDS ds;
839: Vec u;
840: PetscInt f = 0, cdim, bs, *Nc;
842: PetscFunctionBegin;
843: PetscCall(DMGetCoordinateDim(dm, &cdim));
844: PetscCall(DMGetDS(dm, &ds));
845: PetscCall(PetscDSGetComponents(ds, &Nc));
846: PetscCall(PetscCitationsRegister(SwarmProjCitation, &SwarmProjcite));
847: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Currently supported only for a single field");
848: PetscCall(DMSwarmVectorDefineFields(sw, Nf, fieldnames));
849: PetscCall(DMSwarmCreateGlobalVectorFromFields(sw, Nf, fieldnames, &u));
850: PetscCall(VecGetBlockSize(u, &bs));
851: PetscCheck(Nc[f] * cdim == bs, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %" PetscInt_FMT " components %" PetscInt_FMT " * %" PetscInt_FMT " coordinate dim != %" PetscInt_FMT " blocksize for swarm field %s", f, Nc[f], cdim, bs, fieldnames[f]);
852: if (mode == SCATTER_FORWARD) {
853: PetscCall(DMSwarmProjectGradientField_Conservative_PLEX(sw, dm, u, vec));
854: } else {
855: PetscCall(DMSwarmProjectGradientParticles_Conservative_PLEX(sw, dm, u, vec));
856: }
857: PetscCall(DMSwarmDestroyGlobalVectorFromFields(sw, Nf, fieldnames, &u));
858: PetscFunctionReturn(PETSC_SUCCESS);
859: }
861: /*@
862: DMSwarmProjectGradientFields - Project the gradient of continuum fields on a mesh onto particle fields in a `DMSWARM`, or the reverse
864: Collective
866: Input Parameters:
867: + sw - the `DMSWARM`
868: . dm - the continuum `DM` (a `DMPLEX`); if `NULL` the swarm's cell `DM` is used
869: . nfields - the number of fields to project
870: . fieldnames - the names of the swarm fields to receive (or supply) the gradient
871: . fields - the corresponding mesh `Vec` objects
872: - mode - `SCATTER_FORWARD` to project mesh field gradients to particles, `SCATTER_REVERSE` to project particle values back to the mesh
874: Level: intermediate
876: Note:
877: Only `DMPLEX` cell DMs and single-field projection are currently supported. The swarm field block size must equal
878: the mesh field component count times the coordinate dimension.
880: .seealso: `DMSWARM`, `DMPLEX`, `DMSwarmProjectFields()`, `DMSwarmVectorDefineFields()`, `DMSwarmCreateGlobalVectorFromField()`
881: @*/
882: PetscErrorCode DMSwarmProjectGradientFields(DM sw, DM dm, PetscInt nfields, const char *fieldnames[], Vec fields[], ScatterMode mode)
883: {
884: PetscBool isPlex;
885: MPI_Comm comm;
887: PetscFunctionBegin;
888: DMSWARMPICVALID(sw);
889: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
890: if (!dm) PetscCall(DMSwarmGetCellDM(sw, &dm));
891: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
892: if (isPlex) {
893: PetscInt Nf;
895: PetscCall(DMGetNumFields(dm, &Nf));
896: PetscCheck(Nf == nfields, comm, PETSC_ERR_ARG_WRONG, "Number of DM fields %" PetscInt_FMT " != %" PetscInt_FMT " number of requested Swarm fields", Nf, nfields);
897: PetscCall(DMSwarmProjectGradientFields_Plex_Internal(sw, dm, nfields, fieldnames, fields[0], mode));
898: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only supported for cell DMs of type DMPLEX");
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: /*
903: InitializeParticles_Regular - Initialize a regular grid of particles in each cell
905: Input Parameters:
906: + sw - The `DMSWARM`
907: - n - The number of particles per dimension per species
909: Notes:
910: This functions sets the species, cellid, and cell DM coordinates.
912: It places n^d particles per species in each cell of the cell DM.
913: */
914: static PetscErrorCode InitializeParticles_Regular(DM sw, PetscInt n)
915: {
916: DM_Swarm *swarm = (DM_Swarm *)sw->data;
917: DM dm;
918: DMSwarmCellDM celldm;
919: PetscInt dim, Ns, Npc, Np, cStart, cEnd, debug;
920: PetscBool flg;
921: MPI_Comm comm;
923: PetscFunctionBegin;
924: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
926: PetscOptionsBegin(comm, "", "DMSwarm Options", "DMSWARM");
927: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
928: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
929: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
930: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_coords", "Debug output level for particle coordinate computations", "InitializeParticles", 0, &swarm->printCoords, NULL, 0));
931: PetscCall(PetscOptionsBoundedInt("-dm_swarm_print_weights", "Debug output level for particle weight computations", "InitializeWeights", 0, &swarm->printWeights, NULL, 0));
932: PetscOptionsEnd();
933: debug = swarm->printCoords;
935: // n^d particle per cell on the grid
936: PetscCall(DMSwarmGetCellDM(sw, &dm));
937: PetscCall(DMGetDimension(dm, &dim));
938: PetscCheck(!(dim % 2), comm, PETSC_ERR_SUP, "We only support even dimension, not %" PetscInt_FMT, dim);
939: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
940: Npc = Ns * PetscPowInt(n, dim);
941: Np = (cEnd - cStart) * Npc;
942: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
943: if (debug) {
944: PetscInt gNp;
945: PetscCallMPI(MPIU_Allreduce(&Np, &gNp, 1, MPIU_INT, MPIU_SUM, comm));
946: PetscCall(PetscPrintf(comm, "Global Np = %" PetscInt_FMT "\n", gNp));
947: }
948: PetscCall(PetscPrintf(comm, "Regular layout using %" PetscInt_FMT " particles per cell\n", Npc));
950: // Set species and cellid
951: {
952: const char *cellidName;
953: PetscInt *species, *cellid;
955: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
956: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellidName));
957: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
958: PetscCall(DMSwarmGetField(sw, cellidName, NULL, NULL, (void **)&cellid));
959: for (PetscInt c = 0, p = 0; c < cEnd - cStart; ++c) {
960: for (PetscInt s = 0; s < Ns; ++s) {
961: for (PetscInt q = 0; q < Npc / Ns; ++q, ++p) {
962: species[p] = s;
963: cellid[p] = c;
964: }
965: }
966: }
967: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
968: PetscCall(DMSwarmRestoreField(sw, cellidName, NULL, NULL, (void **)&cellid));
969: }
971: // Set particle coordinates
972: {
973: PetscReal *x, *v;
974: const char **coordNames;
975: PetscInt Ncoord;
976: const PetscInt xdim = dim / 2, vdim = dim / 2;
978: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncoord, &coordNames));
979: PetscCheck(Ncoord == 2, comm, PETSC_ERR_SUP, "We only support regular layout for 2 coordinate fields, not %" PetscInt_FMT, Ncoord);
980: PetscCall(DMSwarmGetField(sw, coordNames[0], NULL, NULL, (void **)&x));
981: PetscCall(DMSwarmGetField(sw, coordNames[1], NULL, NULL, (void **)&v));
982: PetscCall(DMSwarmSortGetAccess(sw));
983: PetscCall(DMGetCoordinatesLocalSetUp(dm));
984: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
985: const PetscInt cell = c + cStart;
986: const PetscScalar *a;
987: PetscScalar *coords;
988: PetscReal lower[6], upper[6];
989: PetscBool isDG;
990: PetscInt *pidx, npc, Nc;
992: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &npc, &pidx));
993: PetscCheck(Npc == npc, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of points per cell %" PetscInt_FMT " != %" PetscInt_FMT, npc, Npc);
994: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &Nc, &a, &coords));
995: for (PetscInt d = 0; d < dim; ++d) {
996: lower[d] = PetscRealPart(coords[0 * dim + d]);
997: upper[d] = PetscRealPart(coords[0 * dim + d]);
998: }
999: for (PetscInt i = 1; i < Nc / dim; ++i) {
1000: for (PetscInt d = 0; d < dim; ++d) {
1001: lower[d] = PetscMin(lower[d], PetscRealPart(coords[i * dim + d]));
1002: upper[d] = PetscMax(upper[d], PetscRealPart(coords[i * dim + d]));
1003: }
1004: }
1005: for (PetscInt s = 0; s < Ns; ++s) {
1006: for (PetscInt q = 0; q < Npc / Ns; ++q) {
1007: const PetscInt p = pidx[q * Ns + s];
1008: PetscInt xi[3], vi[3];
1010: xi[0] = q % n;
1011: xi[1] = (q / n) % n;
1012: xi[2] = (q / PetscSqr(n)) % n;
1013: for (PetscInt d = 0; d < xdim; ++d) x[p * xdim + d] = lower[d] + (xi[d] + 0.5) * (upper[d] - lower[d]) / n;
1014: vi[0] = (q / PetscPowInt(n, xdim)) % n;
1015: vi[1] = (q / PetscPowInt(n, xdim + 1)) % n;
1016: vi[2] = (q / PetscPowInt(n, xdim + 2));
1017: for (PetscInt d = 0; d < vdim; ++d) v[p * vdim + d] = lower[xdim + d] + (vi[d] + 0.5) * (upper[xdim + d] - lower[xdim + d]) / n;
1018: if (debug > 1) {
1019: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Particle %4" PetscInt_FMT " ", p));
1020: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x: ("));
1021: for (PetscInt d = 0; d < xdim; ++d) {
1022: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1023: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(x[p * xdim + d])));
1024: }
1025: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") v:("));
1026: for (PetscInt d = 0; d < vdim; ++d) {
1027: if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
1028: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(v[p * vdim + d])));
1029: }
1030: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
1031: }
1032: }
1033: }
1034: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &Nc, &a, &coords));
1035: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1036: }
1037: PetscCall(DMSwarmSortRestoreAccess(sw));
1038: PetscCall(DMSwarmRestoreField(sw, coordNames[0], NULL, NULL, (void **)&x));
1039: PetscCall(DMSwarmRestoreField(sw, coordNames[1], NULL, NULL, (void **)&v));
1040: }
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: /*
1045: @article{MyersColellaVanStraalen2017,
1046: title = {A 4th-order particle-in-cell method with phase-space remapping for the {Vlasov-Poisson} equation},
1047: author = {Andrew Myers and Phillip Colella and Brian Van Straalen},
1048: journal = {SIAM Journal on Scientific Computing},
1049: volume = {39},
1050: issue = {3},
1051: pages = {B467-B485},
1052: doi = {10.1137/16M105962X},
1053: issn = {10957197},
1054: year = {2017},
1055: }
1056: */
1057: static PetscErrorCode W_3_Interpolation_Private(PetscReal x, PetscReal *w)
1058: {
1059: const PetscReal ax = PetscAbsReal(x);
1061: PetscFunctionBegin;
1062: *w = 0.;
1063: // W_3(x) = 1 - 5/2 |x|^2 + 3/2 |x|^3 0 \le |x| \e 1
1064: if (ax <= 1.) *w = 1. - 2.5 * PetscSqr(ax) + 1.5 * PetscSqr(ax) * ax;
1065: // 1/2 (2 - |x|)^2 (1 - |x|) 1 \le |x| \le 2
1066: else if (ax <= 2.) *w = 0.5 * PetscSqr(2. - ax) * (1. - ax);
1067: //PetscCall(PetscPrintf(PETSC_COMM_SELF, " W_3 %g --> %g\n", x, *w));
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }
1071: // Right now, we will assume that the spatial and velocity grids are regular, which will speed up point location immensely
1072: static PetscErrorCode DMSwarmRemap_Colella_Internal(DM sw, DM *rsw)
1073: {
1074: DM xdm, vdm;
1075: DMSwarmCellDM celldm;
1076: PetscReal xmin[3], xmax[3], vmin[3], vmax[3];
1077: PetscInt xend[3], vend[3];
1078: PetscReal *x, *v, *w, *rw;
1079: PetscReal hx[3], hv[3];
1080: PetscInt dim, xcdim, vcdim, xcStart, xcEnd, vcStart, vcEnd, Np, Nfc;
1081: PetscInt debug = ((DM_Swarm *)sw->data)->printWeights;
1082: const char **coordFields;
1084: PetscFunctionBegin;
1085: PetscCall(DMGetDimension(sw, &dim));
1086: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1087: PetscCall(DMGetCoordinateDim(xdm, &xcdim));
1088: // Create a new centroid swarm without weights
1089: PetscCall(DMSwarmDuplicate(sw, rsw));
1090: PetscCall(DMSwarmGetCellDMActive(*rsw, &celldm));
1091: PetscCall(DMSwarmSetCellDMActive(*rsw, "remap"));
1092: PetscCall(InitializeParticles_Regular(*rsw, 1));
1093: PetscCall(DMSwarmSetCellDMActive(*rsw, ((PetscObject)celldm)->name));
1094: PetscCall(DMSwarmGetLocalSize(*rsw, &Np));
1095: // Assume quad mesh and calculate cell diameters (TODO this could be more robust)
1096: {
1097: const PetscScalar *array;
1098: PetscScalar *coords;
1099: PetscBool isDG;
1100: PetscInt Nc;
1102: PetscCall(DMGetBoundingBox(xdm, xmin, xmax));
1103: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1104: PetscCall(DMPlexGetCellCoordinates(xdm, xcStart, &isDG, &Nc, &array, &coords));
1105: hx[0] = PetscRealPart(coords[1 * xcdim + 0] - coords[0 * xcdim + 0]);
1106: hx[1] = xcdim > 1 ? PetscRealPart(coords[2 * xcdim + 1] - coords[1 * xcdim + 1]) : 1.;
1107: PetscCall(DMPlexRestoreCellCoordinates(xdm, xcStart, &isDG, &Nc, &array, &coords));
1108: PetscCall(PetscObjectQuery((PetscObject)sw, "__vdm__", (PetscObject *)&vdm));
1109: PetscCall(DMGetCoordinateDim(vdm, &vcdim));
1110: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
1111: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1112: PetscCall(DMPlexGetCellCoordinates(vdm, vcStart, &isDG, &Nc, &array, &coords));
1113: hv[0] = PetscRealPart(coords[1 * vcdim + 0] - coords[0 * vcdim + 0]);
1114: hv[1] = vcdim > 1 ? PetscRealPart(coords[2 * vcdim + 1] - coords[1 * vcdim + 1]) : 1.;
1115: PetscCall(DMPlexRestoreCellCoordinates(vdm, vcStart, &isDG, &Nc, &array, &coords));
1117: PetscCheck(dim == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Only support 1D distributions at this time");
1118: xend[0] = xcEnd - xcStart;
1119: xend[1] = 1;
1120: vend[0] = vcEnd - vcStart;
1121: vend[1] = 1;
1122: if (debug > 1)
1123: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Phase Grid (%g, %g, %g, %g) (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", (double)PetscRealPart(hx[0]), (double)PetscRealPart(hx[1]), (double)PetscRealPart(hv[0]), (double)PetscRealPart(hv[1]), xend[0], xend[1], vend[0], vend[1]));
1124: }
1125: // Iterate over particles in the original swarm
1126: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1127: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1128: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
1129: PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&x));
1130: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1131: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
1132: PetscCall(DMSwarmGetField(*rsw, "w_q", NULL, NULL, (void **)&rw));
1133: PetscCall(DMSwarmSortGetAccess(sw));
1134: PetscCall(DMSwarmSortGetAccess(*rsw));
1135: PetscCall(DMGetBoundingBox(vdm, vmin, vmax));
1136: PetscCall(DMGetCoordinatesLocalSetUp(xdm));
1137: for (PetscInt i = 0; i < Np; ++i) rw[i] = 0.;
1138: for (PetscInt c = 0; c < xcEnd - xcStart; ++c) {
1139: PetscInt *pidx, Npc;
1140: PetscInt *rpidx, rNpc;
1142: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1143: for (PetscInt q = 0; q < Npc; ++q) {
1144: const PetscInt p = pidx[q];
1145: const PetscReal wp = w[p];
1146: PetscReal Wx[3], Wv[3];
1147: PetscInt xs[3], vs[3];
1149: // Determine the containing cell
1150: for (PetscInt d = 0; d < dim; ++d) {
1151: const PetscReal xp = x[p * dim + d];
1152: const PetscReal vp = v[p * dim + d];
1154: xs[d] = PetscFloorReal((xp - xmin[d]) / hx[d]);
1155: vs[d] = PetscFloorReal((vp - vmin[d]) / hv[d]);
1156: }
1157: // Loop over all grid points within 2 spacings of the particle
1158: if (debug > 2) {
1159: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Interpolating particle %" PetscInt_FMT " wt %g (%g, %g, %g, %g) (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", p, (double)wp, (double)PetscRealPart(x[p * dim + 0]), xcdim > 1 ? (double)PetscRealPart(x[p * xcdim + 1]) : 0., (double)PetscRealPart(v[p * vcdim + 0]), vcdim > 1 ? (double)PetscRealPart(v[p * vcdim + 1]) : 0., xs[0], xs[1], vs[0], vs[1]));
1160: }
1161: for (PetscInt xi = xs[0] - 1; xi < xs[0] + 3; ++xi) {
1162: // Treat xi as periodic
1163: const PetscInt xip = xi < 0 ? xi + xend[0] : (xi >= xend[0] ? xi - xend[0] : xi);
1164: PetscCall(W_3_Interpolation_Private((xmin[0] + (xi + 0.5) * hx[0] - x[p * dim + 0]) / hx[0], &Wx[0]));
1165: for (PetscInt xj = PetscMax(xs[1] - 1, 0); xj < PetscMin(xs[1] + 3, xend[1]); ++xj) {
1166: if (xcdim > 1) PetscCall(W_3_Interpolation_Private((xmin[1] + (xj + 0.5) * hx[1] - x[p * dim + 1]) / hx[1], &Wx[1]));
1167: else Wx[1] = 1.;
1168: for (PetscInt vi = PetscMax(vs[0] - 1, 0); vi < PetscMin(vs[0] + 3, vend[0]); ++vi) {
1169: PetscCall(W_3_Interpolation_Private((vmin[0] + (vi + 0.5) * hv[0] - v[p * dim + 0]) / hv[0], &Wv[0]));
1170: for (PetscInt vj = PetscMax(vs[1] - 1, 0); vj < PetscMin(vs[1] + 3, vend[1]); ++vj) {
1171: const PetscInt rc = xip * xend[1] + xj;
1172: const PetscInt rv = vi * vend[1] + vj;
1174: PetscCall(DMSwarmSortGetPointsPerCell(*rsw, rc, &rNpc, &rpidx));
1175: if (vcdim > 1) PetscCall(W_3_Interpolation_Private((vmin[1] + (vj + 0.5) * hv[1] - v[p * dim + 1]) / hv[1], &Wv[1]));
1176: else Wv[1] = 1.;
1177: if (debug > 2)
1178: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Depositing on particle (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") w = %g (%g, %g, %g, %g)\n", xi, xj, vi, vj, (double)(wp * Wx[0] * Wx[1] * Wv[0] * Wv[1]), (double)Wx[0], (double)Wx[1], (double)Wv[0], (double)Wv[1]));
1179: // Add weight to new particles from original particle using interpolation function
1180: PetscCheck(rNpc == vend[0] * vend[1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid particle velocity binning");
1181: const PetscInt rp = rpidx[rv];
1182: PetscCheck(rp >= 0 && rp < Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Particle index %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", rp, Np);
1183: rw[rp] += wp * Wx[0] * Wx[1] * Wv[0] * Wv[1];
1184: if (debug > 2) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Adding weight %g (%g) to particle %" PetscInt_FMT "\n", (double)(wp * Wx[0] * Wx[1] * Wv[0] * Wv[1]), (double)PetscRealPart(rw[rp]), rp));
1185: PetscCall(DMSwarmSortRestorePointsPerCell(*rsw, rc, &rNpc, &rpidx));
1186: }
1187: }
1188: }
1189: }
1190: }
1191: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1192: }
1193: PetscCall(DMSwarmSortRestoreAccess(sw));
1194: PetscCall(DMSwarmSortRestoreAccess(*rsw));
1195: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&x));
1196: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1197: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
1198: PetscCall(DMSwarmRestoreField(*rsw, "w_q", NULL, NULL, (void **)&rw));
1200: if (debug) {
1201: Vec w;
1203: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, coordFields[0], &w));
1204: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1205: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, coordFields[0], &w));
1206: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &w));
1207: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1208: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &w));
1209: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &w));
1210: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1211: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &w));
1212: PetscCall(DMSwarmCreateGlobalVectorFromField(*rsw, coordFields[0], &w));
1213: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1214: PetscCall(DMSwarmDestroyGlobalVectorFromField(*rsw, coordFields[0], &w));
1215: PetscCall(DMSwarmCreateGlobalVectorFromField(*rsw, "velocity", &w));
1216: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1217: PetscCall(DMSwarmDestroyGlobalVectorFromField(*rsw, "velocity", &w));
1218: PetscCall(DMSwarmCreateGlobalVectorFromField(*rsw, "w_q", &w));
1219: PetscCall(VecViewFromOptions(w, NULL, "-remap_view"));
1220: PetscCall(DMSwarmDestroyGlobalVectorFromField(*rsw, "w_q", &w));
1221: }
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: static void f0_v2(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[])
1226: {
1227: f0[0] = 0.0;
1228: for (PetscInt d = dim / 2; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
1229: }
1231: static PetscErrorCode DMSwarmRemap_PFAK_Internal(DM sw, DM *rsw)
1232: {
1233: DM xdm, vdm, rdm;
1234: DMSwarmCellDM rcelldm;
1235: Mat M_p, rM_p, rPM_p;
1236: Vec w, rw, rhs;
1237: PetscInt Nf;
1238: const char **fields;
1239: PetscBool debug = PETSC_FALSE;
1241: PetscFunctionBegin;
1242: // Create a new centroid swarm without weights
1243: PetscCall(DMSwarmGetCellDM(sw, &xdm));
1244: PetscCall(DMSwarmSetCellDMActive(sw, "velocity"));
1245: PetscCall(DMSwarmGetCellDMActive(sw, &rcelldm));
1246: PetscCall(DMSwarmCellDMGetDM(rcelldm, &vdm));
1247: PetscCall(DMSwarmDuplicate(sw, rsw));
1248: // Set remap cell DM
1249: PetscCall(DMSwarmSetCellDMActive(sw, "remap"));
1250: PetscCall(DMSwarmGetCellDMActive(sw, &rcelldm));
1251: PetscCall(DMSwarmCellDMGetFields(rcelldm, &Nf, &fields));
1252: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "We only allow a single weight field, not %" PetscInt_FMT, Nf);
1253: PetscCall(DMSwarmGetCellDM(sw, &rdm));
1254: PetscCall(DMGetGlobalVector(rdm, &rhs));
1255: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE)); // Bin particles in remap mesh
1256: // Compute rhs = M_p w_p
1257: PetscCall(DMCreateMassMatrix(sw, rdm, &M_p));
1258: if (debug) {
1259: Vec col;
1260: PetscInt m, n;
1261: PetscBool rankDeficient = PETSC_FALSE;
1263: PetscCall(MatGetSize(M_p, &m, &n));
1264: PetscCall(MatCreateVecs(M_p, NULL, &col));
1265: for (PetscInt c = 0; c < n; ++c) {
1266: const PetscScalar *a;
1267: PetscInt num = 0;
1269: PetscCall(MatGetColumnVector(M_p, col, c));
1270: PetscCall(VecGetArrayRead(col, &a));
1271: for (PetscInt r = 0; r < m; ++r)
1272: if (a[r] != 0.) ++num;
1273: PetscCall(VecRestoreArrayRead(col, &a));
1274: if (num < 2) {
1275: rankDeficient = PETSC_TRUE;
1276: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Basis function %" PetscInt_FMT " has only %" PetscInt_FMT " particles in its support\n", c, num));
1277: }
1278: }
1279: PetscCall(VecDestroy(&col));
1280: if (rankDeficient) PetscCall(DMViewFromOptions(sw, NULL, "-rank_def_view"));
1281: }
1282: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fields[0], &w));
1283: PetscCall(VecViewFromOptions(w, NULL, "-remap_w_view"));
1284: PetscCall(MatMultTranspose(M_p, w, rhs));
1285: PetscCall(VecViewFromOptions(rhs, NULL, "-remap_rhs_view"));
1286: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fields[0], &w));
1287: PetscCall(MatDestroy(&M_p));
1288: {
1289: KSP ksp;
1290: Mat M_f;
1291: Vec u_f;
1292: PetscReal mom[4];
1293: PetscInt cdim;
1294: const char *prefix;
1296: PetscCall(DMGetCoordinateDim(rdm, &cdim));
1297: PetscCall(DMCreateMassMatrix(rdm, rdm, &M_f));
1298: PetscCall(DMGetGlobalVector(rdm, &u_f));
1300: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
1301: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1302: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1303: PetscCall(KSPAppendOptionsPrefix(ksp, "ptof_"));
1304: PetscCall(KSPSetFromOptions(ksp));
1306: PetscCall(KSPSetOperators(ksp, M_f, M_f));
1307: PetscCall(KSPSolve(ksp, rhs, u_f));
1308: PetscCall(KSPDestroy(&ksp));
1309: PetscCall(VecViewFromOptions(u_f, NULL, "-remap_uf_view"));
1311: PetscCall(DMPlexComputeMoments(rdm, u_f, mom));
1312: // Energy is not correct since it uses (x^2 + v^2)
1313: PetscDS rds;
1314: PetscScalar rmom;
1315: void *ctx;
1317: PetscCall(DMGetDS(rdm, &rds));
1318: PetscCall(DMGetApplicationContext(rdm, &ctx));
1319: PetscCall(PetscDSSetObjective(rds, 0, &f0_v2));
1320: PetscCall(DMPlexComputeIntegralFEM(rdm, u_f, &rmom, ctx));
1321: mom[1 + cdim] = PetscRealPart(rmom);
1323: PetscCall(DMRestoreGlobalVector(rdm, &u_f));
1324: PetscCall(PetscPrintf(PETSC_COMM_SELF, "========== PFAK u_f ==========\n"));
1325: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom 0: %g\n", (double)mom[0]));
1326: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom x: %g\n", (double)mom[1 + 0]));
1327: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom v: %g\n", (double)mom[1 + 1]));
1328: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom 2: %g\n", (double)mom[1 + cdim]));
1329: PetscCall(MatDestroy(&M_f));
1330: }
1331: // Create Remap particle mass matrix M_p
1332: PetscInt xcStart, xcEnd, vcStart, vcEnd, cStart, cEnd, r;
1334: PetscCall(DMSwarmSetCellDMActive(*rsw, "remap"));
1335: PetscCall(DMPlexGetHeightStratum(xdm, 0, &xcStart, &xcEnd));
1336: PetscCall(DMPlexGetHeightStratum(vdm, 0, &vcStart, &vcEnd));
1337: PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStart, &cEnd));
1338: r = (PetscInt)PetscSqrtReal(((xcEnd - xcStart) * (vcEnd - vcStart)) / (cEnd - cStart));
1339: PetscCall(InitializeParticles_Regular(*rsw, r));
1340: PetscCall(DMSwarmMigrate(*rsw, PETSC_FALSE)); // Bin particles in remap mesh
1341: PetscCall(DMCreateMassMatrix(*rsw, rdm, &rM_p));
1342: PetscCall(MatViewFromOptions(rM_p, NULL, "-rM_p_view"));
1343: // Solve M_p
1344: {
1345: KSP ksp;
1346: PC pc;
1347: const char *prefix;
1348: PetscBool isBjacobi;
1350: PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
1351: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1352: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1353: PetscCall(KSPAppendOptionsPrefix(ksp, "ftop_"));
1354: PetscCall(KSPSetFromOptions(ksp));
1356: PetscCall(KSPGetPC(ksp, &pc));
1357: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &isBjacobi));
1358: if (isBjacobi) {
1359: PetscCall(DMSwarmCreateMassMatrixSquare(sw, rdm, &rPM_p));
1360: } else {
1361: rPM_p = rM_p;
1362: PetscCall(PetscObjectReference((PetscObject)rPM_p));
1363: }
1364: PetscCall(KSPSetOperators(ksp, rM_p, rPM_p));
1365: PetscCall(DMSwarmCreateGlobalVectorFromField(*rsw, fields[0], &rw));
1366: PetscCall(KSPSolveTranspose(ksp, rhs, rw));
1367: PetscCall(VecViewFromOptions(rw, NULL, "-remap_rw_view"));
1368: PetscCall(DMSwarmDestroyGlobalVectorFromField(*rsw, fields[0], &rw));
1369: PetscCall(KSPDestroy(&ksp));
1370: PetscCall(MatDestroy(&rPM_p));
1371: PetscCall(MatDestroy(&rM_p));
1372: }
1373: PetscCall(DMRestoreGlobalVector(rdm, &rhs));
1375: // Restore original cell DM
1376: PetscCall(DMSwarmSetCellDMActive(sw, "space"));
1377: PetscCall(DMSwarmSetCellDMActive(*rsw, "space"));
1378: PetscCall(DMSwarmMigrate(*rsw, PETSC_FALSE)); // Bin particles in spatial mesh
1379: PetscFunctionReturn(PETSC_SUCCESS);
1380: }
1382: static PetscErrorCode DMSwarmRemapMonitor_Internal(DM sw, DM rsw)
1383: {
1384: PetscReal mom[4], rmom[4];
1385: PetscInt cdim;
1387: PetscFunctionBegin;
1388: PetscCall(DMGetCoordinateDim(sw, &cdim));
1389: PetscCall(DMSwarmComputeMoments(sw, "velocity", "w_q", mom));
1390: PetscCall(DMSwarmComputeMoments(rsw, "velocity", "w_q", rmom));
1391: PetscCall(PetscPrintf(PETSC_COMM_SELF, "========== Remapped ==========\n"));
1392: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom 0: %g --> %g\n", (double)mom[0], (double)rmom[0]));
1393: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom 1: %g --> %g\n", (double)mom[1], (double)rmom[1]));
1394: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Mom 2: %g --> %g\n", (double)mom[1 + cdim], (double)rmom[1 + cdim]));
1395: PetscFunctionReturn(PETSC_SUCCESS);
1396: }
1398: /*@
1399: DMSwarmRemap - Project the swarm fields onto a new set of particles
1401: Collective
1403: Input Parameter:
1404: . sw - The `DMSWARM` object
1406: Level: beginner
1408: .seealso: [](ch_dmbase), `DMSWARM`, `DMSwarmMigrate()`, `DMSwarmCrate()`
1409: @*/
1410: PetscErrorCode DMSwarmRemap(DM sw)
1411: {
1412: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1413: DM rsw;
1415: PetscFunctionBegin;
1416: switch (swarm->remap_type) {
1417: case DMSWARM_REMAP_NONE:
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: case DMSWARM_REMAP_COLELLA:
1420: PetscCall(DMSwarmRemap_Colella_Internal(sw, &rsw));
1421: break;
1422: case DMSWARM_REMAP_PFAK:
1423: PetscCall(DMSwarmRemap_PFAK_Internal(sw, &rsw));
1424: break;
1425: default:
1426: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No remap algorithm %s", DMSwarmRemapTypeNames[swarm->remap_type]);
1427: }
1428: PetscCall(DMSwarmRemapMonitor_Internal(sw, rsw));
1429: PetscCall(DMSwarmReplace(sw, &rsw));
1430: PetscFunctionReturn(PETSC_SUCCESS);
1431: }