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