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