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