Actual source code: convest.c

petsc-3.10.5 2019-03-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:   PetscDS        prob;
295:   PetscObject    disc;
296:   MPI_Comm       comm;
297:   const char    *uname, *dmname;
298:   void          *ctx;
299:   Vec            u;
300:   PetscReal      t = 0.0, *x, *y, slope, intercept;
301:   PetscInt      *dof, dim, Nr = ce->Nr, r, f;
302:   PetscLogEvent  event;

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

320:     PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
321:     PetscLogStageRegister(stageName, &stage);
322:     PetscLogStagePush(stage);
323:     if (r > 0) {
324:       DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
325:       DMSetCoarseDM(dm[r], dm[r-1]);
326:       DMSetDS(dm[r], prob);
327:       PetscObjectGetName((PetscObject) dm[r-1], &dmname);
328:       PetscObjectSetName((PetscObject) dm[r], dmname);
329:       for (f = 0; f <= ce->Nf; ++f) {
330:         PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
331:         DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
332:         DMSetNullSpaceConstructor(dm[r],   f,  nspconstr);
333:       }
334:     }
335:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
336:     /* Create solution */
337:     DMCreateGlobalVector(dm[r], &u);
338:     PetscDSGetDiscretization(prob, 0, &disc);
339:     PetscObjectGetName(disc, &uname);
340:     PetscObjectSetName((PetscObject) u, uname);
341:     /* Setup solver */
342:     SNESReset(ce->snes);
343:     SNESSetDM(ce->snes, dm[r]);
344:     DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
345:     SNESSetFromOptions(ce->snes);
346:     /* Create initial guess */
347:     DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
348:     SNESSolve(ce->snes, NULL, u);
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:       DMPlexGetHeightStratum(dm[r], 0, NULL, &dof[r*ce->Nf+f]);
353:     }
354:     /* Monitor */
355:     if (ce->monitor) {
356:       PetscReal *errors = &ce->errors[r*ce->Nf];

358:       PetscPrintf(comm, "L_2 Error: ");
359:       if (ce->Nf > 1) {PetscPrintf(comm, "[");}
360:       for (f = 0; f < ce->Nf; ++f) {
361:         if (f > 0) {PetscPrintf(comm, ", ");}
362:         if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
363:         else                     {PetscPrintf(comm, "%g", (double)errors[f]);}
364:       }
365:       if (ce->Nf > 1) {PetscPrintf(comm, "]");}
366:       PetscPrintf(comm, "\n");
367:     }
368:     /* Cleanup */
369:     VecDestroy(&u);
370:     PetscLogStagePop();
371:   }
372:   for (r = 1; r <= Nr; ++r) {
373:     DMDestroy(&dm[r]);
374:   }
375:   /* Fit convergence rate */
376:   PetscMalloc2(Nr+1, &x, Nr+1, &y);
377:   for (f = 0; f < ce->Nf; ++f) {
378:     for (r = 0; r <= Nr; ++r) {
379:       x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
380:       y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
381:     }
382:     PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
383:     /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
384:     alpha[f] = -slope * dim;
385:   }
386:   PetscFree2(x, y);
387:   PetscFree2(dm, dof);
388:   /* Restore solver */
389:   SNESReset(ce->snes);
390:   SNESSetDM(ce->snes, ce->idm);
391:   DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
392:   SNESSetFromOptions(ce->snes);
393:   return(0);
394: }

396: /*@
397:   PetscConvEstRateView - Displays the convergence rate to a viewer

399:    Collective on SNES

401:    Parameter:
402: +  snes - iterative context obtained from SNESCreate()
403: .  alpha - the convergence rate for each field
404: -  viewer - the viewer to display the reason

406:    Options Database Keys:
407: .  -snes_convergence_estimate - print the convergence rate

409:    Level: developer

411: .seealso: PetscConvEstGetRate()
412: @*/
413: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, PetscReal alpha[], PetscViewer viewer)
414: {
415:   PetscBool      isAscii;
416:   PetscInt       f;

420:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
421:   if (isAscii) {
422:     PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
423:     PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: ");
424:     if (ce->Nf > 1) {PetscViewerASCIIPrintf(viewer, "[");}
425:     for (f = 0; f < ce->Nf; ++f) {
426:       if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
427:       PetscViewerASCIIPrintf(viewer, "%g", (double) alpha[f]);
428:     }
429:     if (ce->Nf > 1) {PetscViewerASCIIPrintf(viewer, "]");}
430:     PetscViewerASCIIPrintf(viewer, "\n");
431:     PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
432:   }
433:   return(0);
434: }