Actual source code: convest.c

petsc-master 2020-06-03
Report Typos and Errors
  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:   PetscFree((*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:   PetscOptionsEnd();
 65:   return(0);
 66: }

 68: /*@
 69:   PetscConvEstView - Views a PetscConvEst object

 71:   Collective on PetscConvEst

 73:   Input Parameters:
 74: + ce     - The PetscConvEst object
 75: - viewer - The PetscViewer object

 77:   Level: beginner

 79: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
 80: @*/
 81: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
 82: {

 86:   PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
 87:   PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
 88:   return(0);
 89: }

 91: /*@
 92:   PetscConvEstGetSolver - Gets the solver used to produce discrete solutions

 94:   Not collective

 96:   Input Parameter:
 97: . ce     - The PetscConvEst object

 99:   Output Parameter:
100: . solver - The solver

102:   Level: intermediate

104: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
105: @*/
106: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, PetscObject *solver)
107: {
111:   *solver = ce->solver;
112:   return(0);
113: }

115: /*@
116:   PetscConvEstSetSolver - Sets the solver used to produce discrete solutions

118:   Not collective

120:   Input Parameters:
121: + ce     - The PetscConvEst object
122: - solver - The solver

124:   Level: intermediate

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

128: .seealso: PetscConvEstGetSNES(), PetscConvEstCreate(), PetscConvEstGetConvRate()
129: @*/
130: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, PetscObject solver)
131: {

137:   ce->solver = solver;
138:   (*ce->ops->setsolver)(ce, solver);
139:   return(0);
140: }

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

145:   Collective on PetscConvEst

147:   Input Parameters:
148: . ce - The PetscConvEst object

150:   Level: beginner

152: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
153: @*/
154: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
155: {
156:   PetscInt       Nf, f, Nds, s;

160:   DMGetNumFields(ce->idm, &Nf);
161:   ce->Nf = PetscMax(Nf, 1);
162:   PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
163:   PetscCalloc3(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol, ce->Nf, &ce->ctxs);
164:   for (f = 0; f < Nf; ++f) ce->initGuess[f] = zero_private;
165:   DMGetNumDS(ce->idm, &Nds);
166:   for (s = 0; s < Nds; ++s) {
167:     PetscDS         ds;
168:     DMLabel         label;
169:     IS              fieldIS;
170:     const PetscInt *fields;
171:     PetscInt        dsNf;

173:     DMGetRegionNumDS(ce->idm, s, &label, &fieldIS, &ds);
174:     PetscDSGetNumFields(ds, &dsNf);
175:     if (fieldIS) {ISGetIndices(fieldIS, &fields);}
176:     for (f = 0; f < dsNf; ++f) {
177:       const PetscInt field = fields[f];
178:       PetscDSGetExactSolution(ds, field, &ce->exactSol[field], &ce->ctxs[field]);
179:     }
180:     if (fieldIS) {ISRestoreIndices(fieldIS, &fields);}
181:   }
182:   for (f = 0; f < Nf; ++f) {
183:     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);
184:   }
185:   return(0);
186: }

188: PetscErrorCode PetscConvEstComputeInitialGuess(PetscConvEst ce, PetscInt r, DM dm, Vec u)
189: {

196:   (*ce->ops->initguess)(ce, r, dm, u);
197:   return(0);
198: }

200: PetscErrorCode PetscConvEstComputeError(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
201: {

209:   (*ce->ops->computeerror)(ce, r, dm, u, errors);
210:   return(0);
211: }

213: static PetscErrorCode PetscConvEstMonitor_Private(PetscConvEst ce, PetscInt r)
214: {
215:   MPI_Comm       comm;
216:   PetscInt       f;

220:   if (ce->monitor) {
221:     PetscReal *errors = &ce->errors[r*ce->Nf];

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

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

243:   PetscObjectGetClassId(ce->solver, &id);
244:   if (id != SNES_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a SNES");
245:   SNESGetDM((SNES) ce->solver, &ce->idm);
246:   return(0);
247: }

249: static PetscErrorCode PetscConvEstInitGuessSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
250: {

254:   DMProjectFunction(dm, 0.0, ce->initGuess, ce->ctxs, INSERT_VALUES, u);
255:   return(0);
256: }

258: static PetscErrorCode PetscConvEstComputeErrorSNES_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
259: {

263:   DMComputeL2FieldDiff(dm, 0.0, ce->exactSol, ce->ctxs, u, errors);
264:   return(0);
265: }

267: static PetscErrorCode PetscConvEstGetConvRateSNES_Private(PetscConvEst ce, PetscReal alpha[])
268: {
269:   SNES           snes = (SNES) ce->solver;
270:   DM            *dm;
271:   PetscObject    disc;
272:   PetscReal     *x, *y, slope, intercept;
273:   PetscInt      *dof, Nr = ce->Nr, r, f, dim, oldlevel, oldnlev;
274:   void          *ctx;

278:   if (ce->r != 2.0) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double) ce->r);
279:   DMGetDimension(ce->idm, &dim);
280:   DMGetApplicationContext(ce->idm, &ctx);
281:   DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
282:   DMGetRefineLevel(ce->idm, &oldlevel);
283:   PetscMalloc2((Nr+1), &dm, (Nr+1)*ce->Nf, &dof);
284:   /* Loop over meshes */
285:   dm[0] = ce->idm;
286:   for (r = 0; r <= Nr; ++r) {
287:     Vec           u;
288: #if defined(PETSC_USE_LOG)
289:     PetscLogStage stage;
290: #endif
291:     char          stageName[PETSC_MAX_PATH_LEN];
292:     const char   *dmname, *uname;

294:     PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
295:     PetscLogStageRegister(stageName, &stage);
296:     PetscLogStagePush(stage);
297:     if (r > 0) {
298:       DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
299:       DMSetCoarseDM(dm[r], dm[r-1]);
300:       DMCopyTransform(ce->idm, dm[r]);
301:       PetscObjectGetName((PetscObject) dm[r-1], &dmname);
302:       PetscObjectSetName((PetscObject) dm[r], dmname);
303:       for (f = 0; f <= ce->Nf; ++f) {
304:         PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);

306:         DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
307:         DMSetNullSpaceConstructor(dm[r],   f,  nspconstr);
308:       }
309:     }
310:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
311:     /* Create solution */
312:     DMCreateGlobalVector(dm[r], &u);
313:     DMGetField(dm[r], 0, NULL, &disc);
314:     PetscObjectGetName(disc, &uname);
315:     PetscObjectSetName((PetscObject) u, uname);
316:     /* Setup solver */
317:     SNESReset(snes);
318:     SNESSetDM(snes, dm[r]);
319:     DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
320:     SNESSetFromOptions(snes);
321:     /* Create initial guess */
322:     PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
323:     SNESSolve(snes, NULL, u);
324:     PetscLogEventBegin(ce->event, ce, 0, 0, 0);
325:     PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*ce->Nf]);
326:     PetscLogEventEnd(ce->event, ce, 0, 0, 0);
327:     for (f = 0; f < ce->Nf; ++f) {
328:       PetscSection s, fs;
329:       PetscInt     lsize;

331:       /* Could use DMGetOutputDM() to add in Dirichlet dofs */
332:       DMGetLocalSection(dm[r], &s);
333:       PetscSectionGetField(s, f, &fs);
334:       PetscSectionGetConstrainedStorageSize(fs, &lsize);
335:       MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) snes));
336:       PetscLogEventSetDof(ce->event, f, dof[r*ce->Nf+f]);
337:       PetscLogEventSetError(ce->event, f, ce->errors[r*ce->Nf+f]);
338:     }
339:     /* Monitor */
340:     PetscConvEstMonitor_Private(ce, r);
341:     if (!r) {
342:       /* PCReset() does not wipe out the level structure */
343:       KSP ksp;
344:       PC  pc;

346:       SNESGetKSP(snes, &ksp);
347:       KSPGetPC(ksp, &pc);
348:       PCMGGetLevels(pc, &oldnlev);
349:     }
350:     /* Cleanup */
351:     VecDestroy(&u);
352:     PetscLogStagePop();
353:   }
354:   for (r = 1; r <= Nr; ++r) {
355:     DMDestroy(&dm[r]);
356:   }
357:   /* Fit convergence rate */
358:   PetscMalloc2(Nr+1, &x, Nr+1, &y);
359:   for (f = 0; f < ce->Nf; ++f) {
360:     for (r = 0; r <= Nr; ++r) {
361:       x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
362:       y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
363:     }
364:     PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
365:     /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
366:     alpha[f] = -slope * dim;
367:   }
368:   PetscFree2(x, y);
369:   PetscFree2(dm, dof);
370:   /* Restore solver */
371:   SNESReset(snes);
372:   {
373:     /* PCReset() does not wipe out the level structure */
374:     KSP ksp;
375:     PC  pc;

377:     SNESGetKSP(snes, &ksp);
378:     KSPGetPC(ksp, &pc);
379:     PCMGSetLevels(pc, oldnlev, NULL);
380:     DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
381:   }
382:   SNESSetDM(snes, ce->idm);
383:   DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
384:   SNESSetFromOptions(snes);
385:   return(0);
386: }

388: /*@
389:   PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization

391:   Not collective

393:   Input Parameter:
394: . ce   - The PetscConvEst object

396:   Output Parameter:
397: . alpha - The convergence rate for each field

399:   Note: The convergence rate alpha is defined by
400: $ || u_\Delta - u_exact || < C \Delta^alpha
401: where u_\Delta is the discrete solution, and Delta is a measure of the discretization size. We usually use h for the
402: spatial resolution and \Delta t for the temporal resolution.

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

407:   Options database keys:
408: + -snes_convergence_estimate : Execute convergence estimation inside SNESSolve() and print out the rate
409: - -ts_convergence_estimate : Execute convergence estimation inside TSSolve() and print out the rate

411:   Level: intermediate

413: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate(), SNESSolve(), TSSolve()
414: @*/
415: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal alpha[])
416: {
417:   PetscInt       f;

421:   if (ce->event < 0) {PetscLogEventRegister("ConvEst Error", PETSC_OBJECT_CLASSID, &ce->event);}
422:   for (f = 0; f < ce->Nf; ++f) alpha[f] = 0.0;
423:   (*ce->ops->getconvrate)(ce, alpha);
424:   return(0);
425: }

427: /*@
428:   PetscConvEstRateView - Displays the convergence rate to a viewer

430:    Collective on SNES

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

437:    Options Database Keys:
438: .  -snes_convergence_estimate - print the convergence rate

440:    Level: developer

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

450:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
451:   if (isAscii) {
452:     PetscInt Nf = ce->Nf, f;

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

468: /*@
469:   PetscConvEstCreate - Create a PetscConvEst object

471:   Collective

473:   Input Parameter:
474: . comm - The communicator for the PetscConvEst object

476:   Output Parameter:
477: . ce   - The PetscConvEst object

479:   Level: beginner

481: .seealso: PetscConvEstDestroy(), PetscConvEstGetConvRate()
482: @*/
483: PetscErrorCode PetscConvEstCreate(MPI_Comm comm, PetscConvEst *ce)
484: {

489:   PetscSysInitializePackage();
490:   PetscHeaderCreate(*ce, PETSC_OBJECT_CLASSID, "PetscConvEst", "ConvergenceEstimator", "SNES", comm, PetscConvEstDestroy, PetscConvEstView);
491:   (*ce)->monitor = PETSC_FALSE;
492:   (*ce)->r       = 2.0;
493:   (*ce)->Nr      = 4;
494:   (*ce)->event   = -1;
495:   (*ce)->ops->setsolver    = PetscConvEstSetSNES_Private;
496:   (*ce)->ops->initguess    = PetscConvEstInitGuessSNES_Private;
497:   (*ce)->ops->computeerror = PetscConvEstComputeErrorSNES_Private;
498:   (*ce)->ops->getconvrate  = PetscConvEstGetConvRateSNES_Private;
499:   return(0);
500: }