Actual source code: convest.c
1: #include "petscsys.h"
2: #include <petscconvest.h>
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: #include <petsc/private/petscconvestimpl.h>
8: static PetscErrorCode zero_private(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
9: {
10: for (PetscInt c = 0; c < Nc; ++c) u[c] = 0.0;
11: return PETSC_SUCCESS;
12: }
14: /*@
15: PetscConvEstDestroy - Destroys a PETSc convergence estimator `PetscConvEst` object
17: Collective
19: Input Parameter:
20: . ce - The `PetscConvEst` object
22: Level: beginner
24: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
25: @*/
26: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
27: {
28: PetscFunctionBegin;
29: if (!*ce) PetscFunctionReturn(PETSC_SUCCESS);
31: if (--((PetscObject)*ce)->refct > 0) {
32: *ce = NULL;
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
35: PetscCall(PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs));
36: PetscCall(PetscFree2((*ce)->dofs, (*ce)->errors));
37: PetscCall(PetscHeaderDestroy(ce));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /*@
42: PetscConvEstSetFromOptions - Sets a convergence estimator `PetscConvEst` object based on values in the options database
44: Collective
46: Input Parameter:
47: . ce - The `PetscConvEst` object
49: Level: beginner
51: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
52: @*/
53: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
54: {
55: PetscFunctionBegin;
56: PetscOptionsBegin(PetscObjectComm((PetscObject)ce), "", "Convergence Estimator Options", "PetscConvEst");
57: PetscCall(PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL));
58: PetscCall(PetscOptionsReal("-convest_refine_factor", "The increase in resolution in each dimension", "PetscConvEst", ce->r, &ce->r, NULL));
59: PetscCall(PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL));
60: PetscCall(PetscOptionsBool("-convest_no_refine", "Debugging flag to run on the same mesh each time", "PetscConvEst", ce->noRefine, &ce->noRefine, NULL));
61: PetscOptionsEnd();
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@
66: PetscConvEstView - Views a `PetscConvEst` object
68: Collective
70: Input Parameters:
71: + ce - The `PetscConvEst` object
72: - viewer - The `PetscViewer`
74: Level: beginner
76: .seealso: `PetscConvEst`, `PetscViewer`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
77: @*/
78: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
79: {
80: PetscFunctionBegin;
81: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ce, viewer));
82: PetscCall(PetscViewerASCIIPrintf(viewer, "ConvEst with %" PetscInt_FMT " levels\n", ce->Nr + 1));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*@
87: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
89: Not Collective
91: Input Parameter:
92: . ce - The `PetscConvEst` object
94: Output Parameter:
95: . solver - The solver
97: Level: intermediate
99: .seealso: `PetscConvEst`, `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
100: @*/
101: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
102: {
103: PetscFunctionBegin;
105: PetscAssertPointer(solver, 2);
106: *solver = ce->solver;
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: /*@
111: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
113: Not Collective
115: Input Parameters:
116: + ce - The `PetscConvEst` object
117: - solver - The solver, must be a `KSP`, `SNES`, or `TS` object with an attached `DM`/`DS`, that can compute an exact solution
119: Level: intermediate
121: .seealso: `PetscConvEst`, `PetscConvEstGetSNES()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
122: @*/
123: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
124: {
125: PetscFunctionBegin;
128: ce->solver = solver;
129: PetscUseTypeMethod(ce, setsolver, solver);
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*@
134: PetscConvEstSetUp - After the solver is specified, create data structures needed for estimating convergence
136: Collective
138: Input Parameter:
139: . ce - The `PetscConvEst` object
141: Level: beginner
143: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
144: @*/
145: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
146: {
147: PetscInt Nf, f, Nds, s;
149: PetscFunctionBegin;
150: PetscCall(DMGetNumFields(ce->idm, &Nf));
151: ce->Nf = PetscMax(Nf, 1);
152: PetscCall(PetscMalloc2((ce->Nr + 1) * ce->Nf, &ce->dofs, (ce->Nr + 1) * ce->Nf, &ce->errors));
153: PetscCall(PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs));
154: for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private;
155: PetscCall(DMGetNumDS(ce->idm, &Nds));
156: for (s = 0; s < Nds; ++s) {
157: PetscDS ds;
158: DMLabel label;
159: IS fieldIS;
160: const PetscInt *fields;
161: PetscInt dsNf;
163: PetscCall(DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds, NULL));
164: PetscCall(PetscDSGetNumFields(ds, &dsNf));
165: if (fieldIS) PetscCall(ISGetIndices(fieldIS, &fields));
166: for (f = 0; f < dsNf; ++f) {
167: const PetscInt field = fields[f];
168: PetscCall(PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]));
169: }
170: if (fieldIS) PetscCall(ISRestoreIndices(fieldIS, &fields));
171: }
172: for (f = 0; f < Nf; ++f) PetscCheck(ce->exactSol[f], PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "DS must contain exact solution functions in order to estimate convergence, missing for field %" PetscInt_FMT, f);
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: PetscConvEstComputeInitialGuess - Fill `u` with the initial guess to use on refinement level `r` of a convergence-estimation run
179: Collective
181: Input Parameters:
182: + ce - the `PetscConvEst` object
183: . r - the refinement level
184: - dm - the `DM` on which `u` is defined (may be `NULL`)
186: Output Parameter:
187: . u - the initial-guess vector
189: Level: developer
191: .seealso: `PetscConvEst`, `PetscConvEstComputeError()`, `PetscConvEstGetConvRate()`
192: @*/
193: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
194: {
195: PetscFunctionBegin;
199: PetscUseTypeMethod(ce, initguess, r, dm, u);
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PetscConvEstComputeError - Compute per-field discretization errors of `u` on refinement level `r`
206: Collective
208: Input Parameters:
209: + ce - the `PetscConvEst` object
210: . r - the refinement level
211: . dm - the `DM` on which `u` is defined (may be `NULL`)
212: - u - the computed solution
214: Output Parameter:
215: . errors - array of length `Nf` (number of fields in the DS) filled with the error in each field
217: Level: developer
219: .seealso: `PetscConvEst`, `PetscConvEstComputeInitialGuess()`, `PetscConvEstGetConvRate()`
220: @*/
221: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
222: {
223: PetscFunctionBegin;
227: PetscAssertPointer(errors, 5);
228: PetscUseTypeMethod(ce, computeerror, r, dm, u, errors);
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*@
233: PetscConvEstMonitorDefault - Monitors the convergence estimation loop
235: Collective
237: Input Parameters:
238: + ce - The `PetscConvEst` object
239: - r - The refinement level
241: Options Database Key:
242: . -convest_monitor - Activate the monitor
244: Level: intermediate
246: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`, `SNESSolve()`, `TSSolve()`
247: @*/
248: PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
249: {
250: MPI_Comm comm;
251: PetscInt f;
253: PetscFunctionBegin;
254: if (ce->monitor) {
255: PetscInt *dofs = &ce->dofs[r * ce->Nf];
256: PetscReal *errors = &ce->errors[r * ce->Nf];
258: PetscCall(PetscObjectGetComm((PetscObject)ce, &comm));
259: PetscCall(PetscPrintf(comm, "N: "));
260: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
261: for (f = 0; f < ce->Nf; ++f) {
262: if (f > 0) PetscCall(PetscPrintf(comm, ", "));
263: PetscCall(PetscPrintf(comm, "%7" PetscInt_FMT, dofs[f]));
264: }
265: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
266: PetscCall(PetscPrintf(comm, " "));
267: PetscCall(PetscPrintf(comm, "L_2 Error: "));
268: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
269: for (f = 0; f < ce->Nf; ++f) {
270: if (f > 0) PetscCall(PetscPrintf(comm, ", "));
271: if (errors[f] < 1.0e-11) PetscCall(PetscPrintf(comm, "< 1e-11"));
272: else PetscCall(PetscPrintf(comm, "%g", (double)errors[f]));
273: }
274: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
275: PetscCall(PetscPrintf(comm, "\n"));
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
281: {
282: PetscClassId id;
284: PetscFunctionBegin;
285: PetscCall(PetscObjectGetClassId(ce->solver, &id));
286: PetscCheck(id == SNES_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES");
287: PetscCall(SNESGetDM((SNES)ce->solver, &ce->idm));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
292: {
293: PetscFunctionBegin;
294: PetscCall(DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
299: {
300: const char *prefix;
301: PetscBool errorView = PETSC_FALSE;
303: PetscFunctionBegin;
304: PetscCall(DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors));
305: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ce, &prefix));
306: PetscCall(PetscOptionsHasName(NULL, prefix, "-convest_error_view", &errorView));
307: if (errorView) {
308: DM dmError;
309: PetscFE feError, fe;
310: PetscQuadrature quad;
311: Vec e;
312: PetscDS ds;
313: PetscSimplePointFn **funcs;
314: void **ctxs;
315: PetscInt dim, Nf;
317: PetscCall(DMGetDimension(dm, &dim));
318: PetscCall(DMGetDS(dm, &ds));
319: PetscCall(PetscDSGetNumFields(ds, &Nf));
320: PetscCall(PetscMalloc2(Nf, &funcs, Nf, &ctxs));
321: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &funcs[f], &ctxs[f]));
322: PetscCall(DMClone(dm, &dmError));
323: PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 0, PETSC_DETERMINE, &feError));
324: PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
325: PetscCall(PetscFEGetQuadrature(fe, &quad));
326: PetscCall(PetscFESetQuadrature(feError, quad));
327: PetscCall(DMSetField(dmError, 0, NULL, (PetscObject)feError));
328: PetscCall(PetscFEDestroy(&feError));
329: PetscCall(DMCreateDS(dmError));
330: PetscCall(DMGetGlobalVector(dmError, &e));
331: PetscCall(PetscObjectSetName((PetscObject)e, "Error"));
332: PetscCall(DMPlexComputeL2DiffVec(dm, 0., funcs, ctxs, u, e));
333: PetscCall(VecViewFromOptions(e, NULL, "-convest_error_view"));
334: PetscCall(DMRestoreGlobalVector(dmError, &e));
335: PetscCall(DMDestroy(&dmError));
336: PetscCall(PetscFree2(funcs, ctxs));
337: }
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: static PetscErrorCode PetscConvEstSetJacobianNullSpace_Private(PetscConvEst ce, SNES snes)
342: {
343: DM dm;
344: PetscInt f;
346: PetscFunctionBegin;
347: PetscCall(SNESGetDM(snes, &dm));
348: for (f = 0; f < ce->Nf; ++f) {
349: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
351: PetscCall(DMGetNullSpaceConstructor(dm, f, &nspconstr));
352: if (nspconstr) {
353: MatNullSpace nullsp;
354: Mat J;
356: PetscCall((*nspconstr)(dm, f, f, &nullsp));
357: PetscCall(SNESSetUp(snes));
358: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
359: PetscCall(MatSetNullSpace(J, nullsp));
360: PetscCall(MatNullSpaceDestroy(&nullsp));
361: break;
362: }
363: }
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
368: {
369: SNES snes = (SNES)ce->solver;
370: DM *dm;
371: PetscObject disc;
372: PetscReal *x, *y, slope, intercept;
373: PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
374: void *ctx;
376: PetscFunctionBegin;
377: PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
378: PetscCall(DMGetDimension(ce->idm, &dim));
379: PetscCall(DMGetApplicationContext(ce->idm, &ctx));
380: PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
381: PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
382: PetscCall(PetscMalloc1(Nr + 1, &dm));
383: /* Loop over meshes */
384: dm[0] = ce->idm;
385: for (r = 0; r <= Nr; ++r) {
386: Vec u;
387: PetscLogStage stage;
388: char stageName[PETSC_MAX_PATH_LEN];
389: const char *dmname, *uname;
391: PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
392: PetscCall(PetscLogStageGetId(stageName, &stage));
393: if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
394: PetscCall(PetscLogStagePush(stage));
395: if (r > 0) {
396: if (!ce->noRefine) {
397: PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
398: PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
399: } else {
400: DM cdm, rcdm;
402: PetscCall(DMClone(dm[r - 1], &dm[r]));
403: PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
404: PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
405: PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
406: PetscCall(DMCopyDisc(cdm, rcdm));
407: }
408: PetscCall(DMCopyTransform(ce->idm, dm[r]));
409: PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
410: PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
411: for (f = 0; f < ce->Nf; ++f) {
412: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
414: PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
415: PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
416: }
417: }
418: PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
419: /* Create solution */
420: PetscCall(DMCreateGlobalVector(dm[r], &u));
421: PetscCall(DMGetField(dm[r], 0, NULL, &disc));
422: PetscCall(PetscObjectGetName(disc, &uname));
423: PetscCall(PetscObjectSetName((PetscObject)u, uname));
424: /* Setup solver */
425: PetscCall(SNESReset(snes));
426: PetscCall(SNESSetDM(snes, dm[r]));
427: PetscCall(DMPlexSetSNESLocalFEM(dm[r], PETSC_FALSE, ctx));
428: PetscCall(DMPlexSetSNESVariableBounds(dm[r], snes));
429: PetscCall(SNESSetFromOptions(snes));
430: /* Set nullspace for Jacobian */
431: PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
432: /* Create initial guess */
433: PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
434: PetscCall(SNESSolve(snes, NULL, u));
435: PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
436: PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * ce->Nf]));
437: PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
438: for (f = 0; f < ce->Nf; ++f) {
439: PetscSection s, fs;
440: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
441: PetscCall(DMGetLocalSection(dm[r], &s));
442: PetscCall(PetscSectionGetField(s, f, &fs));
443: PetscCall(PetscSectionGetConstrainedStorageSize(fs, &ce->dofs[r * ce->Nf + f]));
444: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ce->dofs[r * ce->Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
445: PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * ce->Nf + f]));
446: PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * ce->Nf + f]));
447: }
448: /* Monitor */
449: PetscCall(PetscConvEstMonitorDefault(ce, r));
450: if (!r) {
451: /* PCReset() does not wipe out the level structure */
452: KSP ksp;
453: PC pc;
455: PetscCall(SNESGetKSP(snes, &ksp));
456: PetscCall(KSPGetPC(ksp, &pc));
457: PetscCall(PCMGGetLevels(pc, &oldnlev));
458: }
459: /* Cleanup */
460: PetscCall(VecDestroy(&u));
461: PetscCall(PetscLogStagePop());
462: }
463: for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
464: /* Fit convergence rate */
465: PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
466: for (f = 0; f < ce->Nf; ++f) {
467: for (r = 0; r <= Nr; ++r) {
468: x[r] = PetscLog10Real(ce->dofs[r * ce->Nf + f]);
469: y[r] = PetscLog10Real(ce->errors[r * ce->Nf + f]);
470: }
471: PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
472: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
473: alpha[f] = -slope * dim;
474: }
475: PetscCall(PetscFree2(x, y));
476: PetscCall(PetscFree(dm));
477: /* Restore solver */
478: PetscCall(SNESReset(snes));
479: {
480: /* PCReset() does not wipe out the level structure */
481: KSP ksp;
482: PC pc;
484: PetscCall(SNESGetKSP(snes, &ksp));
485: PetscCall(KSPGetPC(ksp, &pc));
486: PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
487: PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
488: }
489: PetscCall(SNESSetDM(snes, ce->idm));
490: PetscCall(DMPlexSetSNESLocalFEM(ce->idm, PETSC_FALSE, ctx));
491: PetscCall(DMPlexSetSNESVariableBounds(ce->idm, snes));
492: PetscCall(SNESSetFromOptions(snes));
493: PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*@
498: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
500: Not Collective
502: Input Parameter:
503: . ce - The `PetscConvEst` object
505: Output Parameter:
506: . alpha - The convergence rate for each field
508: Options Database Keys:
509: + -snes_convergence_estimate - Execute convergence estimation inside `SNESSolve()` and print out the rate
510: - -ts_convergence_estimate - Execute convergence estimation inside `TSSolve()` and print out the rate
512: Level: intermediate
514: Notes:
515: The convergence rate alpha is defined by
517: $$
518: || u_\Delta - u_{exact} || < C \Delta^\alpha
519: $$
521: where $u_{\Delta} $ is the discrete solution, and $\Delta$ is a measure of the discretization size. We usually use $h$ for the
522: spatial resolution and $\Delta t $ for the temporal resolution.
524: We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
525: based upon the exact solution in the `PetscDS`, and then fit the result to our model above using linear regression.
527: .seealso: `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `SNESSolve()`, `TSSolve()`
528: @*/
529: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
530: {
531: PetscInt f;
533: PetscFunctionBegin;
534: if (ce->event < 0) PetscCall(PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event));
535: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
536: PetscUseTypeMethod(ce, getconvrate, alpha);
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*@
541: PetscConvEstRateView - Displays the convergence rate obtained from `PetscConvEstGetConvRate()` using a `PetscViewer`
543: Collective
545: Input Parameters:
546: + ce - iterative context obtained from `SNESCreate()`
547: . alpha - the convergence rate for each field
548: - viewer - the viewer to display the reason
550: Options Database Key:
551: . -snes_convergence_estimate - print the convergence rate
553: Level: developer
555: .seealso: `PetscConvEst`, `PetscConvEstGetConvRate()`
556: @*/
557: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
558: {
559: PetscBool isAscii;
561: PetscFunctionBegin;
562: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
563: if (isAscii) {
564: PetscInt Nf = ce->Nf, f;
566: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ce)->tablevel));
567: PetscCall(PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: "));
568: if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "["));
569: for (f = 0; f < Nf; ++f) {
570: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
571: PetscCall(PetscViewerASCIIPrintf(viewer, "%#.2g", (double)alpha[f]));
572: }
573: if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "]"));
574: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
575: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ce)->tablevel));
576: }
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: PetscConvEstCreate - Create a `PetscConvEst` object. This is used to study the convergence rate of approximations on grids to a continuum solution
583: Collective
585: Input Parameter:
586: . comm - The communicator for the `PetscConvEst` object
588: Output Parameter:
589: . ce - The `PetscConvEst` object
591: Level: beginner
593: .seealso: `PetscConvEst`, `PetscConvEstDestroy()`, `PetscConvEstGetConvRate()`, `DMAdaptorCreate()`, `DMAdaptor`
594: @*/
595: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
596: {
597: PetscFunctionBegin;
598: PetscAssertPointer(ce, 2);
599: PetscCall(PetscSysInitializePackage());
600: PetscCall(PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView));
601: (*ce)->monitor = PETSC_FALSE;
602: (*ce)->r = 2.0;
603: (*ce)->Nr = 4;
604: (*ce)->event = -1;
605: (*ce)->ops->setsolver = PetscConvEstSetSNES_Private;
606: (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private;
607: (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
608: (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private;
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }