Actual source code: convest.c
petsc-3.11.4 2019-09-28
1: #include <petscconvest.h>
2: #include <petscdmplex.h>
3: #include <petscds.h>
4: #include <petscblaslapack.h>
6: #include <petsc/private/petscconvestimpl.h>
8: static PetscErrorCode zero_private(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
9: {
10: PetscInt c;
11: for (c = 0; c < Nc; ++c) u[c] = 0.0;
12: return 0;
13: }
16: /*@
17: PetscConvEstCreate - Create a PetscConvEst object
19: Collective on MPI_Comm
21: Input Parameter:
22: . comm - The communicator for the PetscConvEst object
24: Output Parameter:
25: . ce - The PetscConvEst object
27: Level: beginner
29: .keywords: PetscConvEst, convergence, create
30: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
31: @*/
32: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
33: {
38: PetscSysInitializePackage();
39: PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
40: (*ce)->monitor = PETSC_FALSE;
41: (*ce)->Nr = 4;
42: return(0);
43: }
45: /*@
46: PetscConvEstDestroy - Destroys a PetscConvEst object
48: Collective on PetscConvEst
50: Input Parameter:
51: . ce - The PetscConvEst object
53: Level: beginner
55: .keywords: PetscConvEst, convergence, destroy
56: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
57: @*/
58: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
59: {
63: if (!*ce) return(0);
65: if (--((PetscObject)(*ce))->refct > 0) {
66: *ce = NULL;
67: return(0);
68: }
69: PetscFree2((*ce)->initGuess, (*ce)->exactSol);
70: PetscFree((*ce)->errors);
71: PetscHeaderDestroy(ce);
72: return(0);
73: }
75: /*@
76: PetscConvEstSetFromOptions - Sets a PetscConvEst object from options
78: Collective on PetscConvEst
80: Input Parameters:
81: . ce - The PetscConvEst object
83: Level: beginner
85: .keywords: PetscConvEst, convergence, options
86: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
87: @*/
88: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
89: {
93: PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");
94: PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
95: PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);
96: PetscOptionsEnd();
97: return(0);
98: }
100: /*@
101: PetscConvEstView - Views a PetscConvEst object
103: Collective on PetscConvEst
105: Input Parameters:
106: + ce - The PetscConvEst object
107: - viewer - The PetscViewer object
109: Level: beginner
111: .keywords: PetscConvEst, convergence, view
112: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
113: @*/
114: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
115: {
119: PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
120: PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
121: return(0);
122: }
124: /*@
125: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
127: Not collective
129: Input Parameter:
130: . ce - The PetscConvEst object
132: Output Parameter:
133: . snes - The solver
135: Level: intermediate
137: .keywords: PetscConvEst, convergence
138: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
139: @*/
140: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes)
141: {
145: *snes = ce->snes;
146: return(0);
147: }
149: /*@
150: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
152: Not collective
154: Input Parameters:
155: + ce - The PetscConvEst object
156: - snes - The solver
158: Level: intermediate
160: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
162: .keywords: PetscConvEst, convergence
163: .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
164: @*/
165: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes)
166: {
172: ce->snes = snes;
173: SNESGetDM(ce->snes, &ce->idm);
174: return(0);
175: }
177: /*@
178: PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence
180: Collective on PetscConvEst
182: Input Parameters:
183: . ce - The PetscConvEst object
185: Level: beginner
187: .keywords: PetscConvEst, convergence, setup
188: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
189: @*/
190: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
191: {
192: PetscDS prob;
193: PetscInt f;
197: DMGetDS(ce->idm, &prob);
198: PetscDSGetNumFields(prob, &ce->Nf);
199: PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
200: PetscMalloc2(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol);
201: for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private;
202: for (f = 0; f < ce->Nf; ++f) {
203: PetscDSGetExactSolution(prob, f, &ce->exactSol[f]);
204: if (!ce->exactSol[f]) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "DS must contain exact solution functions in order to estimate convergence, missing for field %D", f);
205: }
206: return(0);
207: }
209: static PetscErrorCode PetscConvEstLinearRegression_Private(PetscConvEst ce, PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept)
210: {
211: PetscScalar H[4];
212: PetscReal *X, *Y, beta[2];
213: PetscInt i, j, k;
217: *slope = *intercept = 0.0;
218: PetscMalloc2(n*2, &X, n*2, &Y);
219: for (k = 0; k < n; ++k) {
220: /* X[n,2] = [1, x] */
221: X[k*2+0] = 1.0;
222: X[k*2+1] = x[k];
223: }
224: /* H = X^T X */
225: for (i = 0; i < 2; ++i) {
226: for (j = 0; j < 2; ++j) {
227: H[i*2+j] = 0.0;
228: for (k = 0; k < n; ++k) {
229: H[i*2+j] += X[k*2+i] * X[k*2+j];
230: }
231: }
232: }
233: /* H = (X^T X)^{-1} */
234: {
235: PetscBLASInt two = 2, ipiv[2], info;
236: PetscScalar work[2];
238: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
239: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
240: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
241: PetscFPTrapPop();
242: }
243: /* Y = H X^T */
244: for (i = 0; i < 2; ++i) {
245: for (k = 0; k < n; ++k) {
246: Y[i*n+k] = 0.0;
247: for (j = 0; j < 2; ++j) {
248: Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j];
249: }
250: }
251: }
252: /* beta = Y error = [y-intercept, slope] */
253: for (i = 0; i < 2; ++i) {
254: beta[i] = 0.0;
255: for (k = 0; k < n; ++k) {
256: beta[i] += Y[i*n+k] * y[k];
257: }
258: }
259: PetscFree2(X, Y);
260: *intercept = beta[0];
261: *slope = beta[1];
262: return(0);
263: }
265: /*@
266: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
268: Not collective
270: Input Parameter:
271: . ce - The PetscConvEst object
273: Output Parameter:
274: . alpha - The convergence rate for each field
276: Note: The convergence rate alpha is defined by
277: $ || u_h - u_exact || < C h^alpha
278: where u_h is the discrete solution, and h is a measure of the discretization size.
280: We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS,
281: and then fit the result to our model above using linear regression.
283: Options database keys:
284: . -snes_convergence_estimate : Execute convergence estimation and print out the rate
286: Level: intermediate
288: .keywords: PetscConvEst, convergence
289: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
290: @*/
291: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
292: {
293: DM *dm;
294: PetscObject disc;
295: MPI_Comm comm;
296: const char *uname, *dmname;
297: void *ctx;
298: Vec u;
299: PetscReal t = 0.0, *x, *y, slope, intercept;
300: PetscInt *dof, dim, Nr = ce->Nr, r, f, oldlevel, oldnlev;
301: PetscLogEvent event;
305: PetscObjectGetComm((PetscObject) ce, &comm);
306: DMGetDimension(ce->idm, &dim);
307: DMGetApplicationContext(ce->idm, &ctx);
308: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
309: DMGetRefineLevel(ce->idm, &oldlevel);
310: PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);
311: dm[0] = ce->idm;
312: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
313: /* Loop over meshes */
314: PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);
315: for (r = 0; r <= Nr; ++r) {
316: PetscLogStage stage;
317: char stageName[PETSC_MAX_PATH_LEN];
319: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
320: PetscLogStageRegister(stageName, &stage);
321: PetscLogStagePush(stage);
322: if (r > 0) {
323: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
324: DMSetCoarseDM(dm[r], dm[r-1]);
325: DMCopyDisc(ce->idm, dm[r]);
326: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
327: PetscObjectSetName((PetscObject) dm[r], dmname);
328: for (f = 0; f <= ce->Nf; ++f) {
329: PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
330: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
331: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
332: }
333: }
334: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
335: /* Create solution */
336: DMCreateGlobalVector(dm[r], &u);
337: DMGetField(dm[r], 0, NULL, &disc);
338: PetscObjectGetName(disc, &uname);
339: PetscObjectSetName((PetscObject) u, uname);
340: /* Setup solver */
341: SNESReset(ce->snes);
342: SNESSetDM(ce->snes, dm[r]);
343: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
344: SNESSetFromOptions(ce->snes);
345: /* Create initial guess */
346: DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
347: SNESSolve(ce->snes, NULL, u);
348: PetscLogEventBegin(event, ce, 0, 0, 0);
349: DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);
350: PetscLogEventEnd(event, ce, 0, 0, 0);
351: for (f = 0; f < ce->Nf; ++f) {
352: PetscSection s, fs;
353: PetscInt lsize;
355: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
356: DMGetSection(dm[r], &s);
357: PetscSectionGetField(s, f, &fs);
358: PetscSectionGetConstrainedStorageSize(fs, &lsize);
359: MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));
360: PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);
361: PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);
362: }
363: /* Monitor */
364: if (ce->monitor) {
365: PetscReal *errors = &ce->errors[r*ce->Nf];
367: PetscPrintf(comm, "L_2 Error: ");
368: if (ce->Nf > 1) {PetscPrintf(comm, "[");}
369: for (f = 0; f < ce->Nf; ++f) {
370: if (f > 0) {PetscPrintf(comm, ", ");}
371: if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
372: else {PetscPrintf(comm, "%g", (double)errors[f]);}
373: }
374: if (ce->Nf > 1) {PetscPrintf(comm, "]");}
375: PetscPrintf(comm, "\n");
376: }
377: if (!r) {
378: /* PCReset() does not wipe out the level structure */
379: KSP ksp;
380: PC pc;
382: SNESGetKSP(ce->snes, &ksp);
383: KSPGetPC(ksp, &pc);
384: PCMGGetLevels(pc, &oldnlev);
385: }
386: /* Cleanup */
387: VecDestroy(&u);
388: PetscLogStagePop();
389: }
390: for (r = 1; r <= Nr; ++r) {
391: DMDestroy(&dm[r]);
392: }
393: /* Fit convergence rate */
394: PetscMalloc2(Nr+1, &x, Nr+1, &y);
395: for (f = 0; f < ce->Nf; ++f) {
396: for (r = 0; r <= Nr; ++r) {
397: x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
398: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
399: }
400: PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
401: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
402: alpha[f] = -slope * dim;
403: }
404: PetscFree2(x, y);
405: PetscFree2(dm, dof);
406: /* Restore solver */
407: SNESReset(ce->snes);
408: {
409: /* PCReset() does not wipe out the level structure */
410: KSP ksp;
411: PC pc;
413: SNESGetKSP(ce->snes, &ksp);
414: KSPGetPC(ksp, &pc);
415: PCMGSetLevels(pc, oldnlev, NULL);
416: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
417: }
418: SNESSetDM(ce->snes, ce->idm);
419: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
420: SNESSetFromOptions(ce->snes);
421: return(0);
422: }
424: /*@
425: PetscConvEstRateView - Displays the convergence rate to a viewer
427: Collective on SNES
429: Parameter:
430: + snes - iterative context obtained from SNESCreate()
431: . alpha - the convergence rate for each field
432: - viewer - the viewer to display the reason
434: Options Database Keys:
435: . -snes_convergence_estimate - print the convergence rate
437: Level: developer
439: .seealso: PetscConvEstGetRate()
440: @*/
441: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
442: {
443: PetscBool isAscii;
447: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
448: if (isAscii) {
449: DM dm;
450: PetscInt Nf, f;
452: SNESGetDM(ce->snes, &dm);
453: DMGetNumFields(dm, &Nf);
454: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
455: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
456: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
457: for (f = 0; f < Nf; ++f) {
458: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
459: PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);
460: }
461: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
462: PetscViewerASCIIPrintf(viewer, "\n");
463: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
464: }
465: return(0);
466: }