Actual source code: convest.c

petsc-master 2019-09-15
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: }


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