Actual source code: dmproject.c
2: #include <petsc/private/dmimpl.h>
3: #include <petscdm.h>
4: #include <petscdmplex.h>
5: #include <petscksp.h>
6: #include <petscblaslapack.h>
8: typedef struct _projectConstraintsCtx {
9: DM dm;
10: Vec mask;
11: } projectConstraintsCtx;
13: PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y)
14: {
15: DM dm;
16: Vec local, mask;
17: projectConstraintsCtx *ctx;
19: MatShellGetContext(CtC, &ctx);
20: dm = ctx->dm;
21: mask = ctx->mask;
22: DMGetLocalVector(dm, &local);
23: DMGlobalToLocalBegin(dm, x, INSERT_VALUES, local);
24: DMGlobalToLocalEnd(dm, x, INSERT_VALUES, local);
25: if (mask) VecPointwiseMult(local, mask, local);
26: VecSet(y, 0.);
27: DMLocalToGlobalBegin(dm, local, ADD_VALUES, y);
28: DMLocalToGlobalEnd(dm, local, ADD_VALUES, y);
29: DMRestoreLocalVector(dm, &local);
30: return 0;
31: }
33: static PetscErrorCode DMGlobalToLocalSolve_project1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
34: {
35: PetscInt f;
37: for (f = 0; f < Nf; f++) u[f] = 1.;
38: return 0;
39: }
41: /*@
42: DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by `DMGlobalToLocalBegin()`/`DMGlobalToLocalEnd()` with mode
43: = `INSERT_VALUES`. 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
44: a least-squares solution. It is also assumed that `DMLocalToGlobalBegin()`/`DMLocalToGlobalEnd()` with mode = `ADD_VALUES` is the adjoint of the
45: global-to-local map, so that the least-squares solution may be found by the normal equations.
47: collective
49: Input Parameters:
50: + dm - The `DM` object
51: . x - The local vector
52: - y - The global vector: the input value of globalVec is used as an initial guess
54: Output Parameters:
55: . y - The least-squares solution
57: Level: advanced
59: Note:
60: If the `DM` is of type `DMPLEX`, then y is the solution of L' * D * L * y = L' * D * x, where D is a diagonal mask that is 1 for every point in
61: 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
62: closure of any local cell (see `DMPlexGetAnchors()`/`DMPlexSetAnchors()`).
64: .seealso: [](chapter_ksp), `DM`, `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToGlobalEnd()`, `DMPlexGetAnchors()`, `DMPlexSetAnchors()`
65: @*/
66: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
67: {
68: Mat CtC;
69: PetscInt n, N, cStart, cEnd, c;
70: PetscBool isPlex;
71: KSP ksp;
72: PC pc;
73: Vec global, mask = NULL;
74: projectConstraintsCtx ctx;
76: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
77: if (isPlex) {
78: /* mark points in the closure */
79: DMCreateLocalVector(dm, &mask);
80: VecSet(mask, 0.0);
81: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
82: if (cEnd > cStart) {
83: PetscScalar *ones;
84: PetscInt numValues, i;
86: DMPlexVecGetClosure(dm, NULL, mask, cStart, &numValues, NULL);
87: PetscMalloc1(numValues, &ones);
88: for (i = 0; i < numValues; i++) ones[i] = 1.;
89: for (c = cStart; c < cEnd; c++) DMPlexVecSetClosure(dm, NULL, mask, c, ones, INSERT_VALUES);
90: PetscFree(ones);
91: }
92: } else {
93: PetscBool hasMask;
95: DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask);
96: if (!hasMask) {
97: PetscErrorCode (**func)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
98: void **ctx;
99: PetscInt Nf, f;
101: DMGetNumFields(dm, &Nf);
102: PetscMalloc2(Nf, &func, Nf, &ctx);
103: for (f = 0; f < Nf; ++f) {
104: func[f] = DMGlobalToLocalSolve_project1;
105: ctx[f] = NULL;
106: }
107: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
108: DMProjectFunctionLocal(dm, 0.0, func, ctx, INSERT_ALL_VALUES, mask);
109: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
110: PetscFree2(func, ctx);
111: }
112: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
113: }
114: ctx.dm = dm;
115: ctx.mask = mask;
116: VecGetSize(y, &N);
117: VecGetLocalSize(y, &n);
118: MatCreate(PetscObjectComm((PetscObject)dm), &CtC);
119: MatSetSizes(CtC, n, n, N, N);
120: MatSetType(CtC, MATSHELL);
121: MatSetUp(CtC);
122: MatShellSetContext(CtC, &ctx);
123: MatShellSetOperation(CtC, MATOP_MULT, (void (*)(void))MatMult_GlobalToLocalNormal);
124: KSPCreate(PetscObjectComm((PetscObject)dm), &ksp);
125: KSPSetOperators(ksp, CtC, CtC);
126: KSPSetType(ksp, KSPCG);
127: KSPGetPC(ksp, &pc);
128: PCSetType(pc, PCNONE);
129: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
130: KSPSetUp(ksp);
131: DMGetGlobalVector(dm, &global);
132: VecSet(global, 0.);
133: if (mask) VecPointwiseMult(x, mask, x);
134: DMLocalToGlobalBegin(dm, x, ADD_VALUES, global);
135: DMLocalToGlobalEnd(dm, x, ADD_VALUES, global);
136: KSPSolve(ksp, global, y);
137: DMRestoreGlobalVector(dm, &global);
138: /* clean up */
139: KSPDestroy(&ksp);
140: MatDestroy(&CtC);
141: if (isPlex) {
142: VecDestroy(&mask);
143: } else {
144: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
145: }
147: return 0;
148: }
150: /*@C
151: DMProjectField - This projects the given function of the input fields into the function space provided, putting the coefficients in a global vector.
153: Collective on dm
155: Input Parameters:
156: + dm - The `DM`
157: . time - The time
158: . U - The input field vector
159: . funcs - The functions to evaluate, one per field
160: - mode - The insertion mode for values
162: Output Parameter:
163: . X - The output vector
165: Calling sequence of func:
166: .vb
167: func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
168: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
169: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
170: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
171: .ve
173: + dim - The spatial dimension
174: . Nf - The number of input fields
175: . NfAux - The number of input auxiliary fields
176: . uOff - The offset of each field in u[]
177: . uOff_x - The offset of each field in u_x[]
178: . u - The field values at this point in space
179: . u_t - The field time derivative at this point in space (or NULL)
180: . u_x - The field derivatives at this point in space
181: . aOff - The offset of each auxiliary field in u[]
182: . aOff_x - The offset of each auxiliary field in u_x[]
183: . a - The auxiliary field values at this point in space
184: . a_t - The auxiliary field time derivative at this point in space (or NULL)
185: . a_x - The auxiliary field derivatives at this point in space
186: . t - The current time
187: . x - The coordinates of this point
188: . numConstants - The number of constants
189: . constants - The value of each constant
190: - f - The value of the function at this point in space
192: Level: advanced
194: Note:
195: There are three different `DM`s that potentially interact in this function. The output `DM`, dm, specifies the layout of the values calculates by funcs.
196: 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
197: a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary `DM`, attached to the
198: auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.
200: .seealso: [](chapter_ksp), `DM`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`, `DMProjectFunction()`, `DMComputeL2Diff()`
201: @*/
202: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec X)
203: {
204: Vec localX, localU;
205: DM dmIn;
208: DMGetLocalVector(dm, &localX);
209: /* We currently check whether locU == locX to see if we need to apply BC */
210: if (U != X) {
211: VecGetDM(U, &dmIn);
212: DMGetLocalVector(dmIn, &localU);
213: } else {
214: dmIn = dm;
215: localU = localX;
216: }
217: DMGlobalToLocalBegin(dmIn, U, INSERT_VALUES, localU);
218: DMGlobalToLocalEnd(dmIn, U, INSERT_VALUES, localU);
219: DMProjectFieldLocal(dm, time, localU, funcs, mode, localX);
220: DMLocalToGlobalBegin(dm, localX, mode, X);
221: DMLocalToGlobalEnd(dm, localX, mode, X);
222: if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
223: Mat cMat;
225: DMGetDefaultConstraints(dm, NULL, &cMat, NULL);
226: if (cMat) DMGlobalToLocalSolve(dm, localX, X);
227: }
228: DMRestoreLocalVector(dm, &localX);
229: if (U != X) DMRestoreLocalVector(dmIn, &localU);
230: return 0;
231: }
233: /********************* Adaptive Interpolation **************************/
235: /* See the discussion of Adaptive Interpolation in manual/high_level_mg.rst */
236: PetscErrorCode DMAdaptInterpolator(DM dmc, DM dmf, Mat In, KSP smoother, Mat MF, Mat MC, Mat *InAdapt, void *user)
237: {
238: Mat globalA, AF;
239: Vec tmp;
240: const PetscScalar *af, *ac;
241: PetscScalar *A, *b, *x, *workscalar;
242: PetscReal *w, *sing, *workreal, rcond = PETSC_SMALL;
243: PetscBLASInt M, N, one = 1, irank, lwrk, info;
244: PetscInt debug = 0, rStart, rEnd, r, maxcols = 0, k, Nc, ldac, ldaf;
245: PetscBool allocVc = PETSC_FALSE;
247: PetscLogEventBegin(DM_AdaptInterpolator, dmc, dmf, 0, 0);
248: PetscOptionsGetInt(NULL, NULL, "-dm_interpolator_adapt_debug", &debug, NULL);
249: MatGetSize(MF, NULL, &Nc);
250: MatDuplicate(In, MAT_SHARE_NONZERO_PATTERN, InAdapt);
251: MatGetOwnershipRange(In, &rStart, &rEnd);
252: #if 0
253: MatGetMaxRowLen(In, &maxcols);
254: #else
255: for (r = rStart; r < rEnd; ++r) {
256: PetscInt ncols;
258: MatGetRow(In, r, &ncols, NULL, NULL);
259: maxcols = PetscMax(maxcols, ncols);
260: MatRestoreRow(In, r, &ncols, NULL, NULL);
261: }
262: #endif
263: if (Nc < maxcols) PetscPrintf(PETSC_COMM_SELF, "The number of input vectors %" PetscInt_FMT " < %" PetscInt_FMT " the maximum number of column entries\n", Nc, maxcols);
264: for (k = 0; k < Nc && debug; ++k) {
265: char name[PETSC_MAX_PATH_LEN];
266: const char *prefix;
267: Vec vc, vf;
269: PetscObjectGetOptionsPrefix((PetscObject)smoother, &prefix);
271: if (MC) {
272: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sCoarse Vector %" PetscInt_FMT, prefix ? prefix : NULL, k);
273: MatDenseGetColumnVecRead(MC, k, &vc);
274: PetscObjectSetName((PetscObject)vc, name);
275: VecViewFromOptions(vc, NULL, "-dm_adapt_interp_view_coarse");
276: MatDenseRestoreColumnVecRead(MC, k, &vc);
277: }
278: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sFine Vector %" PetscInt_FMT, prefix ? prefix : NULL, k);
279: MatDenseGetColumnVecRead(MF, k, &vf);
280: PetscObjectSetName((PetscObject)vf, name);
281: VecViewFromOptions(vf, NULL, "-dm_adapt_interp_view_fine");
282: MatDenseRestoreColumnVecRead(MF, k, &vf);
283: }
284: PetscBLASIntCast(3 * PetscMin(Nc, maxcols) + PetscMax(2 * PetscMin(Nc, maxcols), PetscMax(Nc, maxcols)), &lwrk);
285: PetscMalloc7(Nc * maxcols, &A, PetscMax(Nc, maxcols), &b, Nc, &w, maxcols, &x, maxcols, &sing, lwrk, &workscalar, 5 * PetscMin(Nc, maxcols), &workreal);
286: /* 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} */
287: KSPGetOperators(smoother, &globalA, NULL);
289: MatMatMult(globalA, MF, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AF);
290: for (k = 0; k < Nc; ++k) {
291: PetscScalar vnorm, vAnorm;
292: Vec vf;
294: w[k] = 1.0;
295: MatDenseGetColumnVecRead(MF, k, &vf);
296: MatDenseGetColumnVecRead(AF, k, &tmp);
297: VecDot(vf, vf, &vnorm);
298: #if 0
299: DMGetGlobalVector(dmf, &tmp2);
300: KSPSolve(smoother, tmp, tmp2);
301: VecDot(vf, tmp2, &vAnorm);
302: DMRestoreGlobalVector(dmf, &tmp2);
303: #else
304: VecDot(vf, tmp, &vAnorm);
305: #endif
306: w[k] = PetscRealPart(vnorm) / PetscRealPart(vAnorm);
307: MatDenseRestoreColumnVecRead(MF, k, &vf);
308: MatDenseRestoreColumnVecRead(AF, k, &tmp);
309: }
310: MatDestroy(&AF);
311: if (!MC) {
312: allocVc = PETSC_TRUE;
313: MatTransposeMatMult(In, MF, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &MC);
314: }
315: /* Solve a LS system for each fine row
316: MATT: Can we generalize to the case where Nc for the fine space
317: is different for Nc for the coarse? */
318: MatDenseGetArrayRead(MF, &af);
319: MatDenseGetLDA(MF, &ldaf);
320: MatDenseGetArrayRead(MC, &ac);
321: MatDenseGetLDA(MC, &ldac);
322: for (r = rStart; r < rEnd; ++r) {
323: PetscInt ncols, c;
324: const PetscInt *cols;
325: const PetscScalar *vals;
327: MatGetRow(In, r, &ncols, &cols, &vals);
328: for (k = 0; k < Nc; ++k) {
329: /* Need to fit lowest mode exactly */
330: const PetscReal wk = ((ncols == 1) && (k > 0)) ? 0.0 : PetscSqrtReal(w[k]);
332: /* b_k = \sqrt{w_k} f^{F,k}_r */
333: b[k] = wk * af[r - rStart + k * ldaf];
334: /* A_{kc} = \sqrt{w_k} f^{C,k}_c */
335: /* TODO Must pull out VecScatter from In, scatter in vc[k] values up front, and access them indirectly just as in MatMult() */
336: for (c = 0; c < ncols; ++c) {
337: /* This is element (k, c) of A */
338: A[c * Nc + k] = wk * ac[cols[c] - rStart + k * ldac];
339: }
340: }
341: PetscBLASIntCast(Nc, &M);
342: PetscBLASIntCast(ncols, &N);
343: if (debug) {
344: #if defined(PETSC_USE_COMPLEX)
345: PetscScalar *tmp;
346: PetscInt j;
348: DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp);
349: for (j = 0; j < Nc; ++j) tmp[j] = w[j];
350: DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, tmp);
351: DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A);
352: for (j = 0; j < Nc; ++j) tmp[j] = b[j];
353: DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, tmp);
354: DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp);
355: #else
356: DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, w);
357: DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A);
358: DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, b);
359: #endif
360: }
361: #if defined(PETSC_USE_COMPLEX)
362: /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
363: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, workreal, &info));
364: #else
365: /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
366: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, &info));
367: #endif
370: if (debug) {
371: PetscPrintf(PETSC_COMM_SELF, "rank %" PetscBLASInt_FMT " rcond %g\n", irank, (double)rcond);
372: #if defined(PETSC_USE_COMPLEX)
373: {
374: PetscScalar *tmp;
375: PetscInt j;
377: DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp);
378: for (j = 0; j < PetscMin(Nc, ncols); ++j) tmp[j] = sing[j];
379: DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, tmp);
380: DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp);
381: }
382: #else
383: DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, sing);
384: #endif
385: DMPrintCellMatrix(r, "Interpolator Row LS old P", ncols, 1, vals);
386: DMPrintCellMatrix(r, "Interpolator Row LS sol", ncols, 1, b);
387: }
388: MatSetValues(*InAdapt, 1, &r, ncols, cols, b, INSERT_VALUES);
389: MatRestoreRow(In, r, &ncols, &cols, &vals);
390: }
391: MatDenseRestoreArrayRead(MF, &af);
392: MatDenseRestoreArrayRead(MC, &ac);
393: PetscFree7(A, b, w, x, sing, workscalar, workreal);
394: if (allocVc) MatDestroy(&MC);
395: MatAssemblyBegin(*InAdapt, MAT_FINAL_ASSEMBLY);
396: MatAssemblyEnd(*InAdapt, MAT_FINAL_ASSEMBLY);
397: PetscLogEventEnd(DM_AdaptInterpolator, dmc, dmf, 0, 0);
398: return 0;
399: }
401: PetscErrorCode DMCheckInterpolator(DM dmf, Mat In, Mat MC, Mat MF, PetscReal tol)
402: {
403: Vec tmp;
404: PetscReal norminf, norm2, maxnorminf = 0.0, maxnorm2 = 0.0;
405: PetscInt k, Nc;
407: DMGetGlobalVector(dmf, &tmp);
408: MatViewFromOptions(In, NULL, "-dm_interpolator_adapt_error");
409: MatGetSize(MF, NULL, &Nc);
410: for (k = 0; k < Nc; ++k) {
411: Vec vc, vf;
413: MatDenseGetColumnVecRead(MC, k, &vc);
414: MatDenseGetColumnVecRead(MF, k, &vf);
415: MatMult(In, vc, tmp);
416: VecAXPY(tmp, -1.0, vf);
417: VecViewFromOptions(vc, NULL, "-dm_interpolator_adapt_error");
418: VecViewFromOptions(vf, NULL, "-dm_interpolator_adapt_error");
419: VecViewFromOptions(tmp, NULL, "-dm_interpolator_adapt_error");
420: VecNorm(tmp, NORM_INFINITY, &norminf);
421: VecNorm(tmp, NORM_2, &norm2);
422: maxnorminf = PetscMax(maxnorminf, norminf);
423: maxnorm2 = PetscMax(maxnorm2, norm2);
424: PetscPrintf(PetscObjectComm((PetscObject)dmf), "Coarse vec %" PetscInt_FMT " ||vf - P vc||_\\infty %g, ||vf - P vc||_2 %g\n", k, (double)norminf, (double)norm2);
425: MatDenseRestoreColumnVecRead(MC, k, &vc);
426: MatDenseRestoreColumnVecRead(MF, k, &vf);
427: }
428: DMRestoreGlobalVector(dmf, &tmp);
430: return 0;
431: }