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