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: }