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: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
177: {
178: PetscFunctionBegin;
182: PetscUseTypeMethod(ce, initguess, r, dm, u);
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
187: {
188: PetscFunctionBegin;
192: PetscAssertPointer(errors, 5);
193: PetscUseTypeMethod(ce, computeerror, r, dm, u, errors);
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@
198: PetscConvEstMonitorDefault - Monitors the convergence estimation loop
200: Collective
202: Input Parameters:
203: + ce - The `PetscConvEst` object
204: - r - The refinement level
206: Options Database Key:
207: . -convest_monitor - Activate the monitor
209: Level: intermediate
211: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`, `SNESSolve()`, `TSSolve()`
212: @*/
213: PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
214: {
215: MPI_Comm comm;
216: PetscInt f;
218: PetscFunctionBegin;
219: if (ce->monitor) {
220: PetscInt *dofs = &ce->dofs[r * ce->Nf];
221: PetscReal *errors = &ce->errors[r * ce->Nf];
223: PetscCall(PetscObjectGetComm((PetscObject)ce, &comm));
224: PetscCall(PetscPrintf(comm, "N: "));
225: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
226: for (f = 0; f < ce->Nf; ++f) {
227: if (f > 0) PetscCall(PetscPrintf(comm, ", "));
228: PetscCall(PetscPrintf(comm, "%7" PetscInt_FMT, dofs[f]));
229: }
230: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
231: PetscCall(PetscPrintf(comm, " "));
232: PetscCall(PetscPrintf(comm, "L_2 Error: "));
233: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
234: for (f = 0; f < ce->Nf; ++f) {
235: if (f > 0) PetscCall(PetscPrintf(comm, ", "));
236: if (errors[f] < 1.0e-11) PetscCall(PetscPrintf(comm, "< 1e-11"));
237: else PetscCall(PetscPrintf(comm, "%g", (double)errors[f]));
238: }
239: if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
240: PetscCall(PetscPrintf(comm, "\n"));
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
246: {
247: PetscClassId id;
249: PetscFunctionBegin;
250: PetscCall(PetscObjectGetClassId(ce->solver, &id));
251: PetscCheck(id == SNES_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES");
252: PetscCall(SNESGetDM((SNES)ce->solver, &ce->idm));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
257: {
258: PetscFunctionBegin;
259: PetscCall(DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
264: {
265: const char *prefix;
266: PetscBool errorView = PETSC_FALSE;
268: PetscFunctionBegin;
269: PetscCall(DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors));
270: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ce, &prefix));
271: PetscCall(PetscOptionsHasName(NULL, prefix, "-convest_error_view", &errorView));
272: if (errorView) {
273: DM dmError;
274: PetscFE feError, fe;
275: PetscQuadrature quad;
276: Vec e;
277: PetscDS ds;
278: PetscSimplePointFn **funcs;
279: void **ctxs;
280: PetscInt dim, Nf;
282: PetscCall(DMGetDimension(dm, &dim));
283: PetscCall(DMGetDS(dm, &ds));
284: PetscCall(PetscDSGetNumFields(ds, &Nf));
285: PetscCall(PetscMalloc2(Nf, &funcs, Nf, &ctxs));
286: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &funcs[f], &ctxs[f]));
287: PetscCall(DMClone(dm, &dmError));
288: PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 0, PETSC_DETERMINE, &feError));
289: PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
290: PetscCall(PetscFEGetQuadrature(fe, &quad));
291: PetscCall(PetscFESetQuadrature(feError, quad));
292: PetscCall(DMSetField(dmError, 0, NULL, (PetscObject)feError));
293: PetscCall(PetscFEDestroy(&feError));
294: PetscCall(DMCreateDS(dmError));
295: PetscCall(DMGetGlobalVector(dmError, &e));
296: PetscCall(PetscObjectSetName((PetscObject)e, "Error"));
297: PetscCall(DMPlexComputeL2DiffVec(dm, 0., funcs, ctxs, u, e));
298: PetscCall(VecViewFromOptions(e, NULL, "-convest_error_view"));
299: PetscCall(DMRestoreGlobalVector(dmError, &e));
300: PetscCall(DMDestroy(&dmError));
301: PetscCall(PetscFree2(funcs, ctxs));
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode PetscConvEstSetJacobianNullSpace_Private(PetscConvEst ce, SNES snes)
307: {
308: DM dm;
309: PetscInt f;
311: PetscFunctionBegin;
312: PetscCall(SNESGetDM(snes, &dm));
313: for (f = 0; f < ce->Nf; ++f) {
314: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
316: PetscCall(DMGetNullSpaceConstructor(dm, f, &nspconstr));
317: if (nspconstr) {
318: MatNullSpace nullsp;
319: Mat J;
321: PetscCall((*nspconstr)(dm, f, f, &nullsp));
322: PetscCall(SNESSetUp(snes));
323: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
324: PetscCall(MatSetNullSpace(J, nullsp));
325: PetscCall(MatNullSpaceDestroy(&nullsp));
326: break;
327: }
328: }
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
333: {
334: SNES snes = (SNES)ce->solver;
335: DM *dm;
336: PetscObject disc;
337: PetscReal *x, *y, slope, intercept;
338: PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
339: void *ctx;
341: PetscFunctionBegin;
342: PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
343: PetscCall(DMGetDimension(ce->idm, &dim));
344: PetscCall(DMGetApplicationContext(ce->idm, &ctx));
345: PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
346: PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
347: PetscCall(PetscMalloc1(Nr + 1, &dm));
348: /* Loop over meshes */
349: dm[0] = ce->idm;
350: for (r = 0; r <= Nr; ++r) {
351: Vec u;
352: PetscLogStage stage;
353: char stageName[PETSC_MAX_PATH_LEN];
354: const char *dmname, *uname;
356: PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
357: PetscCall(PetscLogStageGetId(stageName, &stage));
358: if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
359: PetscCall(PetscLogStagePush(stage));
360: if (r > 0) {
361: if (!ce->noRefine) {
362: PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
363: PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
364: } else {
365: DM cdm, rcdm;
367: PetscCall(DMClone(dm[r - 1], &dm[r]));
368: PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
369: PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
370: PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
371: PetscCall(DMCopyDisc(cdm, rcdm));
372: }
373: PetscCall(DMCopyTransform(ce->idm, dm[r]));
374: PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
375: PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
376: for (f = 0; f < ce->Nf; ++f) {
377: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
379: PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
380: PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
381: }
382: }
383: PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
384: /* Create solution */
385: PetscCall(DMCreateGlobalVector(dm[r], &u));
386: PetscCall(DMGetField(dm[r], 0, NULL, &disc));
387: PetscCall(PetscObjectGetName(disc, &uname));
388: PetscCall(PetscObjectSetName((PetscObject)u, uname));
389: /* Setup solver */
390: PetscCall(SNESReset(snes));
391: PetscCall(SNESSetDM(snes, dm[r]));
392: PetscCall(DMPlexSetSNESLocalFEM(dm[r], PETSC_FALSE, ctx));
393: PetscCall(DMPlexSetSNESVariableBounds(dm[r], snes));
394: PetscCall(SNESSetFromOptions(snes));
395: /* Set nullspace for Jacobian */
396: PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
397: /* Create initial guess */
398: PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
399: PetscCall(SNESSolve(snes, NULL, u));
400: PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
401: PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * ce->Nf]));
402: PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
403: for (f = 0; f < ce->Nf; ++f) {
404: PetscSection s, fs;
405: PetscInt lsize;
407: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
408: PetscCall(DMGetLocalSection(dm[r], &s));
409: PetscCall(PetscSectionGetField(s, f, &fs));
410: PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize));
411: PetscCallMPI(MPIU_Allreduce(&lsize, &ce->dofs[r * ce->Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
412: PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * ce->Nf + f]));
413: PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * ce->Nf + f]));
414: }
415: /* Monitor */
416: PetscCall(PetscConvEstMonitorDefault(ce, r));
417: if (!r) {
418: /* PCReset() does not wipe out the level structure */
419: KSP ksp;
420: PC pc;
422: PetscCall(SNESGetKSP(snes, &ksp));
423: PetscCall(KSPGetPC(ksp, &pc));
424: PetscCall(PCMGGetLevels(pc, &oldnlev));
425: }
426: /* Cleanup */
427: PetscCall(VecDestroy(&u));
428: PetscCall(PetscLogStagePop());
429: }
430: for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
431: /* Fit convergence rate */
432: PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
433: for (f = 0; f < ce->Nf; ++f) {
434: for (r = 0; r <= Nr; ++r) {
435: x[r] = PetscLog10Real(ce->dofs[r * ce->Nf + f]);
436: y[r] = PetscLog10Real(ce->errors[r * ce->Nf + f]);
437: }
438: PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
439: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
440: alpha[f] = -slope * dim;
441: }
442: PetscCall(PetscFree2(x, y));
443: PetscCall(PetscFree(dm));
444: /* Restore solver */
445: PetscCall(SNESReset(snes));
446: {
447: /* PCReset() does not wipe out the level structure */
448: KSP ksp;
449: PC pc;
451: PetscCall(SNESGetKSP(snes, &ksp));
452: PetscCall(KSPGetPC(ksp, &pc));
453: PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
454: PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
455: }
456: PetscCall(SNESSetDM(snes, ce->idm));
457: PetscCall(DMPlexSetSNESLocalFEM(ce->idm, PETSC_FALSE, ctx));
458: PetscCall(DMPlexSetSNESVariableBounds(ce->idm, snes));
459: PetscCall(SNESSetFromOptions(snes));
460: PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
467: Not Collective
469: Input Parameter:
470: . ce - The `PetscConvEst` object
472: Output Parameter:
473: . alpha - The convergence rate for each field
475: Options Database Keys:
476: + -snes_convergence_estimate - Execute convergence estimation inside `SNESSolve()` and print out the rate
477: - -ts_convergence_estimate - Execute convergence estimation inside `TSSolve()` and print out the rate
479: Level: intermediate
481: Notes:
482: The convergence rate alpha is defined by
484: $$
485: || u_\Delta - u_{exact} || < C \Delta^\alpha
486: $$
488: where $u_{\Delta} $ is the discrete solution, and $\Delta$ is a measure of the discretization size. We usually use $h$ for the
489: spatial resolution and $\Delta t $ for the temporal resolution.
491: We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
492: based upon the exact solution in the `PetscDS`, and then fit the result to our model above using linear regression.
494: .seealso: `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `SNESSolve()`, `TSSolve()`
495: @*/
496: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
497: {
498: PetscInt f;
500: PetscFunctionBegin;
501: if (ce->event < 0) PetscCall(PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event));
502: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
503: PetscUseTypeMethod(ce, getconvrate, alpha);
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@
508: PetscConvEstRateView - Displays the convergence rate obtained from `PetscConvEstGetConvRate()` using a `PetscViewer`
510: Collective
512: Input Parameters:
513: + ce - iterative context obtained from `SNESCreate()`
514: . alpha - the convergence rate for each field
515: - viewer - the viewer to display the reason
517: Options Database Key:
518: . -snes_convergence_estimate - print the convergence rate
520: Level: developer
522: .seealso: `PetscConvEst`, `PetscConvEstGetConvRate()`
523: @*/
524: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
525: {
526: PetscBool isAscii;
528: PetscFunctionBegin;
529: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
530: if (isAscii) {
531: PetscInt Nf = ce->Nf, f;
533: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ce)->tablevel));
534: PetscCall(PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: "));
535: if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "["));
536: for (f = 0; f < Nf; ++f) {
537: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
538: PetscCall(PetscViewerASCIIPrintf(viewer, "%#.2g", (double)alpha[f]));
539: }
540: if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "]"));
541: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
542: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ce)->tablevel));
543: }
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: PetscConvEstCreate - Create a `PetscConvEst` object. This is used to study the convergence rate of approximations on grids to a continuum solution
550: Collective
552: Input Parameter:
553: . comm - The communicator for the `PetscConvEst` object
555: Output Parameter:
556: . ce - The `PetscConvEst` object
558: Level: beginner
560: .seealso: `PetscConvEst`, `PetscConvEstDestroy()`, `PetscConvEstGetConvRate()`, `DMAdaptorCreate()`, `DMAdaptor`
561: @*/
562: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
563: {
564: PetscFunctionBegin;
565: PetscAssertPointer(ce, 2);
566: PetscCall(PetscSysInitializePackage());
567: PetscCall(PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView));
568: (*ce)->monitor = PETSC_FALSE;
569: (*ce)->r = 2.0;
570: (*ce)->Nr = 4;
571: (*ce)->event = -1;
572: (*ce)->ops->setsolver = PetscConvEstSetSNES_Private;
573: (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private;
574: (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
575: (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }