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