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