Actual source code: convest.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1:  #include <petscconvest.h>
  2:  #include <petscdmplex.h>
  3:  #include <petscds.h>
  4:  #include <petscblaslapack.h>

  6:  #include <petsc/private/petscconvestimpl.h>

  8: static PetscErrorCode zero_private(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
  9: {
 10:   PetscInt c;
 11:   for (c = 0; c < Nc; ++c) u[c] = 0.0;
 12:   return 0;
 13: }


 16: /*@
 17:   PetscConvEstCreate - Create a PetscConvEst object

 19:   Collective on MPI_Comm

 21:   Input Parameter:
 22: . comm - The communicator for the PetscConvEst object

 24:   Output Parameter:
 25: . ce   - The PetscConvEst object

 27:   Level: beginner

 29: .keywords: PetscConvEst, convergence, create
 30: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
 31: @*/
 32: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
 33: {

 38:   PetscSysInitializePackage();
 39:   PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
 40:   (*ce)->monitor = PETSC_FALSE;
 41:   (*ce)->Nr      = 4;
 42:   return(0);
 43: }

 45: /*@
 46:   PetscConvEstDestroy - Destroys a PetscConvEst object

 48:   Collective on PetscConvEst

 50:   Input Parameter:
 51: . ce - The PetscConvEst object

 53:   Level: beginner

 55: .keywords: PetscConvEst, convergence, destroy
 56: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
 57: @*/
 58: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
 59: {

 63:   if (!*ce) return(0);
 65:   if (--((PetscObject)(*ce))->refct > 0) {
 66:     *ce = NULL;
 67:     return(0);
 68:   }
 69:   PetscFree2((*ce)->initGuess, (*ce)->exactSol);
 70:   PetscFree((*ce)->errors);
 71:   PetscHeaderDestroy(ce);
 72:   return(0);
 73: }

 75: /*@
 76:   PetscConvEstSetFromOptions - Sets a PetscConvEst object from options

 78:   Collective on PetscConvEst

 80:   Input Parameters:
 81: . ce - The PetscConvEst object

 83:   Level: beginner

 85: .keywords: PetscConvEst, convergence, options
 86: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
 87: @*/
 88: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
 89: {

 93:   PetscOptionsBegin(PetscObjectComm((PetscObject) ce), "", "Convergence Estimator Options", "PetscConvEst");
 94:   PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
 95:   PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);
 96:   PetscOptionsEnd();
 97:   return(0);
 98: }

100: /*@
101:   PetscConvEstView - Views a PetscConvEst object

103:   Collective on PetscConvEst

105:   Input Parameters:
106: + ce     - The PetscConvEst object
107: - viewer - The PetscViewer object

109:   Level: beginner

111: .keywords: PetscConvEst, convergence, view
112: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
113: @*/
114: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
115: {

119:   PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
120:   PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
121:   return(0);
122: }

124: /*@
125:   PetscConvEstGetSolver - Gets the solver used to produce discrete solutions

127:   Not collective

129:   Input Parameter:
130: . ce   - The PetscConvEst object

132:   Output Parameter:
133: . snes - The solver

135:   Level: intermediate

137: .keywords: PetscConvEst, convergence
138: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
139: @*/
140: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes)
141: {
145:   *snes = ce->snes;
146:   return(0);
147: }

149: /*@
150:   PetscConvEstSetSolver - Sets the solver used to produce discrete solutions

152:   Not collective

154:   Input Parameters:
155: + ce   - The PetscConvEst object
156: - snes - The solver

158:   Level: intermediate

160:   Note: The solver MUST have an attached DM/DS, so that we know the exact solution

162: .keywords: PetscConvEst, convergence
163: .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
164: @*/
165: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes)
166: {

172:   ce->snes = snes;
173:   SNESGetDM(ce->snes, &ce->idm);
174:   return(0);
175: }

177: /*@
178:   PetscConvEstSetUp - After the solver is specified, we create structures for estimating convergence

180:   Collective on PetscConvEst

182:   Input Parameters:
183: . ce - The PetscConvEst object

185:   Level: beginner

187: .keywords: PetscConvEst, convergence, setup
188: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
189: @*/
190: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
191: {
192:   PetscDS        prob;
193:   PetscInt       f;

197:   DMGetDS(ce->idm, &prob);
198:   PetscDSGetNumFields(prob, &ce->Nf);
199:   PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
200:   PetscMalloc2(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol);
201:   for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private;
202:   for (f = 0; f < ce->Nf; ++f) {
203:     PetscDSGetExactSolution(prob, f, &ce->exactSol[f]);
204:     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);
205:   }
206:   return(0);
207: }

209: static PetscErrorCode PetscConvEstLinearRegression_Private(PetscConvEst ce, PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept)
210: {
211:   PetscScalar    H[4];
212:   PetscReal     *X, *Y, beta[2];
213:   PetscInt       i, j, k;

217:   *slope = *intercept = 0.0;
218:   PetscMalloc2(n*2, &X, n*2, &Y);
219:   for (k = 0; k < n; ++k) {
220:     /* X[n,2] = [1, x] */
221:     X[k*2+0] = 1.0;
222:     X[k*2+1] = x[k];
223:   }
224:   /* H = X^T X */
225:   for (i = 0; i < 2; ++i) {
226:     for (j = 0; j < 2; ++j) {
227:       H[i*2+j] = 0.0;
228:       for (k = 0; k < n; ++k) {
229:         H[i*2+j] += X[k*2+i] * X[k*2+j];
230:       }
231:     }
232:   }
233:   /* H = (X^T X)^{-1} */
234:   {
235:     PetscBLASInt two = 2, ipiv[2], info;
236:     PetscScalar  work[2];

238:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
239:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
240:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
241:     PetscFPTrapPop();
242:   }
243:     /* Y = H X^T */
244:   for (i = 0; i < 2; ++i) {
245:     for (k = 0; k < n; ++k) {
246:       Y[i*n+k] = 0.0;
247:       for (j = 0; j < 2; ++j) {
248:         Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j];
249:       }
250:     }
251:   }
252:   /* beta = Y error = [y-intercept, slope] */
253:   for (i = 0; i < 2; ++i) {
254:     beta[i] = 0.0;
255:     for (k = 0; k < n; ++k) {
256:       beta[i] += Y[i*n+k] * y[k];
257:     }
258:   }
259:   PetscFree2(X, Y);
260:   *intercept = beta[0];
261:   *slope     = beta[1];
262:   return(0);
263: }

265: /*@
266:   PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization

268:   Not collective

270:   Input Parameter:
271: . ce   - The PetscConvEst object

273:   Output Parameter:
274: . alpha - The convergence rate for each field

276:   Note: The convergence rate alpha is defined by
277: $ || u_h - u_exact || < C h^alpha
278: where u_h is the discrete solution, and h is a measure of the discretization size.

280: We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS,
281: and then fit the result to our model above using linear regression.

283:   Options database keys:
284: . -snes_convergence_estimate : Execute convergence estimation and print out the rate

286:   Level: intermediate

288: .keywords: PetscConvEst, convergence
289: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
290: @*/
291: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
292: {
293:   DM            *dm;
294:   PetscObject    disc;
295:   MPI_Comm       comm;
296:   const char    *uname, *dmname;
297:   void          *ctx;
298:   Vec            u;
299:   PetscReal      t = 0.0, *x, *y, slope, intercept;
300:   PetscInt      *dof, dim, Nr = ce->Nr, r, f, oldlevel, oldnlev;
301:   PetscLogEvent  event;

305:   PetscObjectGetComm((PetscObject) ce, &comm);
306:   DMGetDimension(ce->idm, &dim);
307:   DMGetApplicationContext(ce->idm, &ctx);
308:   DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
309:   DMGetRefineLevel(ce->idm, &oldlevel);
310:   PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);
311:   dm[0]  = ce->idm;
312:   for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
313:   /* Loop over meshes */
314:   PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &event);
315:   for (r = 0; r <= Nr; ++r) {
316:     PetscLogStage stage;
317:     char          stageName[PETSC_MAX_PATH_LEN];

319:     PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
320:     PetscLogStageRegister(stageName, &stage);
321:     PetscLogStagePush(stage);
322:     if (r > 0) {
323:       DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
324:       DMSetCoarseDM(dm[r], dm[r-1]);
325:       DMCopyDisc(ce->idm, dm[r]);
326:       PetscObjectGetName((PetscObject) dm[r-1], &dmname);
327:       PetscObjectSetName((PetscObject) dm[r], dmname);
328:       for (f = 0; f <= ce->Nf; ++f) {
329:         PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
330:         DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
331:         DMSetNullSpaceConstructor(dm[r],   f,  nspconstr);
332:       }
333:     }
334:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
335:     /* Create solution */
336:     DMCreateGlobalVector(dm[r], &u);
337:     DMGetField(dm[r], 0, NULL, &disc);
338:     PetscObjectGetName(disc, &uname);
339:     PetscObjectSetName((PetscObject) u, uname);
340:     /* Setup solver */
341:     SNESReset(ce->snes);
342:     SNESSetDM(ce->snes, dm[r]);
343:     DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
344:     SNESSetFromOptions(ce->snes);
345:     /* Create initial guess */
346:     DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
347:     SNESSolve(ce->snes, NULL, u);
348:     PetscLogEventBegin(event, ce, 0, 0, 0);
349:     DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);
350:     PetscLogEventEnd(event, ce, 0, 0, 0);
351:     for (f = 0; f < ce->Nf; ++f) {
352:       PetscSection s, fs;
353:       PetscInt     lsize;

355:       /* Could use DMGetOutputDM() to add in Dirichlet dofs */
356:       DMGetSection(dm[r], &s);
357:       PetscSectionGetField(s, f, &fs);
358:       PetscSectionGetConstrainedStorageSize(fs, &lsize);
359:       MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));
360:       PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);
361:       PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);
362:     }
363:     /* Monitor */
364:     if (ce->monitor) {
365:       PetscReal *errors = &ce->errors[r*ce->Nf];

367:       PetscPrintf(comm, "L_2 Error: ");
368:       if (ce->Nf > 1) {PetscPrintf(comm, "[");}
369:       for (f = 0; f < ce->Nf; ++f) {
370:         if (f > 0) {PetscPrintf(comm, ", ");}
371:         if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
372:         else                     {PetscPrintf(comm, "%g", (double)errors[f]);}
373:       }
374:       if (ce->Nf > 1) {PetscPrintf(comm, "]");}
375:       PetscPrintf(comm, "\n");
376:     }
377:     if (!r) {
378:       /* PCReset() does not wipe out the level structure */
379:       KSP ksp;
380:       PC  pc;

382:       SNESGetKSP(ce->snes, &ksp);
383:       KSPGetPC(ksp, &pc);
384:       PCMGGetLevels(pc, &oldnlev);
385:     }
386:     /* Cleanup */
387:     VecDestroy(&u);
388:     PetscLogStagePop();
389:   }
390:   for (r = 1; r <= Nr; ++r) {
391:     DMDestroy(&dm[r]);
392:   }
393:   /* Fit convergence rate */
394:   PetscMalloc2(Nr+1, &x, Nr+1, &y);
395:   for (f = 0; f < ce->Nf; ++f) {
396:     for (r = 0; r <= Nr; ++r) {
397:       x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
398:       y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
399:     }
400:     PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
401:     /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
402:     alpha[f] = -slope * dim;
403:   }
404:   PetscFree2(x, y);
405:   PetscFree2(dm, dof);
406:   /* Restore solver */
407:   SNESReset(ce->snes);
408:   {
409:     /* PCReset() does not wipe out the level structure */
410:     KSP ksp;
411:     PC  pc;

413:     SNESGetKSP(ce->snes, &ksp);
414:     KSPGetPC(ksp, &pc);
415:     PCMGSetLevels(pc, oldnlev, NULL);
416:     DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
417:   }
418:   SNESSetDM(ce->snes, ce->idm);
419:   DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
420:   SNESSetFromOptions(ce->snes);
421:   return(0);
422: }

424: /*@
425:   PetscConvEstRateView - Displays the convergence rate to a viewer

427:    Collective on SNES

429:    Parameter:
430: +  snes - iterative context obtained from SNESCreate()
431: .  alpha - the convergence rate for each field
432: -  viewer - the viewer to display the reason

434:    Options Database Keys:
435: .  -snes_convergence_estimate - print the convergence rate

437:    Level: developer

439: .seealso: PetscConvEstGetRate()
440: @*/
441: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
442: {
443:   PetscBool      isAscii;

447:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
448:   if (isAscii) {
449:     DM       dm;
450:     PetscInt Nf, f;

452:     SNESGetDM(ce->snes, &dm);
453:     DMGetNumFields(dm, &Nf);
454:     PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
455:     PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
456:     if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
457:     for (f = 0; f < Nf; ++f) {
458:       if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
459:       PetscViewerASCIIPrintf(viewer, "%#.2g", (double) alpha[f]);
460:     }
461:     if (Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
462:     PetscViewerASCIIPrintf(viewer, "\n");
463:     PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
464:   }
465:   return(0);
466: }