Actual source code: convest.c
1: #include <petscconvest.h>
2: #include <petscdmplex.h>
3: #include <petscds.h>
5: #include <petsc/private/petscconvestimpl.h>
7: static PetscErrorCode zero_private(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
8: {
9: PetscInt c;
10: for (c = 0; c < Nc; ++c) u[c] = 0.0;
11: return 0;
12: }
14: /*@
15: PetscConvEstDestroy - Destroys a PetscConvEst object
17: Collective on PetscConvEst
19: Input Parameter:
20: . ce - The PetscConvEst object
22: Level: beginner
24: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
25: @*/
26: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
27: {
28: if (!*ce) return 0;
30: if (--((PetscObject)(*ce))->refct > 0) {
31: *ce = NULL;
32: return 0;
33: }
34: PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs);
35: PetscFree2((*ce)->dofs, (*ce)->errors);
36: PetscHeaderDestroy(ce);
37: return 0;
38: }
40: /*@
41: PetscConvEstSetFromOptions - Sets a PetscConvEst object from options
43: Collective on PetscConvEst
45: Input Parameters:
46: . ce - The PetscConvEst object
48: Level: beginner
50: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
51: @*/
52: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
53: {
56: PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");
57: PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
58: PetscOptionsReal("-convest_refine_factor", "The increase in resolution in each dimension", "PetscConvEst", ce->r, &ce->r, NULL);
59: PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);
60: PetscOptionsBool("-convest_no_refine", "Debugging flag to run on the same mesh each time", "PetscConvEst", ce->noRefine, &ce->noRefine, NULL);
61: PetscOptionsEnd();
62: return 0;
63: }
65: /*@
66: PetscConvEstView - Views a PetscConvEst object
68: Collective on PetscConvEst
70: Input Parameters:
71: + ce - The PetscConvEst object
72: - viewer - The PetscViewer object
74: Level: beginner
76: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
77: @*/
78: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
79: {
80: PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
81: PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
82: return 0;
83: }
85: /*@
86: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
88: Not collective
90: Input Parameter:
91: . ce - The PetscConvEst object
93: Output Parameter:
94: . solver - The solver
96: Level: intermediate
98: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
99: @*/
100: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
101: {
104: *solver = ce->solver;
105: return 0;
106: }
108: /*@
109: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
111: Not collective
113: Input Parameters:
114: + ce - The PetscConvEst object
115: - solver - The solver
117: Level: intermediate
119: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
121: .seealso: PetscConvEstGetSNES(), PetscConvEstCreate(), PetscConvEstGetConvRate()
122: @*/
123: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
124: {
127: ce->solver = solver;
128: (*ce->ops->setsolver)(ce, solver);
129: return 0;
130: }
132: /*@
133: PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence
135: Collective on PetscConvEst
137: Input Parameters:
138: . ce - The PetscConvEst object
140: Level: beginner
142: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
143: @*/
144: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
145: {
146: PetscInt Nf, f, Nds, s;
148: DMGetNumFields(ce->idm, &Nf);
149: ce->Nf = PetscMax(Nf, 1);
150: PetscMalloc2((ce->Nr+1)*ce->Nf, &ce->dofs, (ce->Nr+1)*ce->Nf, &ce->errors);
151: PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);
152: for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private;
153: DMGetNumDS(ce->idm, &Nds);
154: for (s = 0; s < Nds; ++s) {
155: PetscDS ds;
156: DMLabel label;
157: IS fieldIS;
158: const PetscInt *fields;
159: PetscInt dsNf;
161: DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds);
162: PetscDSGetNumFields(ds, &dsNf);
163: if (fieldIS) ISGetIndices(fieldIS, &fields);
164: for (f = 0; f < dsNf; ++f) {
165: const PetscInt field = fields[f];
166: PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]);
167: }
168: if (fieldIS) ISRestoreIndices(fieldIS, &fields);
169: }
170: for (f = 0; f < Nf; ++f) {
172: }
173: return 0;
174: }
176: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
177: {
181: (*ce->ops->initguess)(ce, r, dm, u);
182: return 0;
183: }
185: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
186: {
191: (*ce->ops->computeerror)(ce, r, dm, u, errors);
192: return 0;
193: }
195: /*@
196: PetscConvEstMonitorDefault - Monitors the convergence estimation loop
198: Collective on PetscConvEst
200: Input Parameters:
201: + ce - The PetscConvEst object
202: - r - The refinement level
204: Options database keys:
205: . -convest_monitor - Activate the monitor
207: Level: intermediate
209: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve()
210: @*/
211: PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
212: {
213: MPI_Comm comm;
214: PetscInt f;
216: if (ce->monitor) {
217: PetscInt *dofs = &ce->dofs[r*ce->Nf];
218: PetscReal *errors = &ce->errors[r*ce->Nf];
220: PetscObjectGetComm((PetscObject) ce, &comm);
221: PetscPrintf(comm, "N: ");
222: if (ce->Nf > 1) PetscPrintf(comm, "[");
223: for (f = 0; f < ce->Nf; ++f) {
224: if (f > 0) PetscPrintf(comm, ", ");
225: PetscPrintf(comm, "%7D", dofs[f]);
226: }
227: if (ce->Nf > 1) PetscPrintf(comm, "]");
228: PetscPrintf(comm, " ");
229: PetscPrintf(comm, "L_2 Error: ");
230: if (ce->Nf > 1) PetscPrintf(comm, "[");
231: for (f = 0; f < ce->Nf; ++f) {
232: if (f > 0) PetscPrintf(comm, ", ");
233: if (errors[f] < 1.0e-11) PetscPrintf(comm, "< 1e-11");
234: else PetscPrintf(comm, "%g", (double) errors[f]);
235: }
236: if (ce->Nf > 1) PetscPrintf(comm, "]");
237: PetscPrintf(comm, "\n");
238: }
239: return 0;
240: }
242: static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
243: {
244: PetscClassId id;
246: PetscObjectGetClassId(ce->solver, &id);
248: SNESGetDM((SNES) ce->solver, &ce->idm);
249: return 0;
250: }
252: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
253: {
254: DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u);
255: return 0;
256: }
258: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
259: {
260: DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors);
261: return 0;
262: }
264: static PetscErrorCode PetscConvEstSetJacobianNullspace_Private(PetscConvEst ce, SNES snes)
265: {
266: DM dm;
267: PetscInt f;
269: SNESGetDM(snes, &dm);
270: for (f = 0; f < ce->Nf; ++f) {
271: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
273: DMGetNullSpaceConstructor(dm, f, &nspconstr);
274: if (nspconstr) {
275: MatNullSpace nullsp;
276: Mat J;
278: (*nspconstr)(dm, f, f,&nullsp);
279: SNESSetUp(snes);
280: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
281: MatSetNullSpace(J, nullsp);
282: MatNullSpaceDestroy(&nullsp);
283: break;
284: }
285: }
286: return 0;
287: }
289: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
290: {
291: SNES snes = (SNES) ce->solver;
292: DM *dm;
293: PetscObject disc;
294: PetscReal *x, *y, slope, intercept;
295: PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
296: void *ctx;
299: DMGetDimension(ce->idm, &dim);
300: DMGetApplicationContext(ce->idm, &ctx);
301: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
302: DMGetRefineLevel(ce->idm, &oldlevel);
303: PetscMalloc1((Nr+1), &dm);
304: /* Loop over meshes */
305: dm[0] = ce->idm;
306: for (r = 0; r <= Nr; ++r) {
307: Vec u;
308: #if defined(PETSC_USE_LOG)
309: PetscLogStage stage;
310: #endif
311: char stageName[PETSC_MAX_PATH_LEN];
312: const char *dmname, *uname;
314: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
315: #if defined(PETSC_USE_LOG)
316: PetscLogStageGetId(stageName, &stage);
317: if (stage < 0) PetscLogStageRegister(stageName, &stage);
318: #endif
319: PetscLogStagePush(stage);
320: if (r > 0) {
321: if (!ce->noRefine) {
322: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
323: DMSetCoarseDM(dm[r], dm[r-1]);
324: } else {
325: DM cdm, rcdm;
327: DMClone(dm[r-1], &dm[r]);
328: DMCopyDisc(dm[r-1], dm[r]);
329: DMGetCoordinateDM(dm[r-1], &cdm);
330: DMGetCoordinateDM(dm[r], &rcdm);
331: DMCopyDisc(cdm, rcdm);
332: }
333: DMCopyTransform(ce->idm, dm[r]);
334: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
335: PetscObjectSetName((PetscObject) dm[r], dmname);
336: for (f = 0; f < ce->Nf; ++f) {
337: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
339: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
340: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
341: }
342: }
343: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
344: /* Create solution */
345: DMCreateGlobalVector(dm[r], &u);
346: DMGetField(dm[r], 0, NULL, &disc);
347: PetscObjectGetName(disc, &uname);
348: PetscObjectSetName((PetscObject) u, uname);
349: /* Setup solver */
350: SNESReset(snes);
351: SNESSetDM(snes, dm[r]);
352: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
353: SNESSetFromOptions(snes);
354: /* Set nullspace for Jacobian */
355: PetscConvEstSetJacobianNullspace_Private(ce, snes);
356: /* Create initial guess */
357: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
358: SNESSolve(snes, NULL, u);
359: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
360: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*ce->Nf]);
361: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
362: for (f = 0; f < ce->Nf; ++f) {
363: PetscSection s, fs;
364: PetscInt lsize;
366: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
367: DMGetLocalSection(dm[r], &s);
368: PetscSectionGetField(s, f, &fs);
369: PetscSectionGetConstrainedStorageSize(fs, &lsize);
370: MPI_Allreduce(&lsize, &ce->dofs[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) snes));
371: PetscLogEventSetDof(ce->event, f, ce->dofs[r*ce->Nf+f]);
372: PetscLogEventSetError(ce->event, f, ce->errors[r*ce->Nf+f]);
373: }
374: /* Monitor */
375: PetscConvEstMonitorDefault(ce, r);
376: if (!r) {
377: /* PCReset() does not wipe out the level structure */
378: KSP ksp;
379: PC pc;
381: SNESGetKSP(snes, &ksp);
382: KSPGetPC(ksp, &pc);
383: PCMGGetLevels(pc, &oldnlev);
384: }
385: /* Cleanup */
386: VecDestroy(&u);
387: PetscLogStagePop();
388: }
389: for (r = 1; r <= Nr; ++r) {
390: DMDestroy(&dm[r]);
391: }
392: /* Fit convergence rate */
393: PetscMalloc2(Nr+1, &x, Nr+1, &y);
394: for (f = 0; f < ce->Nf; ++f) {
395: for (r = 0; r <= Nr; ++r) {
396: x[r] = PetscLog10Real(ce->dofs[r*ce->Nf+f]);
397: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
398: }
399: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
400: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
401: alpha[f] = -slope * dim;
402: }
403: PetscFree2(x, y);
404: PetscFree(dm);
405: /* Restore solver */
406: SNESReset(snes);
407: {
408: /* PCReset() does not wipe out the level structure */
409: KSP ksp;
410: PC pc;
412: SNESGetKSP(snes, &ksp);
413: KSPGetPC(ksp, &pc);
414: PCMGSetLevels(pc, oldnlev, NULL);
415: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
416: }
417: SNESSetDM(snes, ce->idm);
418: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
419: SNESSetFromOptions(snes);
420: PetscConvEstSetJacobianNullspace_Private(ce, snes);
421: return 0;
422: }
424: /*@
425: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
427: Not collective
429: Input Parameter:
430: . ce - The PetscConvEst object
432: Output Parameter:
433: . alpha - The convergence rate for each field
435: Note: The convergence rate alpha is defined by
436: $ || u_\Delta - u_exact || < C \Delta^alpha
437: where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the
438: spatial resolution and \Delta t for the temporal resolution.
440: We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
441: based upon the exact solution in the DS, and then fit the result to our model above using linear regression.
443: Options database keys:
444: + -snes_convergence_estimate - Execute convergence estimation inside SNESSolve() and print out the rate
445: - -ts_convergence_estimate - Execute convergence estimation inside TSSolve() and print out the rate
447: Level: intermediate
449: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve()
450: @*/
451: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
452: {
453: PetscInt f;
455: if (ce->event < 0) PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);
456: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
457: (*ce->ops->getconvrate)(ce, alpha);
458: return 0;
459: }
461: /*@
462: PetscConvEstRateView - Displays the convergence rate to a viewer
464: Collective on SNES
466: Parameter:
467: + snes - iterative context obtained from SNESCreate()
468: . alpha - the convergence rate for each field
469: - viewer - the viewer to display the reason
471: Options Database Keys:
472: . -snes_convergence_estimate - print the convergence rate
474: Level: developer
476: .seealso: PetscConvEstGetRate()
477: @*/
478: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
479: {
480: PetscBool isAscii;
482: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
483: if (isAscii) {
484: PetscInt Nf = ce->Nf, f;
486: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
487: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
488: if (Nf > 1) PetscViewerASCIIPrintf(viewer, "[");
489: for (f = 0; f < Nf; ++f) {
490: if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
491: PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);
492: }
493: if (Nf > 1) PetscViewerASCIIPrintf(viewer, "]");
494: PetscViewerASCIIPrintf(viewer, "\n");
495: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
496: }
497: return 0;
498: }
500: /*@
501: PetscConvEstCreate - Create a PetscConvEst object
503: Collective
505: Input Parameter:
506: . comm - The communicator for the PetscConvEst object
508: Output Parameter:
509: . ce - The PetscConvEst object
511: Level: beginner
513: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
514: @*/
515: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
516: {
518: PetscSysInitializePackage();
519: PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
520: (*ce)->monitor = PETSC_FALSE;
521: (*ce)->r = 2.0;
522: (*ce)->Nr = 4;
523: (*ce)->event = -1;
524: (*ce)->ops->setsolver = PetscConvEstSetSNES_Private;
525: (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private;
526: (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
527: (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private;
528: return 0;
529: }