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