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:   PetscFunctionBegin;
 20:   PetscCall(MatShellGetContext(CtC, &ctx));
 21:   dm   = ctx->dm;
 22:   mask = ctx->mask;
 23:   PetscCall(DMGetLocalVector(dm, &local));
 24:   PetscCall(DMGlobalToLocalBegin(dm, x, INSERT_VALUES, local));
 25:   PetscCall(DMGlobalToLocalEnd(dm, x, INSERT_VALUES, local));
 26:   if (mask) PetscCall(VecPointwiseMult(local, mask, local));
 27:   PetscCall(VecSet(y, 0.));
 28:   PetscCall(DMLocalToGlobalBegin(dm, local, ADD_VALUES, y));
 29:   PetscCall(DMLocalToGlobalEnd(dm, local, ADD_VALUES, y));
 30:   PetscCall(DMRestoreLocalVector(dm, &local));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: static PetscErrorCode DMGlobalToLocalSolve_project1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
 35: {
 36:   PetscInt f;

 38:   PetscFunctionBegin;
 39:   for (f = 0; f < Nf; f++) u[f] = 1.;
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: /*@
 44:   DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by `DMGlobalToLocalBegin()`/`DMGlobalToLocalEnd()` with mode
 45:   = `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
 46:   a least-squares solution.  It is also assumed that `DMLocalToGlobalBegin()`/`DMLocalToGlobalEnd()` with mode = `ADD_VALUES` is the adjoint of the
 47:   global-to-local map, so that the least-squares solution may be found by the normal equations.

 49:   Collective

 51:   Input Parameters:
 52: + dm - The `DM` object
 53: . x - The local vector
 54: - y - The global vector: the input value of globalVec is used as an initial guess

 56:   Output Parameter:
 57: . y - The least-squares solution

 59:   Level: advanced

 61:   Note:
 62:   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
 63:   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
 64:   closure of any local cell (see `DMPlexGetAnchors()`/`DMPlexSetAnchors()`).

 66: .seealso: [](ch_ksp), `DM`, `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToGlobalEnd()`, `DMPlexGetAnchors()`, `DMPlexSetAnchors()`
 67: @*/
 68: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
 69: {
 70:   Mat                   CtC;
 71:   PetscInt              n, N, cStart, cEnd, c;
 72:   PetscBool             isPlex;
 73:   KSP                   ksp;
 74:   PC                    pc;
 75:   Vec                   global, mask = NULL;
 76:   projectConstraintsCtx ctx;

 78:   PetscFunctionBegin;
 79:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
 80:   if (isPlex) {
 81:     /* mark points in the closure */
 82:     PetscCall(DMCreateLocalVector(dm, &mask));
 83:     PetscCall(VecSet(mask, 0.0));
 84:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
 85:     if (cEnd > cStart) {
 86:       PetscScalar *ones;
 87:       PetscInt     numValues, i;

 89:       PetscCall(DMPlexVecGetClosure(dm, NULL, mask, cStart, &numValues, NULL));
 90:       PetscCall(PetscMalloc1(numValues, &ones));
 91:       for (i = 0; i < numValues; i++) ones[i] = 1.;
 92:       for (c = cStart; c < cEnd; c++) PetscCall(DMPlexVecSetClosure(dm, NULL, mask, c, ones, INSERT_VALUES));
 93:       PetscCall(PetscFree(ones));
 94:     }
 95:   } else {
 96:     PetscBool hasMask;

 98:     PetscCall(DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask));
 99:     if (!hasMask) {
100:       PetscErrorCode (**func)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
101:       void   **ctx;
102:       PetscInt Nf, f;

104:       PetscCall(DMGetNumFields(dm, &Nf));
105:       PetscCall(PetscMalloc2(Nf, &func, Nf, &ctx));
106:       for (f = 0; f < Nf; ++f) {
107:         func[f] = DMGlobalToLocalSolve_project1;
108:         ctx[f]  = NULL;
109:       }
110:       PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
111:       PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctx, INSERT_ALL_VALUES, mask));
112:       PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
113:       PetscCall(PetscFree2(func, ctx));
114:     }
115:     PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
116:   }
117:   ctx.dm   = dm;
118:   ctx.mask = mask;
119:   PetscCall(VecGetSize(y, &N));
120:   PetscCall(VecGetLocalSize(y, &n));
121:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &CtC));
122:   PetscCall(MatSetSizes(CtC, n, n, N, N));
123:   PetscCall(MatSetType(CtC, MATSHELL));
124:   PetscCall(MatSetUp(CtC));
125:   PetscCall(MatShellSetContext(CtC, &ctx));
126:   PetscCall(MatShellSetOperation(CtC, MATOP_MULT, (void (*)(void))MatMult_GlobalToLocalNormal));
127:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
128:   PetscCall(KSPSetOperators(ksp, CtC, CtC));
129:   PetscCall(KSPSetType(ksp, KSPCG));
130:   PetscCall(KSPGetPC(ksp, &pc));
131:   PetscCall(PCSetType(pc, PCNONE));
132:   PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
133:   PetscCall(KSPSetUp(ksp));
134:   PetscCall(DMGetGlobalVector(dm, &global));
135:   PetscCall(VecSet(global, 0.));
136:   if (mask) PetscCall(VecPointwiseMult(x, mask, x));
137:   PetscCall(DMLocalToGlobalBegin(dm, x, ADD_VALUES, global));
138:   PetscCall(DMLocalToGlobalEnd(dm, x, ADD_VALUES, global));
139:   PetscCall(KSPSolve(ksp, global, y));
140:   PetscCall(DMRestoreGlobalVector(dm, &global));
141:   /* clean up */
142:   PetscCall(KSPDestroy(&ksp));
143:   PetscCall(MatDestroy(&CtC));
144:   if (isPlex) {
145:     PetscCall(VecDestroy(&mask));
146:   } else {
147:     PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
148:   }

150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: /*@C
154:   DMProjectField - This projects the given function of the input fields into the function space provided, putting the coefficients in a global vector.

156:   Collective

158:   Input Parameters:
159: + dm      - The `DM`
160: . time    - The time
161: . U       - The input field vector
162: . funcs   - The functions to evaluate, one per field
163: - mode    - The insertion mode for values

165:   Output Parameter:
166: . X       - The output vector

168:    Calling sequence of `func`:
169: .vb
170:    void funcs(PetscInt dim, PetscInt Nf, PetscInt NfAux,
171:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
172:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
173:               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
174: .ve
175: +  dim          - The spatial dimension
176: .  Nf           - The number of input fields
177: .  NfAux        - The number of input auxiliary fields
178: .  uOff         - The offset of each field in u[]
179: .  uOff_x       - The offset of each field in u_x[]
180: .  u            - The field values at this point in space
181: .  u_t          - The field time derivative at this point in space (or `NULL`)
182: .  u_x          - The field derivatives at this point in space
183: .  aOff         - The offset of each auxiliary field in u[]
184: .  aOff_x       - The offset of each auxiliary field in u_x[]
185: .  a            - The auxiliary field values at this point in space
186: .  a_t          - The auxiliary field time derivative at this point in space (or `NULL`)
187: .  a_x          - The auxiliary field derivatives at this point in space
188: .  t            - The current time
189: .  x            - The coordinates of this point
190: .  numConstants - The number of constants
191: .  constants    - The value of each constant
192: -  f            - The value of the function at this point in space

194:   Level: advanced

196:   Note:
197:   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.
198:   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
199:   a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary `DM`, attached to the
200:   auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.

202: .seealso: [](ch_ksp), `DM`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`, `DMProjectFunction()`, `DMComputeL2Diff()`
203: @*/
204: 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)
205: {
206:   Vec localX, localU;
207:   DM  dmIn;

209:   PetscFunctionBegin;
211:   PetscCall(DMGetLocalVector(dm, &localX));
212:   /* We currently check whether locU == locX to see if we need to apply BC */
213:   if (U != X) {
214:     PetscCall(VecGetDM(U, &dmIn));
215:     PetscCall(DMGetLocalVector(dmIn, &localU));
216:   } else {
217:     dmIn   = dm;
218:     localU = localX;
219:   }
220:   PetscCall(DMGlobalToLocalBegin(dmIn, U, INSERT_VALUES, localU));
221:   PetscCall(DMGlobalToLocalEnd(dmIn, U, INSERT_VALUES, localU));
222:   PetscCall(DMProjectFieldLocal(dm, time, localU, funcs, mode, localX));
223:   PetscCall(DMLocalToGlobalBegin(dm, localX, mode, X));
224:   PetscCall(DMLocalToGlobalEnd(dm, localX, mode, X));
225:   if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
226:     Mat cMat;

228:     PetscCall(DMGetDefaultConstraints(dm, NULL, &cMat, NULL));
229:     if (cMat) PetscCall(DMGlobalToLocalSolve(dm, localX, X));
230:   }
231:   PetscCall(DMRestoreLocalVector(dm, &localX));
232:   if (U != X) PetscCall(DMRestoreLocalVector(dmIn, &localU));
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: /********************* Adaptive Interpolation **************************/

238: /* See the discussion of Adaptive Interpolation in manual/high_level_mg.rst */
239: PetscErrorCode DMAdaptInterpolator(DM dmc, DM dmf, Mat In, KSP smoother, Mat MF, Mat MC, Mat *InAdapt, void *user)
240: {
241:   Mat                globalA, AF;
242:   Vec                tmp;
243:   const PetscScalar *af, *ac;
244:   PetscScalar       *A, *b, *x, *workscalar;
245:   PetscReal         *w, *sing, *workreal, rcond = PETSC_SMALL;
246:   PetscBLASInt       M, N, one = 1, irank, lwrk, info;
247:   PetscInt           debug = 0, rStart, rEnd, r, maxcols = 0, k, Nc, ldac, ldaf;
248:   PetscBool          allocVc = PETSC_FALSE;

250:   PetscFunctionBegin;
251:   PetscCall(PetscLogEventBegin(DM_AdaptInterpolator, dmc, dmf, 0, 0));
252:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dm_interpolator_adapt_debug", &debug, NULL));
253:   PetscCall(MatGetSize(MF, NULL, &Nc));
254:   PetscCall(MatDuplicate(In, MAT_SHARE_NONZERO_PATTERN, InAdapt));
255:   PetscCall(MatGetOwnershipRange(In, &rStart, &rEnd));
256: #if 0
257:   PetscCall(MatGetMaxRowLen(In, &maxcols));
258: #else
259:   for (r = rStart; r < rEnd; ++r) {
260:     PetscInt ncols;

262:     PetscCall(MatGetRow(In, r, &ncols, NULL, NULL));
263:     maxcols = PetscMax(maxcols, ncols);
264:     PetscCall(MatRestoreRow(In, r, &ncols, NULL, NULL));
265:   }
266: #endif
267:   if (Nc < maxcols) PetscCall(PetscPrintf(PETSC_COMM_SELF, "The number of input vectors %" PetscInt_FMT " < %" PetscInt_FMT " the maximum number of column entries\n", Nc, maxcols));
268:   for (k = 0; k < Nc && debug; ++k) {
269:     char        name[PETSC_MAX_PATH_LEN];
270:     const char *prefix;
271:     Vec         vc, vf;

273:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)smoother, &prefix));

275:     if (MC) {
276:       PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sCoarse Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
277:       PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
278:       PetscCall(PetscObjectSetName((PetscObject)vc, name));
279:       PetscCall(VecViewFromOptions(vc, NULL, "-dm_adapt_interp_view_coarse"));
280:       PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
281:     }
282:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sFine Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
283:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
284:     PetscCall(PetscObjectSetName((PetscObject)vf, name));
285:     PetscCall(VecViewFromOptions(vf, NULL, "-dm_adapt_interp_view_fine"));
286:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
287:   }
288:   PetscCall(PetscBLASIntCast(3 * PetscMin(Nc, maxcols) + PetscMax(2 * PetscMin(Nc, maxcols), PetscMax(Nc, maxcols)), &lwrk));
289:   PetscCall(PetscMalloc7(Nc * maxcols, &A, PetscMax(Nc, maxcols), &b, Nc, &w, maxcols, &x, maxcols, &sing, lwrk, &workscalar, 5 * PetscMin(Nc, maxcols), &workreal));
290:   /* 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} */
291:   PetscCall(KSPGetOperators(smoother, &globalA, NULL));

293:   PetscCall(MatMatMult(globalA, MF, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AF));
294:   for (k = 0; k < Nc; ++k) {
295:     PetscScalar vnorm, vAnorm;
296:     Vec         vf;

298:     w[k] = 1.0;
299:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
300:     PetscCall(MatDenseGetColumnVecRead(AF, k, &tmp));
301:     PetscCall(VecDot(vf, vf, &vnorm));
302: #if 0
303:     PetscCall(DMGetGlobalVector(dmf, &tmp2));
304:     PetscCall(KSPSolve(smoother, tmp, tmp2));
305:     PetscCall(VecDot(vf, tmp2, &vAnorm));
306:     PetscCall(DMRestoreGlobalVector(dmf, &tmp2));
307: #else
308:     PetscCall(VecDot(vf, tmp, &vAnorm));
309: #endif
310:     w[k] = PetscRealPart(vnorm) / PetscRealPart(vAnorm);
311:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
312:     PetscCall(MatDenseRestoreColumnVecRead(AF, k, &tmp));
313:   }
314:   PetscCall(MatDestroy(&AF));
315:   if (!MC) {
316:     allocVc = PETSC_TRUE;
317:     PetscCall(MatTransposeMatMult(In, MF, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &MC));
318:   }
319:   /* Solve a LS system for each fine row
320:      MATT: Can we generalize to the case where Nc for the fine space
321:      is different for Nc for the coarse? */
322:   PetscCall(MatDenseGetArrayRead(MF, &af));
323:   PetscCall(MatDenseGetLDA(MF, &ldaf));
324:   PetscCall(MatDenseGetArrayRead(MC, &ac));
325:   PetscCall(MatDenseGetLDA(MC, &ldac));
326:   for (r = rStart; r < rEnd; ++r) {
327:     PetscInt           ncols, c;
328:     const PetscInt    *cols;
329:     const PetscScalar *vals;

331:     PetscCall(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:       b[k] = wk * af[r - rStart + k * ldaf];
338:       /* A_{kc} = \sqrt{w_k} f^{C,k}_c */
339:       /* TODO Must pull out VecScatter from In, scatter in vc[k] values up front, and access them indirectly just as in MatMult() */
340:       for (c = 0; c < ncols; ++c) {
341:         /* This is element (k, c) of A */
342:         A[c * Nc + k] = wk * ac[cols[c] - rStart + k * ldac];
343:       }
344:     }
345:     PetscCall(PetscBLASIntCast(Nc, &M));
346:     PetscCall(PetscBLASIntCast(ncols, &N));
347:     if (debug) {
348: #if defined(PETSC_USE_COMPLEX)
349:       PetscScalar *tmp;
350:       PetscInt     j;

352:       PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
353:       for (j = 0; j < Nc; ++j) tmp[j] = w[j];
354:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, tmp));
355:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
356:       for (j = 0; j < Nc; ++j) tmp[j] = b[j];
357:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, tmp));
358:       PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
359: #else
360:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, w));
361:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
362:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, b));
363: #endif
364:     }
365: #if defined(PETSC_USE_COMPLEX)
366:     /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
367:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, workreal, &info));
368: #else
369:     /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
370:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, &info));
371: #endif
372:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS");
373:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge");
374:     if (debug) {
375:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "rank %" PetscBLASInt_FMT " rcond %g\n", irank, (double)rcond));
376: #if defined(PETSC_USE_COMPLEX)
377:       {
378:         PetscScalar *tmp;
379:         PetscInt     j;

381:         PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
382:         for (j = 0; j < PetscMin(Nc, ncols); ++j) tmp[j] = sing[j];
383:         PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, tmp));
384:         PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
385:       }
386: #else
387:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, sing));
388: #endif
389:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS old P", ncols, 1, vals));
390:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS sol", ncols, 1, b));
391:     }
392:     PetscCall(MatSetValues(*InAdapt, 1, &r, ncols, cols, b, INSERT_VALUES));
393:     PetscCall(MatRestoreRow(In, r, &ncols, &cols, &vals));
394:   }
395:   PetscCall(MatDenseRestoreArrayRead(MF, &af));
396:   PetscCall(MatDenseRestoreArrayRead(MC, &ac));
397:   PetscCall(PetscFree7(A, b, w, x, sing, workscalar, workreal));
398:   if (allocVc) PetscCall(MatDestroy(&MC));
399:   PetscCall(MatAssemblyBegin(*InAdapt, MAT_FINAL_ASSEMBLY));
400:   PetscCall(MatAssemblyEnd(*InAdapt, MAT_FINAL_ASSEMBLY));
401:   PetscCall(PetscLogEventEnd(DM_AdaptInterpolator, dmc, dmf, 0, 0));
402:   PetscFunctionReturn(PETSC_SUCCESS);
403: }

405: PetscErrorCode DMCheckInterpolator(DM dmf, Mat In, Mat MC, Mat MF, PetscReal tol)
406: {
407:   Vec       tmp;
408:   PetscReal norminf, norm2, maxnorminf = 0.0, maxnorm2 = 0.0;
409:   PetscInt  k, Nc;

411:   PetscFunctionBegin;
412:   PetscCall(DMGetGlobalVector(dmf, &tmp));
413:   PetscCall(MatViewFromOptions(In, NULL, "-dm_interpolator_adapt_error"));
414:   PetscCall(MatGetSize(MF, NULL, &Nc));
415:   for (k = 0; k < Nc; ++k) {
416:     Vec vc, vf;

418:     PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
419:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
420:     PetscCall(MatMult(In, vc, tmp));
421:     PetscCall(VecAXPY(tmp, -1.0, vf));
422:     PetscCall(VecViewFromOptions(vc, NULL, "-dm_interpolator_adapt_error"));
423:     PetscCall(VecViewFromOptions(vf, NULL, "-dm_interpolator_adapt_error"));
424:     PetscCall(VecViewFromOptions(tmp, NULL, "-dm_interpolator_adapt_error"));
425:     PetscCall(VecNorm(tmp, NORM_INFINITY, &norminf));
426:     PetscCall(VecNorm(tmp, NORM_2, &norm2));
427:     maxnorminf = PetscMax(maxnorminf, norminf);
428:     maxnorm2   = PetscMax(maxnorm2, norm2);
429:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmf), "Coarse vec %" PetscInt_FMT " ||vf - P vc||_\\infty %g, ||vf - P vc||_2 %g\n", k, (double)norminf, (double)norm2));
430:     PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
431:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
432:   }
433:   PetscCall(DMRestoreGlobalVector(dmf, &tmp));
434:   PetscCheck(maxnorm2 <= tol, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_WRONG, "max_k ||vf_k - P vc_k||_2 %g > tol %g", (double)maxnorm2, (double)tol);
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }