Actual source code: convest.c
petsc-3.10.5 2019-03-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: PetscDS prob;
295: PetscObject disc;
296: MPI_Comm comm;
297: const char *uname, *dmname;
298: void *ctx;
299: Vec u;
300: PetscReal t = 0.0, *x, *y, slope, intercept;
301: PetscInt *dof, dim, Nr = ce->Nr, r, f;
302: PetscLogEvent event;
306: PetscObjectGetComm((PetscObject) ce, &comm);
307: DMGetDimension(ce->idm, &dim);
308: DMGetApplicationContext(ce->idm, &ctx);
309: DMGetDS(ce->idm, &prob);
310: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
311: PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);
312: dm[0] = ce->idm;
313: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
314: /* Loop over meshes */
315: PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);
316: for (r = 0; r <= Nr; ++r) {
317: PetscLogStage stage;
318: char stageName[PETSC_MAX_PATH_LEN];
320: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
321: PetscLogStageRegister(stageName, &stage);
322: PetscLogStagePush(stage);
323: if (r > 0) {
324: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
325: DMSetCoarseDM(dm[r], dm[r-1]);
326: DMSetDS(dm[r], prob);
327: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
328: PetscObjectSetName((PetscObject) dm[r], dmname);
329: for (f = 0; f <= ce->Nf; ++f) {
330: PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
331: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
332: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
333: }
334: }
335: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
336: /* Create solution */
337: DMCreateGlobalVector(dm[r], &u);
338: PetscDSGetDiscretization(prob, 0, &disc);
339: PetscObjectGetName(disc, &uname);
340: PetscObjectSetName((PetscObject) u, uname);
341: /* Setup solver */
342: SNESReset(ce->snes);
343: SNESSetDM(ce->snes, dm[r]);
344: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
345: SNESSetFromOptions(ce->snes);
346: /* Create initial guess */
347: DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
348: SNESSolve(ce->snes, NULL, u);
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: DMPlexGetHeightStratum(dm[r], 0, NULL, &dof[r*ce->Nf+f]);
353: }
354: /* Monitor */
355: if (ce->monitor) {
356: PetscReal *errors = &ce->errors[r*ce->Nf];
358: PetscPrintf(comm, "L_2 Error: ");
359: if (ce->Nf > 1) {PetscPrintf(comm, "[");}
360: for (f = 0; f < ce->Nf; ++f) {
361: if (f > 0) {PetscPrintf(comm, ", ");}
362: if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
363: else {PetscPrintf(comm, "%g", (double)errors[f]);}
364: }
365: if (ce->Nf > 1) {PetscPrintf(comm, "]");}
366: PetscPrintf(comm, "\n");
367: }
368: /* Cleanup */
369: VecDestroy(&u);
370: PetscLogStagePop();
371: }
372: for (r = 1; r <= Nr; ++r) {
373: DMDestroy(&dm[r]);
374: }
375: /* Fit convergence rate */
376: PetscMalloc2(Nr+1, &x, Nr+1, &y);
377: for (f = 0; f < ce->Nf; ++f) {
378: for (r = 0; r <= Nr; ++r) {
379: x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
380: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
381: }
382: PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
383: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
384: alpha[f] = -slope * dim;
385: }
386: PetscFree2(x, y);
387: PetscFree2(dm, dof);
388: /* Restore solver */
389: SNESReset(ce->snes);
390: SNESSetDM(ce->snes, ce->idm);
391: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
392: SNESSetFromOptions(ce->snes);
393: return(0);
394: }
396: /*@
397: PetscConvEstRateView - Displays the convergence rate to a viewer
399: Collective on SNES
401: Parameter:
402: + snes - iterative context obtained from SNESCreate()
403: . alpha - the convergence rate for each field
404: - viewer - the viewer to display the reason
406: Options Database Keys:
407: . -snes_convergence_estimate - print the convergence rate
409: Level: developer
411: .seealso: PetscConvEstGetRate()
412: @*/
413: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, PetscReal alpha[], PetscViewer viewer)
414: {
415: PetscBool isAscii;
416: PetscInt f;
420: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
421: if (isAscii) {
422: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
423: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
424: if (ce->Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
425: for (f = 0; f < ce->Nf; ++f) {
426: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
427: PetscViewerASCIIPrintf(viewer, "%g", (double) alpha[f]);
428: }
429: if (ce->Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
430: PetscViewerASCIIPrintf(viewer, "\n");
431: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
432: }
433: return(0);
434: }