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: PetscOptionsBool("-convest_no_refine", "Debugging flag to run on the same mesh each time", "PetscConvEst", ce->noRefine, &ce->noRefine, NULL);
65: PetscOptionsEnd();
66: return(0);
67: }
69: /*@
70: PetscConvEstView - Views a PetscConvEst object
72: Collective on PetscConvEst
74: Input Parameters:
75: + ce - The PetscConvEst object
76: - viewer - The PetscViewer object
78: Level: beginner
80: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
81: @*/
82: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
83: {
87: PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
88: PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
89: return(0);
90: }
92: /*@
93: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
95: Not collective
97: Input Parameter:
98: . ce - The PetscConvEst object
100: Output Parameter:
101: . solver - The solver
103: Level: intermediate
105: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
106: @*/
107: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
108: {
112: *solver = ce->solver;
113: return(0);
114: }
116: /*@
117: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
119: Not collective
121: Input Parameters:
122: + ce - The PetscConvEst object
123: - solver - The solver
125: Level: intermediate
127: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
129: .seealso: PetscConvEstGetSNES(), PetscConvEstCreate(), PetscConvEstGetConvRate()
130: @*/
131: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
132: {
138: ce->solver = solver;
139: (*ce->ops->setsolver)(ce, solver);
140: return(0);
141: }
143: /*@
144: PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence
146: Collective on PetscConvEst
148: Input Parameters:
149: . ce - The PetscConvEst object
151: Level: beginner
153: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
154: @*/
155: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
156: {
157: PetscInt Nf, f, Nds, s;
161: DMGetNumFields(ce->idm, &Nf);
162: ce->Nf = PetscMax(Nf, 1);
163: PetscMalloc2((ce->Nr+1)*ce->Nf, &ce->dofs, (ce->Nr+1)*ce->Nf, &ce->errors);
164: PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);
165: for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private;
166: DMGetNumDS(ce->idm, &Nds);
167: for (s = 0; s < Nds; ++s) {
168: PetscDS ds;
169: DMLabel label;
170: IS fieldIS;
171: const PetscInt *fields;
172: PetscInt dsNf;
174: DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds);
175: PetscDSGetNumFields(ds, &dsNf);
176: if (fieldIS) {ISGetIndices(fieldIS, &fields);}
177: for (f = 0; f < dsNf; ++f) {
178: const PetscInt field = fields[f];
179: PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]);
180: }
181: if (fieldIS) {ISRestoreIndices(fieldIS, &fields);}
182: }
183: for (f = 0; f < Nf; ++f) {
184: 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);
185: }
186: return(0);
187: }
189: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
190: {
197: (*ce->ops->initguess)(ce, r, dm, u);
198: return(0);
199: }
201: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
202: {
210: (*ce->ops->computeerror)(ce, r, dm, u, errors);
211: return(0);
212: }
214: /*@
215: PetscConvEstMonitorDefault - Monitors the convergence estimation loop
217: Collective on PetscConvEst
219: Input Parameters:
220: + ce - The PetscConvEst object
221: - r - The refinement level
223: Options database keys:
224: . -convest_monitor - Activate the monitor
226: Level: intermediate
228: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve()
229: @*/
230: PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
231: {
232: MPI_Comm comm;
233: PetscInt f;
237: if (ce->monitor) {
238: PetscInt *dofs = &ce->dofs[r*ce->Nf];
239: PetscReal *errors = &ce->errors[r*ce->Nf];
241: PetscObjectGetComm((PetscObject) ce, &comm);
242: PetscPrintf(comm, "N: ");
243: if (ce->Nf > 1) {PetscPrintf(comm, "[");}
244: for (f = 0; f < ce->Nf; ++f) {
245: if (f > 0) {PetscPrintf(comm, ", ");}
246: PetscPrintf(comm, "%7D", dofs[f]);
247: }
248: if (ce->Nf > 1) {PetscPrintf(comm, "]");}
249: PetscPrintf(comm, " ");
250: PetscPrintf(comm, "L_2 Error: ");
251: if (ce->Nf > 1) {PetscPrintf(comm, "[");}
252: for (f = 0; f < ce->Nf; ++f) {
253: if (f > 0) {PetscPrintf(comm, ", ");}
254: if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
255: else {PetscPrintf(comm, "%g", (double) errors[f]);}
256: }
257: if (ce->Nf > 1) {PetscPrintf(comm, "]");}
258: PetscPrintf(comm, "\n");
259: }
260: return(0);
261: }
263: static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
264: {
265: PetscClassId id;
269: PetscObjectGetClassId(ce->solver, &id);
270: if (id != SNES_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES");
271: SNESGetDM((SNES) ce->solver, &ce->idm);
272: return(0);
273: }
275: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
276: {
280: DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u);
281: return(0);
282: }
284: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
285: {
289: DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors);
290: return(0);
291: }
293: static PetscErrorCode PetscConvEstSetJacobianNullspace_Private(PetscConvEst ce, SNES snes)
294: {
295: DM dm;
296: PetscInt f;
300: SNESGetDM(snes, &dm);
301: for (f = 0; f < ce->Nf; ++f) {
302: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
304: DMGetNullSpaceConstructor(dm, f, &nspconstr);
305: if (nspconstr) {
306: MatNullSpace nullsp;
307: Mat J;
309: (*nspconstr)(dm, f, f,&nullsp);
310: SNESSetUp(snes);
311: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
312: MatSetNullSpace(J, nullsp);
313: MatNullSpaceDestroy(&nullsp);
314: break;
315: }
316: }
317: return(0);
318: }
320: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
321: {
322: SNES snes = (SNES) ce->solver;
323: DM *dm;
324: PetscObject disc;
325: PetscReal *x, *y, slope, intercept;
326: PetscInt Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
327: void *ctx;
331: if (ce->r != 2.0) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double) ce->r);
332: DMGetDimension(ce->idm, &dim);
333: DMGetApplicationContext(ce->idm, &ctx);
334: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
335: DMGetRefineLevel(ce->idm, &oldlevel);
336: PetscMalloc1((Nr+1), &dm);
337: /* Loop over meshes */
338: dm[0] = ce->idm;
339: for (r = 0; r <= Nr; ++r) {
340: Vec u;
341: #if defined(PETSC_USE_LOG)
342: PetscLogStage stage;
343: #endif
344: char stageName[PETSC_MAX_PATH_LEN];
345: const char *dmname, *uname;
347: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
348: #if defined(PETSC_USE_LOG)
349: PetscLogStageGetId(stageName, &stage);
350: if (stage < 0) {PetscLogStageRegister(stageName, &stage);}
351: #endif
352: PetscLogStagePush(stage);
353: if (r > 0) {
354: if (!ce->noRefine) {
355: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
356: DMSetCoarseDM(dm[r], dm[r-1]);
357: } else {
358: DM cdm, rcdm;
360: DMClone(dm[r-1], &dm[r]);
361: DMCopyDisc(dm[r-1], dm[r]);
362: DMGetCoordinateDM(dm[r-1], &cdm);
363: DMGetCoordinateDM(dm[r], &rcdm);
364: DMCopyDisc(cdm, rcdm);
365: }
366: DMCopyTransform(ce->idm, dm[r]);
367: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
368: PetscObjectSetName((PetscObject) dm[r], dmname);
369: for (f = 0; f < ce->Nf; ++f) {
370: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
372: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
373: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
374: }
375: }
376: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
377: /* Create solution */
378: DMCreateGlobalVector(dm[r], &u);
379: DMGetField(dm[r], 0, NULL, &disc);
380: PetscObjectGetName(disc, &uname);
381: PetscObjectSetName((PetscObject) u, uname);
382: /* Setup solver */
383: SNESReset(snes);
384: SNESSetDM(snes, dm[r]);
385: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
386: SNESSetFromOptions(snes);
387: /* Set nullspace for Jacobian */
388: PetscConvEstSetJacobianNullspace_Private(ce, snes);
389: /* Create initial guess */
390: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
391: SNESSolve(snes, NULL, u);
392: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
393: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*ce->Nf]);
394: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
395: for (f = 0; f < ce->Nf; ++f) {
396: PetscSection s, fs;
397: PetscInt lsize;
399: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
400: DMGetLocalSection(dm[r], &s);
401: PetscSectionGetField(s, f, &fs);
402: PetscSectionGetConstrainedStorageSize(fs, &lsize);
403: MPI_Allreduce(&lsize, &ce->dofs[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) snes));
404: PetscLogEventSetDof(ce->event, f, ce->dofs[r*ce->Nf+f]);
405: PetscLogEventSetError(ce->event, f, ce->errors[r*ce->Nf+f]);
406: }
407: /* Monitor */
408: PetscConvEstMonitorDefault(ce, r);
409: if (!r) {
410: /* PCReset() does not wipe out the level structure */
411: KSP ksp;
412: PC pc;
414: SNESGetKSP(snes, &ksp);
415: KSPGetPC(ksp, &pc);
416: PCMGGetLevels(pc, &oldnlev);
417: }
418: /* Cleanup */
419: VecDestroy(&u);
420: PetscLogStagePop();
421: }
422: for (r = 1; r <= Nr; ++r) {
423: DMDestroy(&dm[r]);
424: }
425: /* Fit convergence rate */
426: PetscMalloc2(Nr+1, &x, Nr+1, &y);
427: for (f = 0; f < ce->Nf; ++f) {
428: for (r = 0; r <= Nr; ++r) {
429: x[r] = PetscLog10Real(ce->dofs[r*ce->Nf+f]);
430: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
431: }
432: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
433: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
434: alpha[f] = -slope * dim;
435: }
436: PetscFree2(x, y);
437: PetscFree(dm);
438: /* Restore solver */
439: SNESReset(snes);
440: {
441: /* PCReset() does not wipe out the level structure */
442: KSP ksp;
443: PC pc;
445: SNESGetKSP(snes, &ksp);
446: KSPGetPC(ksp, &pc);
447: PCMGSetLevels(pc, oldnlev, NULL);
448: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
449: }
450: SNESSetDM(snes, ce->idm);
451: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
452: SNESSetFromOptions(snes);
453: PetscConvEstSetJacobianNullspace_Private(ce, snes);
454: return(0);
455: }
457: /*@
458: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
460: Not collective
462: Input Parameter:
463: . ce - The PetscConvEst object
465: Output Parameter:
466: . alpha - The convergence rate for each field
468: Note: The convergence rate alpha is defined by
469: $ || u_\Delta - u_exact || < C \Delta^alpha
470: where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the
471: spatial resolution and \Delta t for the temporal resolution.
473: We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
474: based upon the exact solution in the DS, and then fit the result to our model above using linear regression.
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
480: Level: intermediate
482: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve()
483: @*/
484: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
485: {
486: PetscInt f;
490: if (ce->event < 0) {PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);}
491: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
492: (*ce->ops->getconvrate)(ce, alpha);
493: return(0);
494: }
496: /*@
497: PetscConvEstRateView - Displays the convergence rate to a viewer
499: Collective on SNES
501: Parameter:
502: + snes - iterative context obtained from SNESCreate()
503: . alpha - the convergence rate for each field
504: - viewer - the viewer to display the reason
506: Options Database Keys:
507: . -snes_convergence_estimate - print the convergence rate
509: Level: developer
511: .seealso: PetscConvEstGetRate()
512: @*/
513: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
514: {
515: PetscBool isAscii;
519: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
520: if (isAscii) {
521: PetscInt Nf = ce->Nf, f;
523: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
524: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
525: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
526: for (f = 0; f < Nf; ++f) {
527: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
528: PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);
529: }
530: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
531: PetscViewerASCIIPrintf(viewer, "\n");
532: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
533: }
534: return(0);
535: }
537: /*@
538: PetscConvEstCreate - Create a PetscConvEst object
540: Collective
542: Input Parameter:
543: . comm - The communicator for the PetscConvEst object
545: Output Parameter:
546: . ce - The PetscConvEst object
548: Level: beginner
550: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
551: @*/
552: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
553: {
558: PetscSysInitializePackage();
559: PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
560: (*ce)->monitor = PETSC_FALSE;
561: (*ce)->r = 2.0;
562: (*ce)->Nr = 4;
563: (*ce)->event = -1;
564: (*ce)->ops->setsolver = PetscConvEstSetSNES_Private;
565: (*ce)->ops->initguess = PetscConvEstInitGuessSNES_Private;
566: (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
567: (*ce)->ops->getconvrate = PetscConvEstGetConvRateSNES_Private;
568: return(0);
569: }