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: }