Actual source code: convest.c
petsc-3.12.5 2020-03-29
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: }
15: /*@
16: PetscConvEstCreate - Create a PetscConvEst object
18: Collective
20: Input Parameter:
21: . comm - The communicator for the PetscConvEst object
23: Output Parameter:
24: . ce - The PetscConvEst object
26: Level: beginner
28: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
29: @*/
30: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
31: {
36: PetscSysInitializePackage();
37: PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
38: (*ce)->monitor = PETSC_FALSE;
39: (*ce)->Nr = 4;
40: return(0);
41: }
43: /*@
44: PetscConvEstDestroy - Destroys a PetscConvEst object
46: Collective on PetscConvEst
48: Input Parameter:
49: . ce - The PetscConvEst object
51: Level: beginner
53: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
54: @*/
55: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
56: {
60: if (!*ce) return(0);
62: if (--((PetscObject)(*ce))->refct > 0) {
63: *ce = NULL;
64: return(0);
65: }
66: PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs);
67: PetscFree((*ce)->errors);
68: PetscHeaderDestroy(ce);
69: return(0);
70: }
72: /*@
73: PetscConvEstSetFromOptions - Sets a PetscConvEst object from options
75: Collective on PetscConvEst
77: Input Parameters:
78: . ce - The PetscConvEst object
80: Level: beginner
82: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
83: @*/
84: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
85: {
89: PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");
90: PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
91: PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);
92: PetscOptionsEnd();
93: return(0);
94: }
96: /*@
97: PetscConvEstView - Views a PetscConvEst object
99: Collective on PetscConvEst
101: Input Parameters:
102: + ce - The PetscConvEst object
103: - viewer - The PetscViewer object
105: Level: beginner
107: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
108: @*/
109: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
110: {
114: PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
115: PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
116: return(0);
117: }
119: /*@
120: PetscConvEstGetSolver - Gets the solver used to produce discrete solutions
122: Not collective
124: Input Parameter:
125: . ce - The PetscConvEst object
127: Output Parameter:
128: . snes - The solver
130: Level: intermediate
132: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
133: @*/
134: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes)
135: {
139: *snes = ce->snes;
140: return(0);
141: }
143: /*@
144: PetscConvEstSetSolver - Sets the solver used to produce discrete solutions
146: Not collective
148: Input Parameters:
149: + ce - The PetscConvEst object
150: - snes - The solver
152: Level: intermediate
154: Note: The solver MUST have an attached DM/DS, so that we know the exact solution
156: .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
157: @*/
158: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes)
159: {
165: ce->snes = snes;
166: SNESGetDM(ce->snes, &ce->idm);
167: return(0);
168: }
170: /*@
171: PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence
173: Collective on PetscConvEst
175: Input Parameters:
176: . ce - The PetscConvEst object
178: Level: beginner
180: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
181: @*/
182: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
183: {
184: PetscDS prob;
185: PetscInt f;
189: DMGetDS(ce->idm, &prob);
190: PetscDSGetNumFields(prob, &ce->Nf);
191: PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
192: PetscMalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);
193: for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private;
194: for (f = 0; f < ce->Nf; ++f) {
195: PetscDSGetExactSolution(prob, f, &ce->exactSol[f], &ce->ctxs[f]);
196: 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);
197: }
198: return(0);
199: }
201: /*@
202: PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization
204: Not collective
206: Input Parameter:
207: . ce - The PetscConvEst object
209: Output Parameter:
210: . alpha - The convergence rate for each field
212: Note: The convergence rate alpha is defined by
213: $ || u_h - u_exact || < C h^alpha
214: where u_h is the discrete solution, and h is a measure of the discretization size.
216: We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS,
217: and then fit the result to our model above using linear regression.
219: Options database keys:
220: . -snes_convergence_estimate : Execute convergence estimation and print out the rate
222: Level: intermediate
224: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
225: @*/
226: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
227: {
228: DM *dm;
229: PetscObject disc;
230: MPI_Comm comm;
231: const char *uname, *dmname;
232: void *ctx;
233: Vec u;
234: PetscReal t = 0.0, *x, *y, slope, intercept;
235: PetscInt *dof, dim, Nr = ce->Nr, r, f, oldlevel, oldnlev;
236: PetscLogEvent event;
240: PetscObjectGetComm((PetscObject) ce, &comm);
241: DMGetDimension(ce->idm, &dim);
242: DMGetApplicationContext(ce->idm, &ctx);
243: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
244: DMGetRefineLevel(ce->idm, &oldlevel);
245: PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);
246: dm[0] = ce->idm;
247: for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
248: /* Loop over meshes */
249: PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);
250: for (r = 0; r <= Nr; ++r) {
251: PetscLogStage stage;
252: char stageName[PETSC_MAX_PATH_LEN];
254: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
255: PetscLogStageRegister(stageName, &stage);
256: PetscLogStagePush(stage);
257: if (r > 0) {
258: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
259: DMSetCoarseDM(dm[r], dm[r-1]);
260: DMCopyDisc(ce->idm, dm[r]);
261: DMCopyTransform(ce->idm, dm[r]);
262: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
263: PetscObjectSetName((PetscObject) dm[r], dmname);
264: for (f = 0; f <= ce->Nf; ++f) {
265: PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
266: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
267: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
268: }
269: }
270: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
271: /* Create solution */
272: DMCreateGlobalVector(dm[r], &u);
273: DMGetField(dm[r], 0, NULL, &disc);
274: PetscObjectGetName(disc, &uname);
275: PetscObjectSetName((PetscObject) u, uname);
276: /* Setup solver */
277: SNESReset(ce->snes);
278: SNESSetDM(ce->snes, dm[r]);
279: DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
280: SNESSetFromOptions(ce->snes);
281: /* Create initial guess */
282: DMProjectFunction(dm[r], t, ce->initGuess, ce->ctxs, INSERT_VALUES, u);
283: SNESSolve(ce->snes, NULL, u);
284: PetscLogEventBegin(event, ce, 0, 0, 0);
285: DMComputeL2FieldDiff(dm[r], t, ce->exactSol, ce->ctxs, u, &ce->errors[r*ce->Nf]);
286: PetscLogEventEnd(event, ce, 0, 0, 0);
287: for (f = 0; f < ce->Nf; ++f) {
288: PetscSection s, fs;
289: PetscInt lsize;
291: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
292: DMGetLocalSection(dm[r], &s);
293: PetscSectionGetField(s, f, &fs);
294: PetscSectionGetConstrainedStorageSize(fs, &lsize);
295: MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));
296: PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);
297: PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);
298: }
299: /* Monitor */
300: if (ce->monitor) {
301: PetscReal *errors = &ce->errors[r*ce->Nf];
303: PetscPrintf(comm, "L_2 Error: ");
304: if (ce->Nf > 1) {PetscPrintf(comm, "[");}
305: for (f = 0; f < ce->Nf; ++f) {
306: if (f > 0) {PetscPrintf(comm, ", ");}
307: if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
308: else {PetscPrintf(comm, "%g", (double)errors[f]);}
309: }
310: if (ce->Nf > 1) {PetscPrintf(comm, "]");}
311: PetscPrintf(comm, "\n");
312: }
313: if (!r) {
314: /* PCReset() does not wipe out the level structure */
315: KSP ksp;
316: PC pc;
318: SNESGetKSP(ce->snes, &ksp);
319: KSPGetPC(ksp, &pc);
320: PCMGGetLevels(pc, &oldnlev);
321: }
322: /* Cleanup */
323: VecDestroy(&u);
324: PetscLogStagePop();
325: }
326: for (r = 1; r <= Nr; ++r) {
327: DMDestroy(&dm[r]);
328: }
329: /* Fit convergence rate */
330: PetscMalloc2(Nr+1, &x, Nr+1, &y);
331: for (f = 0; f < ce->Nf; ++f) {
332: for (r = 0; r <= Nr; ++r) {
333: x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
334: y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
335: }
336: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
337: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
338: alpha[f] = -slope * dim;
339: }
340: PetscFree2(x, y);
341: PetscFree2(dm, dof);
342: /* Restore solver */
343: SNESReset(ce->snes);
344: {
345: /* PCReset() does not wipe out the level structure */
346: KSP ksp;
347: PC pc;
349: SNESGetKSP(ce->snes, &ksp);
350: KSPGetPC(ksp, &pc);
351: PCMGSetLevels(pc, oldnlev, NULL);
352: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
353: }
354: SNESSetDM(ce->snes, ce->idm);
355: DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
356: SNESSetFromOptions(ce->snes);
357: return(0);
358: }
360: /*@
361: PetscConvEstRateView - Displays the convergence rate to a viewer
363: Collective on SNES
365: Parameter:
366: + snes - iterative context obtained from SNESCreate()
367: . alpha - the convergence rate for each field
368: - viewer - the viewer to display the reason
370: Options Database Keys:
371: . -snes_convergence_estimate - print the convergence rate
373: Level: developer
375: .seealso: PetscConvEstGetRate()
376: @*/
377: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
378: {
379: PetscBool isAscii;
383: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
384: if (isAscii) {
385: DM dm;
386: PetscInt Nf, f;
388: SNESGetDM(ce->snes, &dm);
389: DMGetNumFields(dm, &Nf);
390: PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
391: PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
392: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
393: for (f = 0; f < Nf; ++f) {
394: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
395: PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);
396: }
397: if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
398: PetscViewerASCIIPrintf(viewer, "\n");
399: PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
400: }
401: return(0);
402: }