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