Actual source code: convest.c
petsc-3.9.4 2018-09-11
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("-num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
95: PetscOptionsEnd();
96: return(0);
97: }
99: /*@
100: PetscConvEstView - Views a PetscConvEst object
102: Collective on PetscConvEst
104: Input Parameters:
105: + ce - The PetscConvEst object
106: - viewer - The PetscViewer object
108: Level: beginner
110: .keywords: PetscConvEst, convergence, view
111: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
112: @*/
113: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
114: {
118: PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
119: PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
120: return(0);
121: }
123: /*@
124: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
126: Not collective
128: Input Parameter:
129: . ce - The PetscConvEst object
131: Output Parameter:
132: . snes - The solver
134: Level: intermediate
136: .keywords: PetscConvEst, convergence
137: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
138: @*/
139: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes)
140: {
144: *snes = ce->snes;
145: return(0);
146: }
148: /*@
149: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
151: Not collective
153: Input Parameters:
154: + ce - The PetscConvEst object
155: - snes - The solver
157: Level: intermediate
159: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
161: .keywords: PetscConvEst, convergence
162: .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
163: @*/
164: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes)
165: {
171: ce->snes = snes;
172: SNESGetDM(ce->snes, &ce->idm);
173: return(0);
174: }
176: /*@
177: PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence
179: Collective on PetscConvEst
181: Input Parameters:
182: . ce - The PetscConvEst object
184: Level: beginner
186: .keywords: PetscConvEst, convergence, setup
187: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
188: @*/
189: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
190: {
191: PetscDS prob;
192: PetscInt f;
196: DMGetDS(ce->idm, &prob);
197: PetscDSGetNumFields(prob, &ce->Nf);
198: PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
199: PetscMalloc2(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol);
200: for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private;
201: for (f = 0; f < ce->Nf; ++f) {
202: PetscDSGetExactSolution(prob, f, &ce->exactSol[f]);
203: 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);
204: }
205: return(0);
206: }
208: static PetscErrorCode PetscConvEstLinearRegression_Private(PetscConvEst ce, PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept)
209: {
210: PetscScalar H[4];
211: PetscReal *X, *Y, beta[2];
212: PetscInt i, j, k;
216: *slope = *intercept = 0.0;
217: PetscMalloc2(n*2, &X, n*2, &Y);
218: for (k = 0; k < n; ++k) {
219: /* X[n,2] = [1, x] */
220: X[k*2+0] = 1.0;
221: X[k*2+1] = x[k];
222: }
223: /* H = X^T X */
224: for (i = 0; i < 2; ++i) {
225: for (j = 0; j < 2; ++j) {
226: H[i*2+j] = 0.0;
227: for (k = 0; k < n; ++k) {
228: H[i*2+j] += X[k*2+i] * X[k*2+j];
229: }
230: }
231: }
232: /* H = (X^T X)^{-1} */
233: {
234: PetscBLASInt two = 2, ipiv[2], info;
235: PetscScalar work[2];
237: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
238: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
239: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
240: PetscFPTrapPop();
241: }
242: /* Y = H X^T */
243: for (i = 0; i < 2; ++i) {
244: for (k = 0; k < n; ++k) {
245: Y[i*n+k] = 0.0;
246: for (j = 0; j < 2; ++j) {
247: Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j];
248: }
249: }
250: }
251: /* beta = Y error = [y-intercept, slope] */
252: for (i = 0; i < 2; ++i) {
253: beta[i] = 0.0;
254: for (k = 0; k < n; ++k) {
255: beta[i] += Y[i*n+k] * y[k];
256: }
257: }
258: PetscFree2(X, Y);
259: *intercept = beta[0];
260: *slope = beta[1];
261: return(0);
262: }
264: /*@
265: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
267: Not collective
269: Input Parameter:
270: . ce - The PetscConvEst object
272: Output Parameter:
273: . alpha - The convergence rate
275: Note: The convergence rate alpha is defined by
276: $ || u_h - u_exact || < C h^alpha
277: where u_h is the discrete solution, and h is a measure of the discretization size.
279: We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS,
280: and then fit the result to our model above using linear regression.
282: Options database keys:
283: . -snes_convergence_estimate : Execute convergence estimation and print out the rate
285: Level: intermediate
287: .keywords: PetscConvEst, convergence
288: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
289: @*/
290: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal *alpha)
291: {
292: DM *dm;
293: PetscDS prob;
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;
304: PetscObjectGetComm((PetscObject) ce, &comm);
305: DMGetDimension(ce->idm, &dim);
306: DMGetApplicationContext(ce->idm, &ctx);
307: DMGetDS(ce->idm, &prob);
308: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
309: PetscMalloc2((Nr+1), &dm, (Nr+1), &dof);
310: dm[0] = ce->idm;
311: *alpha = 0.0;
312: /* Loop over meshes */
313: for (r = 0; r <= Nr; ++r) {
314: if (r > 0) {
315: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
316: DMSetCoarseDM(dm[r], dm[r-1]);
317: DMSetDS(dm[r], prob);
318: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
319: PetscObjectSetName((PetscObject) dm[r], dmname);
320: }
321: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
322: DMPlexGetHeightStratum(dm[r], 0, NULL, &dof[r]);
323: /* Create solution */
324: DMCreateGlobalVector(dm[r], &u);
325: PetscDSGetDiscretization(prob, 0, &disc);
326: PetscObjectGetName(disc, &uname);
327: PetscObjectSetName((PetscObject) u, uname);
328: /* Setup solver */
329: SNESReset(ce->snes);
330: SNESSetDM(ce->snes, dm[r]);
331: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
332: SNESSetFromOptions(ce->snes);
333: /* Create initial guess */
334: DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
335: SNESSolve(ce->snes, NULL, u);
336: DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);
337: /* Monitor */
338: if (ce->monitor) {
339: PetscReal *errors = &ce->errors[r*ce->Nf];
340: PetscInt f;
342: PetscPrintf(comm, "L_2 Error: [");
343: for (f = 0; f < ce->Nf; ++f) {
344: if (f > 0) {PetscPrintf(comm, ", ");}
345: if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
346: else {PetscPrintf(comm, "%g", (double)errors[f]);}
347: }
348: PetscPrintf(comm, "]\n");
349: }
350: /* Cleanup */
351: VecDestroy(&u);
352: }
353: for (r = 1; r <= Nr; ++r) {
354: DMDestroy(&dm[r]);
355: }
356: /* Fit convergence rate */
357: PetscMalloc2(Nr+1, &x, Nr+1, &y);
358: for (r = 0; r <= Nr; ++r) {
359: x[r] = PetscLog10Real(dof[r]);
360: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+0]);
361: }
362: PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
363: PetscFree2(x, y);
364: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
365: *alpha = -slope * dim;
366: PetscFree2(dm, dof);
367: /* Restore solver */
368: SNESReset(ce->snes);
369: SNESSetDM(ce->snes, ce->idm);
370: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
371: SNESSetFromOptions(ce->snes);
372: return(0);
373: }
375: /*@
376: PetscConvEstRateView - Displays the convergence rate to a viewer
378: Collective on SNES
380: Parameter:
381: + snes - iterative context obtained from SNESCreate()
382: . alpha - the convergence rate
383: - viewer - the viewer to display the reason
385: Options Database Keys:
386: . -snes_convergence_estimate - print the convergence rate
388: Level: developer
390: .seealso: PetscConvEstGetRate()
391: @*/
392: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, PetscReal alpha, PetscViewer viewer)
393: {
394: PetscBool isAscii;
398: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
399: if (isAscii) {
400: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
401: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: %g\n", (double) alpha);
402: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
403: }
404: return(0);
405: }