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

 19:   Input Parameter:
 20: . ce - The PetscConvEst object

 22:   Level: beginner

 24: .seealso: `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
 25: @*/
 26: PetscErrorCode PetscConvEstDestroy(PetscConvEst *ce)
 27: {
 28:   if (!*ce) return 0;
 30:   if (--((PetscObject)(*ce))->refct > 0) {
 31:     *ce = NULL;
 32:     return 0;
 33:   }
 34:   PetscFree3((*ce)->initGuess, (*ce)->exactSol, (*ce)->ctxs);
 35:   PetscFree2((*ce)->dofs, (*ce)->errors);
 36:   PetscHeaderDestroy(ce);
 37:   return 0;
 38: }

 40: /*@
 41:   PetscConvEstSetFromOptions - Sets a PetscConvEst object from options

 43:   Collective

 45:   Input Parameters:
 46: . ce - The PetscConvEst object

 48:   Level: beginner

 50: .seealso: `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
 51: @*/
 52: PetscErrorCode PetscConvEstSetFromOptions(PetscConvEst ce)
 53: {
 54:   PetscOptionsBegin(PetscObjectComm((PetscObject)ce), "", "Convergence Estimator Options", "PetscConvEst");
 55:   PetscOptionsInt("-convest_num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
 56:   PetscOptionsReal("-convest_refine_factor", "The increase in resolution in each dimension", "PetscConvEst", ce->r, &ce->r, NULL);
 57:   PetscOptionsBool("-convest_monitor", "Monitor the error for each convergence check", "PetscConvEst", ce->monitor, &ce->monitor, NULL);
 58:   PetscOptionsBool("-convest_no_refine", "Debugging flag to run on the same mesh each time", "PetscConvEst", ce->noRefine, &ce->noRefine, NULL);
 59:   PetscOptionsEnd();
 60:   return 0;
 61: }

 63: /*@
 64:   PetscConvEstView - Views a PetscConvEst object

 66:   Collective

 68:   Input Parameters:
 69: + ce     - The PetscConvEst object
 70: - viewer - The PetscViewer object

 72:   Level: beginner

 74: .seealso: `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
 75: @*/
 76: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
 77: {
 78:   PetscObjectPrintClassNamePrefixType((PetscObject)ce, viewer);
 79:   PetscViewerASCIIPrintf(viewer, "ConvEst with %" PetscInt_FMT " levels\n", ce->Nr + 1);
 80:   return 0;
 81: }

 83: /*@
 84:   PetscConvEstGetSolver - Gets the solver used to produce discrete solutions

 86:   Not collective

 88:   Input Parameter:
 89: . ce     - The PetscConvEst object

 91:   Output Parameter:
 92: . solver - The solver

 94:   Level: intermediate

 96: .seealso: `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
 97: @*/
 98: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
 99: {
102:   *solver = ce->solver;
103:   return 0;
104: }

106: /*@
107:   PetscConvEstSetSolver - Sets the solver used to produce discrete solutions

109:   Not collective

111:   Input Parameters:
112: + ce     - The PetscConvEst object
113: - solver - The solver

115:   Level: intermediate

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

119: .seealso: `PetscConvEstGetSNES()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
120: @*/
121: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
122: {
125:   ce->solver = solver;
126:   PetscUseTypeMethod(ce, setsolver, solver);
127:   return 0;
128: }

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

133:   Collective

135:   Input Parameters:
136: . ce - The PetscConvEst object

138:   Level: beginner

140: .seealso: `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`
141: @*/
142: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
143: {
144:   PetscInt Nf, f, Nds, s;

146:   DMGetNumFields(ce->idm, &Nf);
147:   ce->Nf = PetscMax(Nf, 1);
148:   PetscMalloc2((ce->Nr + 1) * ce->Nf, &ce->dofs, (ce->Nr + 1) * ce->Nf, &ce->errors);
149:   PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);
150:   for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private;
151:   DMGetNumDS(ce->idm, &Nds);
152:   for (s = 0; s < Nds; ++s) {
153:     PetscDS         ds;
154:     DMLabel         label;
155:     IS              fieldIS;
156:     const PetscInt *fields;
157:     PetscInt        dsNf;

159:     DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds);
160:     PetscDSGetNumFields(ds, &dsNf);
161:     if (fieldIS) ISGetIndices(fieldIS, &fields);
162:     for (f = 0; f < dsNf; ++f) {
163:       const PetscInt field = fields[f];
164:       PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]);
165:     }
166:     if (fieldIS) ISRestoreIndices(fieldIS, &fields);
167:   }
169:   return 0;
170: }

172: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
173: {
177:   PetscUseTypeMethod(ce, initguess, r, dm, u);
178:   return 0;
179: }

181: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
182: {
187:   PetscUseTypeMethod(ce, computeerror, r, dm, u, errors);
188:   return 0;
189: }

191: /*@
192:   PetscConvEstMonitorDefault - Monitors the convergence estimation loop

194:   Collective

196:   Input Parameters:
197: + ce - The `PetscConvEst` object
198: - r  - The refinement level

200:   Options database keys:
201: . -convest_monitor - Activate the monitor

203:   Level: intermediate

205: .seealso: `PetscConvEst`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`, `SNESSolve()`, `TSSolve()`
206: @*/
207: PetscErrorCode PetscConvEstMonitorDefault(PetscConvEst ce, PetscInt r)
208: {
209:   MPI_Comm comm;
210:   PetscInt f;

212:   if (ce->monitor) {
213:     PetscInt  *dofs   = &ce->dofs[r * ce->Nf];
214:     PetscReal *errors = &ce->errors[r * ce->Nf];

216:     PetscObjectGetComm((PetscObject)ce, &comm);
217:     PetscPrintf(comm, "N: ");
218:     if (ce->Nf > 1) PetscPrintf(comm, "[");
219:     for (f = 0; f < ce->Nf; ++f) {
220:       if (f > 0) PetscPrintf(comm, ", ");
221:       PetscPrintf(comm, "%7" PetscInt_FMT, dofs[f]);
222:     }
223:     if (ce->Nf > 1) PetscPrintf(comm, "]");
224:     PetscPrintf(comm, "  ");
225:     PetscPrintf(comm, "L_2 Error: ");
226:     if (ce->Nf > 1) PetscPrintf(comm, "[");
227:     for (f = 0; f < ce->Nf; ++f) {
228:       if (f > 0) PetscPrintf(comm, ", ");
229:       if (errors[f] < 1.0e-11) PetscPrintf(comm, "< 1e-11");
230:       else PetscPrintf(comm, "%g", (double)errors[f]);
231:     }
232:     if (ce->Nf > 1) PetscPrintf(comm, "]");
233:     PetscPrintf(comm, "\n");
234:   }
235:   return 0;
236: }

238: static PetscErrorCode PetscConvEstSetSNES_Private(PetscConvEst ce, PetscObject solver)
239: {
240:   PetscClassId id;

242:   PetscObjectGetClassId(ce->solver, &id);
244:   SNESGetDM((SNES)ce->solver, &ce->idm);
245:   return 0;
246: }

248: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
249: {
250:   DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u);
251:   return 0;
252: }

254: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
255: {
256:   DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors);
257:   return 0;
258: }

260: static PetscErrorCode PetscConvEstSetJacobianNullSpace_Private(PetscConvEst ce, SNES snes)
261: {
262:   DM       dm;
263:   PetscInt f;

265:   SNESGetDM(snes, &dm);
266:   for (f = 0; f < ce->Nf; ++f) {
267:     PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);

269:     DMGetNullSpaceConstructor(dm, f, &nspconstr);
270:     if (nspconstr) {
271:       MatNullSpace nullsp;
272:       Mat          J;

274:       (*nspconstr)(dm, f, f, &nullsp);
275:       SNESSetUp(snes);
276:       SNESGetJacobian(snes, &J, NULL, NULL, NULL);
277:       MatSetNullSpace(J, nullsp);
278:       MatNullSpaceDestroy(&nullsp);
279:       break;
280:     }
281:   }
282:   return 0;
283: }

285: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
286: {
287:   SNES        snes = (SNES)ce->solver;
288:   DM         *dm;
289:   PetscObject disc;
290:   PetscReal  *x, *y, slope, intercept;
291:   PetscInt    Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
292:   void       *ctx;

295:   DMGetDimension(ce->idm, &dim);
296:   DMGetApplicationContext(ce->idm, &ctx);
297:   DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
298:   DMGetRefineLevel(ce->idm, &oldlevel);
299:   PetscMalloc1((Nr + 1), &dm);
300:   /* Loop over meshes */
301:   dm[0] = ce->idm;
302:   for (r = 0; r <= Nr; ++r) {
303:     Vec u;
304: #if defined(PETSC_USE_LOG)
305:     PetscLogStage stage;
306: #endif
307:     char        stageName[PETSC_MAX_PATH_LEN];
308:     const char *dmname, *uname;

310:     PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r);
311: #if defined(PETSC_USE_LOG)
312:     PetscLogStageGetId(stageName, &stage);
313:     if (stage < 0) PetscLogStageRegister(stageName, &stage);
314: #endif
315:     PetscLogStagePush(stage);
316:     if (r > 0) {
317:       if (!ce->noRefine) {
318:         DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]);
319:         DMSetCoarseDM(dm[r], dm[r - 1]);
320:       } else {
321:         DM cdm, rcdm;

323:         DMClone(dm[r - 1], &dm[r]);
324:         DMCopyDisc(dm[r - 1], dm[r]);
325:         DMGetCoordinateDM(dm[r - 1], &cdm);
326:         DMGetCoordinateDM(dm[r], &rcdm);
327:         DMCopyDisc(cdm, rcdm);
328:       }
329:       DMCopyTransform(ce->idm, dm[r]);
330:       PetscObjectGetName((PetscObject)dm[r - 1], &dmname);
331:       PetscObjectSetName((PetscObject)dm[r], dmname);
332:       for (f = 0; f < ce->Nf; ++f) {
333:         PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);

335:         DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr);
336:         DMSetNullSpaceConstructor(dm[r], f, nspconstr);
337:       }
338:     }
339:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
340:     /* Create solution */
341:     DMCreateGlobalVector(dm[r], &u);
342:     DMGetField(dm[r], 0, NULL, &disc);
343:     PetscObjectGetName(disc, &uname);
344:     PetscObjectSetName((PetscObject)u, uname);
345:     /* Setup solver */
346:     SNESReset(snes);
347:     SNESSetDM(snes, dm[r]);
348:     DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
349:     SNESSetFromOptions(snes);
350:     /* Set nullspace for Jacobian */
351:     PetscConvEstSetJacobianNullSpace_Private(ce, snes);
352:     /* Create initial guess */
353:     PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
354:     SNESSolve(snes, NULL, u);
355:     PetscLogEventBegin(ce->event, ce, 0, 0, 0);
356:     PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * ce->Nf]);
357:     PetscLogEventEnd(ce->event, ce, 0, 0, 0);
358:     for (f = 0; f < ce->Nf; ++f) {
359:       PetscSection s, fs;
360:       PetscInt     lsize;

362:       /* Could use DMGetOutputDM() to add in Dirichlet dofs */
363:       DMGetLocalSection(dm[r], &s);
364:       PetscSectionGetField(s, f, &fs);
365:       PetscSectionGetConstrainedStorageSize(fs, &lsize);
366:       MPI_Allreduce(&lsize, &ce->dofs[r * ce->Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes));
367:       PetscLogEventSetDof(ce->event, f, ce->dofs[r * ce->Nf + f]);
368:       PetscLogEventSetError(ce->event, f, ce->errors[r * ce->Nf + f]);
369:     }
370:     /* Monitor */
371:     PetscConvEstMonitorDefault(ce, r);
372:     if (!r) {
373:       /* PCReset() does not wipe out the level structure */
374:       KSP ksp;
375:       PC  pc;

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

406:     SNESGetKSP(snes, &ksp);
407:     KSPGetPC(ksp, &pc);
408:     PCMGSetLevels(pc, oldnlev, NULL);
409:     DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
410:   }
411:   SNESSetDM(snes, ce->idm);
412:   DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
413:   SNESSetFromOptions(snes);
414:   PetscConvEstSetJacobianNullSpace_Private(ce, snes);
415:   return 0;
416: }

418: /*@
419:   PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization

421:   Not collective

423:   Input Parameter:
424: . ce   - The `PetscConvEst` object

426:   Output Parameter:
427: . alpha - The convergence rate for each field

429:   Options database keys:
430: + -snes_convergence_estimate - Execute convergence estimation inside `SNESSolve()` and print out the rate
431: - -ts_convergence_estimate - Execute convergence estimation inside `TSSolve()` and print out the rate

433:   Notes:
434:   The convergence rate alpha is defined by
435: $ || u_\Delta - u_exact || < C \Delta^alpha
436:   where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the
437:   spatial resolution and \Delta t for the temporal resolution.

439:   We solve a series of problems using increasing resolution (refined meshes or decreased timesteps), calculate an error
440:   based upon the exact solution in the DS, and then fit the result to our model above using linear regression.

442:   Level: intermediate

444: .seealso: `PetscConvEstSetSolver()`, `PetscConvEstCreate()`, `PetscConvEstGetConvRate()`, `SNESSolve()`, `TSSolve()`
445: @*/
446: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
447: {
448:   PetscInt f;

450:   if (ce->event < 0) PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);
451:   for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
452:   PetscUseTypeMethod(ce, getconvrate, alpha);
453:   return 0;
454: }

456: /*@
457:   PetscConvEstRateView - Displays the convergence rate to a viewer

459:    Collective

461:    Parameter:
462: +  snes - iterative context obtained from `SNESCreate()`
463: .  alpha - the convergence rate for each field
464: -  viewer - the viewer to display the reason

466:    Options Database Keys:
467: .  -snes_convergence_estimate - print the convergence rate

469:    Level: developer

471: .seealso: `PetscConvEst`, `PetscConvEstGetRate()`
472: @*/
473: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, const PetscReal alpha[], PetscViewer viewer)
474: {
475:   PetscBool isAscii;

477:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii);
478:   if (isAscii) {
479:     PetscInt Nf = ce->Nf, f;

481:     PetscViewerASCIIAddTab(viewer, ((PetscObject)ce)->tablevel);
482:     PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
483:     if (Nf > 1) PetscViewerASCIIPrintf(viewer, "[");
484:     for (f = 0; f < Nf; ++f) {
485:       if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
486:       PetscViewerASCIIPrintf(viewer, "%#.2g", (double)alpha[f]);
487:     }
488:     if (Nf > 1) PetscViewerASCIIPrintf(viewer, "]");
489:     PetscViewerASCIIPrintf(viewer, "\n");
490:     PetscViewerASCIISubtractTab(viewer, ((PetscObject)ce)->tablevel);
491:   }
492:   return 0;
493: }

495: /*@
496:   PetscConvEstCreate - Create a `PetscConvEst` object

498:   Collective

500:   Input Parameter:
501: . comm - The communicator for the `PetscConvEst` object

503:   Output Parameter:
504: . ce   - The `PetscConvEst` object

506:   Level: beginner

508: .seealso: `PetscConvEst`, `PetscConvEstDestroy()`, `PetscConvEstGetConvRate()`
509: @*/
510: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
511: {
513:   PetscSysInitializePackage();
514:   PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
515:   (*ce)->monitor           = PETSC_FALSE;
516:   (*ce)->r                 = 2.0;
517:   (*ce)->Nr                = 4;
518:   (*ce)->event             = -1;
519:   (*ce)->ops->setsolver    = PetscConvEstSetSNES_Private;
520:   (*ce)->ops->initguess    = PetscConvEstInitGuessSNES_Private;
521:   (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
522:   (*ce)->ops->getconvrate  = PetscConvEstGetConvRateSNES_Private;
523:   return 0;
524: }