1 #include "petscsys.h"
2 #include <petscconvest.h> /*I "petscconvest.h" I*/
3 #include <petscdmplex.h>
4 #include <petscds.h>
5
6 #include <petsc/private/petscconvestimpl.h>
7
zero_private(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)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 }
14
15 /*@
16 PetscConvEstDestroy - Destroys a PETSc convergence estimator `PetscConvEst` object
17
18 Collective
19
20 Input Parameter:
21 . ce - The `PetscConvEst` object
22
23 Level: beginner
24
25 .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
26 @*/
PetscConvEstDestroy(PetscConvEst * ce)27 PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
28 {
29 PetscFunctionBegin;
30 if (!*ce) PetscFunctionReturn(PETSC_SUCCESS);
31 PetscValidHeaderSpecific(*ce, PETSC_OBJECT_CLASSID, 1);
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 }
41
42 /*@
43 PetscConvEstSetFromOptions - Sets a convergence estimator `PetscConvEst` object based on values in the options database
44
45 Collective
46
47 Input Parameter:
48 . ce - The `PetscConvEst` object
49
50 Level: beginner
51
52 .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
53 @*/
PetscConvEstSetFromOptions(PetscConvEst ce)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 }
65
66 /*@
67 PetscConvEstView - Views a `PetscConvEst` object
68
69 Collective
70
71 Input Parameters:
72 + ce - The `PetscConvEst` object
73 - viewer - The `PetscViewer`
74
75 Level: beginner
76
77 .seealso: `PetscConvEst`, `PetscViewer`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
78 @*/
PetscConvEstView(PetscConvEst ce,PetscViewer viewer)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 }
86
87 /*@
88 PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
89
90 Not Collective
91
92 Input Parameter:
93 . ce - The `PetscConvEst` object
94
95 Output Parameter:
96 . solver - The solver
97
98 Level: intermediate
99
100 .seealso: `PetscConvEst`, `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
101 @*/
PetscConvEstGetSolver(PetscConvEst ce,PetscObject * solver)102 PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
103 {
104 PetscFunctionBegin;
105 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
106 PetscAssertPointer(solver, 2);
107 *solver = ce->solver;
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
111 /*@
112 PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
113
114 Not Collective
115
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
119
120 Level: intermediate
121
122 .seealso: `PetscConvEst`, `PetscConvEstGetSNES()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
123 @*/
PetscConvEstSetSolver(PetscConvEst ce,PetscObject solver)124 PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
125 {
126 PetscFunctionBegin;
127 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
128 PetscValidHeader(solver, 2);
129 ce->solver = solver;
130 PetscUseTypeMethod(ce, setsolver, solver);
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
134 /*@
135 PetscConvEstSetUp - After the solver is specified, create data structures needed for estimating convergence
136
137 Collective
138
139 Input Parameter:
140 . ce - The `PetscConvEst` object
141
142 Level: beginner
143
144 .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
145 @*/
PetscConvEstSetUp(PetscConvEst ce)146 PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
147 {
148 PetscInt Nf, f, Nds, s;
149
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;
163
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 }
176
PetscConvEstComputeInitialGuess(PetscConvEst ce,PetscInt r,DM dm,Vec u)177 PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
178 {
179 PetscFunctionBegin;
180 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
181 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 3);
182 PetscValidHeaderSpecific(u, VEC_CLASSID, 4);
183 PetscUseTypeMethod(ce, initguess, r, dm, u);
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
PetscConvEstComputeError(PetscConvEst ce,PetscInt r,DM dm,Vec u,PetscReal errors[])187 PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
188 {
189 PetscFunctionBegin;
190 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
191 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 3);
192 PetscValidHeaderSpecific(u, VEC_CLASSID, 4);
193 PetscAssertPointer(errors, 5);
194 PetscUseTypeMethod(ce, computeerror, r, dm, u, errors);
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
198 /*@
199 PetscConvEstMonitorDefault - Monitors the convergence estimation loop
200
201 Collective
202
203 Input Parameters:
204 + ce - The `PetscConvEst` object
205 - r - The refinement level
206
207 Options Database Key:
208 . -convest_monitor - Activate the monitor
209
210 Level: intermediate
211
212 .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`, `SNESSolve()`, `TSSolve()`
213 @*/
PetscConvEstMonitorDefault(PetscConvEst ce,PetscInt r)214 PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
215 {
216 MPI_Comm comm;
217 PetscInt f;
218
219 PetscFunctionBegin;
220 if (ce->monitor) {
221 PetscInt *dofs = &ce->dofs[r * ce->Nf];
222 PetscReal *errors = &ce->errors[r * ce->Nf];
223
224 PetscCall(PetscObjectGetComm((PetscObject)ce, &comm));
225 PetscCall(PetscPrintf(comm, "N: "));
226 if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
227 for (f = 0; f < ce->Nf; ++f) {
228 if (f > 0) PetscCall(PetscPrintf(comm, ", "));
229 PetscCall(PetscPrintf(comm, "%7" PetscInt_FMT, dofs[f]));
230 }
231 if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
232 PetscCall(PetscPrintf(comm, " "));
233 PetscCall(PetscPrintf(comm, "L_2 Error: "));
234 if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "["));
235 for (f = 0; f < ce->Nf; ++f) {
236 if (f > 0) PetscCall(PetscPrintf(comm, ", "));
237 if (errors[f] < 1.0e-11) PetscCall(PetscPrintf(comm, "< 1e-11"));
238 else PetscCall(PetscPrintf(comm, "%g", (double)errors[f]));
239 }
240 if (ce->Nf > 1) PetscCall(PetscPrintf(comm, "]"));
241 PetscCall(PetscPrintf(comm, "\n"));
242 }
243 PetscFunctionReturn(PETSC_SUCCESS);
244 }
245
PetscConvEstSetSNES_Private(PetscConvEst ce,PetscObject solver)246 static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
247 {
248 PetscClassId id;
249
250 PetscFunctionBegin;
251 PetscCall(PetscObjectGetClassId(ce->solver, &id));
252 PetscCheck(id == SNES_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES");
253 PetscCall(SNESGetDM((SNES)ce->solver, &ce->idm));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
PetscConvEstInitGuessSNES_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u)257 static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
258 {
259 PetscFunctionBegin;
260 PetscCall(DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u));
261 PetscFunctionReturn(PETSC_SUCCESS);
262 }
263
PetscConvEstComputeErrorSNES_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u,PetscReal errors[])264 static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
265 {
266 const char *prefix;
267 PetscBool errorView = PETSC_FALSE;
268
269 PetscFunctionBegin;
270 PetscCall(DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors));
271 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ce, &prefix));
272 PetscCall(PetscOptionsHasName(NULL, prefix, "-convest_error_view", &errorView));
273 if (errorView) {
274 DM dmError;
275 PetscFE feError, fe;
276 PetscQuadrature quad;
277 Vec e;
278 PetscDS ds;
279 PetscSimplePointFn **funcs;
280 void **ctxs;
281 PetscInt dim, Nf;
282
283 PetscCall(DMGetDimension(dm, &dim));
284 PetscCall(DMGetDS(dm, &ds));
285 PetscCall(PetscDSGetNumFields(ds, &Nf));
286 PetscCall(PetscMalloc2(Nf, &funcs, Nf, &ctxs));
287 for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &funcs[f], &ctxs[f]));
288 PetscCall(DMClone(dm, &dmError));
289 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 0, PETSC_DETERMINE, &feError));
290 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
291 PetscCall(PetscFEGetQuadrature(fe, &quad));
292 PetscCall(PetscFESetQuadrature(feError, quad));
293 PetscCall(DMSetField(dmError, 0, NULL, (PetscObject)feError));
294 PetscCall(PetscFEDestroy(&feError));
295 PetscCall(DMCreateDS(dmError));
296 PetscCall(DMGetGlobalVector(dmError, &e));
297 PetscCall(PetscObjectSetName((PetscObject)e, "Error"));
298 PetscCall(DMPlexComputeL2DiffVec(dm, 0., funcs, ctxs, u, e));
299 PetscCall(VecViewFromOptions(e, NULL, "-convest_error_view"));
300 PetscCall(DMRestoreGlobalVector(dmError, &e));
301 PetscCall(DMDestroy(&dmError));
302 PetscCall(PetscFree2(funcs, ctxs));
303 }
304 PetscFunctionReturn(PETSC_SUCCESS);
305 }
306
PetscConvEstSetJacobianNullSpace_Private(PetscConvEst ce,SNES snes)307 static PetscErrorCode PetscConvEstSetJacobianNullSpace_Private(PetscConvEst ce, SNES snes)
308 {
309 DM dm;
310 PetscInt f;
311
312 PetscFunctionBegin;
313 PetscCall(SNESGetDM(snes, &dm));
314 for (f = 0; f < ce->Nf; ++f) {
315 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
316
317 PetscCall(DMGetNullSpaceConstructor(dm, f, &nspconstr));
318 if (nspconstr) {
319 MatNullSpace nullsp;
320 Mat J;
321
322 PetscCall((*nspconstr)(dm, f, f, &nullsp));
323 PetscCall(SNESSetUp(snes));
324 PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
325 PetscCall(MatSetNullSpace(J, nullsp));
326 PetscCall(MatNullSpaceDestroy(&nullsp));
327 break;
328 }
329 }
330 PetscFunctionReturn(PETSC_SUCCESS);
331 }
332
PetscConvEstGetConvRateSNES_Private(PetscConvEst ce,PetscReal alpha[])333 static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
334 {
335 SNES snes = (SNES)ce->solver;
336 DM *dm;
337 PetscObject disc;
338 PetscReal *x, *y, slope, intercept;
339 PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
340 void *ctx;
341
342 PetscFunctionBegin;
343 PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
344 PetscCall(DMGetDimension(ce->idm, &dim));
345 PetscCall(DMGetApplicationContext(ce->idm, &ctx));
346 PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
347 PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
348 PetscCall(PetscMalloc1(Nr + 1, &dm));
349 /* Loop over meshes */
350 dm[0] = ce->idm;
351 for (r = 0; r <= Nr; ++r) {
352 Vec u;
353 PetscLogStage stage;
354 char stageName[PETSC_MAX_PATH_LEN];
355 const char *dmname, *uname;
356
357 PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
358 PetscCall(PetscLogStageGetId(stageName, &stage));
359 if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
360 PetscCall(PetscLogStagePush(stage));
361 if (r > 0) {
362 if (!ce->noRefine) {
363 PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
364 PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
365 } else {
366 DM cdm, rcdm;
367
368 PetscCall(DMClone(dm[r - 1], &dm[r]));
369 PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
370 PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
371 PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
372 PetscCall(DMCopyDisc(cdm, rcdm));
373 }
374 PetscCall(DMCopyTransform(ce->idm, dm[r]));
375 PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
376 PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
377 for (f = 0; f < ce->Nf; ++f) {
378 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
379
380 PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
381 PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
382 }
383 }
384 PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
385 /* Create solution */
386 PetscCall(DMCreateGlobalVector(dm[r], &u));
387 PetscCall(DMGetField(dm[r], 0, NULL, &disc));
388 PetscCall(PetscObjectGetName(disc, &uname));
389 PetscCall(PetscObjectSetName((PetscObject)u, uname));
390 /* Setup solver */
391 PetscCall(SNESReset(snes));
392 PetscCall(SNESSetDM(snes, dm[r]));
393 PetscCall(DMPlexSetSNESLocalFEM(dm[r], PETSC_FALSE, ctx));
394 PetscCall(DMPlexSetSNESVariableBounds(dm[r], snes));
395 PetscCall(SNESSetFromOptions(snes));
396 /* Set nullspace for Jacobian */
397 PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
398 /* Create initial guess */
399 PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
400 PetscCall(SNESSolve(snes, NULL, u));
401 PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
402 PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * ce->Nf]));
403 PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
404 for (f = 0; f < ce->Nf; ++f) {
405 PetscSection s, fs;
406 PetscInt lsize;
407
408 /* Could use DMGetOutputDM() to add in Dirichlet dofs */
409 PetscCall(DMGetLocalSection(dm[r], &s));
410 PetscCall(PetscSectionGetField(s, f, &fs));
411 PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize));
412 PetscCallMPI(MPIU_Allreduce(&lsize, &ce->dofs[r * ce->Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
413 PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * ce->Nf + f]));
414 PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * ce->Nf + f]));
415 }
416 /* Monitor */
417 PetscCall(PetscConvEstMonitorDefault(ce, r));
418 if (!r) {
419 /* PCReset() does not wipe out the level structure */
420 KSP ksp;
421 PC pc;
422
423 PetscCall(SNESGetKSP(snes, &ksp));
424 PetscCall(KSPGetPC(ksp, &pc));
425 PetscCall(PCMGGetLevels(pc, &oldnlev));
426 }
427 /* Cleanup */
428 PetscCall(VecDestroy(&u));
429 PetscCall(PetscLogStagePop());
430 }
431 for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
432 /* Fit convergence rate */
433 PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
434 for (f = 0; f < ce->Nf; ++f) {
435 for (r = 0; r <= Nr; ++r) {
436 x[r] = PetscLog10Real(ce->dofs[r * ce->Nf + f]);
437 y[r] = PetscLog10Real(ce->errors[r * ce->Nf + f]);
438 }
439 PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
440 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
441 alpha[f] = -slope * dim;
442 }
443 PetscCall(PetscFree2(x, y));
444 PetscCall(PetscFree(dm));
445 /* Restore solver */
446 PetscCall(SNESReset(snes));
447 {
448 /* PCReset() does not wipe out the level structure */
449 KSP ksp;
450 PC pc;
451
452 PetscCall(SNESGetKSP(snes, &ksp));
453 PetscCall(KSPGetPC(ksp, &pc));
454 PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
455 PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
456 }
457 PetscCall(SNESSetDM(snes, ce->idm));
458 PetscCall(DMPlexSetSNESLocalFEM(ce->idm, PETSC_FALSE, ctx));
459 PetscCall(DMPlexSetSNESVariableBounds(ce->idm, snes));
460 PetscCall(SNESSetFromOptions(snes));
461 PetscCall(PetscConvEstSetJacobianNullSpace_Private(ce, snes));
462 PetscFunctionReturn(PETSC_SUCCESS);
463 }
464
465 /*@
466 PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
467
468 Not Collective
469
470 Input Parameter:
471 . ce - The `PetscConvEst` object
472
473 Output Parameter:
474 . alpha - The convergence rate for each field
475
476 Options Database Keys:
477 + -snes_convergence_estimate - Execute convergence estimation inside `SNESSolve()` and print out the rate
478 - -ts_convergence_estimate - Execute convergence estimation inside `TSSolve()` and print out the rate
479
480 Level: intermediate
481
482 Notes:
483 The convergence rate alpha is defined by
484
485 $$
486 || u_\Delta - u_{exact} || < C \Delta^\alpha
487 $$
488
489 where $u_{\Delta} $ is the discrete solution, and $\Delta$ is a measure of the discretization size. We usually use $h$ for the
490 spatial resolution and $\Delta t $ for the temporal resolution.
491
492 We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
493 based upon the exact solution in the `PetscDS`, and then fit the result to our model above using linear regression.
494
495 .seealso: `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `SNESSolve()`, `TSSolve()`
496 @*/
PetscConvEstGetConvRate(PetscConvEst ce,PetscReal alpha[])497 PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
498 {
499 PetscInt f;
500
501 PetscFunctionBegin;
502 if (ce->event < 0) PetscCall(PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event));
503 for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
504 PetscUseTypeMethod(ce, getconvrate, alpha);
505 PetscFunctionReturn(PETSC_SUCCESS);
506 }
507
508 /*@
509 PetscConvEstRateView - Displays the convergence rate obtained from `PetscConvEstGetConvRate()` using a `PetscViewer`
510
511 Collective
512
513 Input Parameters:
514 + ce - iterative context obtained from `SNESCreate()`
515 . alpha - the convergence rate for each field
516 - viewer - the viewer to display the reason
517
518 Options Database Key:
519 . -snes_convergence_estimate - print the convergence rate
520
521 Level: developer
522
523 .seealso: `PetscConvEst`, `PetscConvEstGetConvRate()`
524 @*/
PetscConvEstRateView(PetscConvEst ce,const PetscReal alpha[],PetscViewer viewer)525 PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
526 {
527 PetscBool isAscii;
528
529 PetscFunctionBegin;
530 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
531 if (isAscii) {
532 PetscInt Nf = ce->Nf, f;
533
534 PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ce)->tablevel));
535 PetscCall(PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: "));
536 if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "["));
537 for (f = 0; f < Nf; ++f) {
538 if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
539 PetscCall(PetscViewerASCIIPrintf(viewer, "%#.2g", (double)alpha[f]));
540 }
541 if (Nf > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "]"));
542 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
543 PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ce)->tablevel));
544 }
545 PetscFunctionReturn(PETSC_SUCCESS);
546 }
547
548 /*@
549 PetscConvEstCreate - Create a `PetscConvEst` object. This is used to study the convergence rate of approximations on grids to a continuum solution
550
551 Collective
552
553 Input Parameter:
554 . comm - The communicator for the `PetscConvEst` object
555
556 Output Parameter:
557 . ce - The `PetscConvEst` object
558
559 Level: beginner
560
561 .seealso: `PetscConvEst`, `PetscConvEstDestroy()`, `PetscConvEstGetConvRate()`, `DMAdaptorCreate()`, `DMAdaptor`
562 @*/
PetscConvEstCreate(MPI_Comm comm,PetscConvEst * ce)563 PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
564 {
565 PetscFunctionBegin;
566 PetscAssertPointer(ce, 2);
567 PetscCall(PetscSysInitializePackage());
568 PetscCall(PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView));
569 (*ce)->monitor = PETSC_FALSE;
570 (*ce)->r = 2.0;
571 (*ce)->Nr = 4;
572 (*ce)->event = -1;
573 (*ce)->ops->setsolver = PetscConvEstSetSNES_Private;
574 (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private;
575 (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
576 (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private;
577 PetscFunctionReturn(PETSC_SUCCESS);
578 }
579