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