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