Actual source code: convest.c

petsc-3.8.4 2018-03-24
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("-num_refine", "The number of refinements for the convergence check", "PetscConvEst", ce->Nr, &ce->Nr, NULL);
 95:   PetscOptionsEnd();
 96:   return(0);
 97: }

 99: /*@
100:   PetscConvEstView - Views a PetscConvEst object

102:   Collective on PetscConvEst

104:   Input Parameters:
105: + ce     - The PetscConvEst object
106: - viewer - The PetscViewer object

108:   Level: beginner

110: .keywords: PetscConvEst, convergence, view
111: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
112: @*/
113: PetscErrorCode PetscConvEstView(PetscConvEst ce, PetscViewer viewer)
114: {

118:   PetscObjectPrintClassNamePrefixType((PetscObject) ce, viewer);
119:   PetscViewerASCIIPrintf(viewer, "ConvEst with %D levels\n", ce->Nr+1);
120:   return(0);
121: }

123: /*@
124:   PetscConvEstGetSolver - Gets the solver used to produce discrete solutions

126:   Not collective

128:   Input Parameter:
129: . ce   - The PetscConvEst object

131:   Output Parameter:
132: . snes - The solver

134:   Level: intermediate

136: .keywords: PetscConvEst, convergence
137: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
138: @*/
139: PetscErrorCode PetscConvEstGetSolver(PetscConvEst ce, SNES *snes)
140: {
144:   *snes = ce->snes;
145:   return(0);
146: }

148: /*@
149:   PetscConvEstSetSolver - Sets the solver used to produce discrete solutions

151:   Not collective

153:   Input Parameters:
154: + ce   - The PetscConvEst object
155: - snes - The solver

157:   Level: intermediate

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

161: .keywords: PetscConvEst, convergence
162: .seealso: PetscConvEstGetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
163: @*/
164: PetscErrorCode PetscConvEstSetSolver(PetscConvEst ce, SNES snes)
165: {

171:   ce->snes = snes;
172:   SNESGetDM(ce->snes, &ce->idm);
173:   return(0);
174: }

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

179:   Collective on PetscConvEst

181:   Input Parameters:
182: . ce - The PetscConvEst object

184:   Level: beginner

186: .keywords: PetscConvEst, convergence, setup
187: .seealso: PetscConvEstCreate(), PetscConvEstGetConvRate()
188: @*/
189: PetscErrorCode PetscConvEstSetUp(PetscConvEst ce)
190: {
191:   PetscDS        prob;
192:   PetscInt       f;

196:   DMGetDS(ce->idm, &prob);
197:   PetscDSGetNumFields(prob, &ce->Nf);
198:   PetscMalloc1((ce->Nr+1)*ce->Nf, &ce->errors);
199:   PetscMalloc2(ce->Nf, &ce->initGuess, ce->Nf, &ce->exactSol);
200:   for (f = 0; f < ce->Nf; ++f) ce->initGuess[f] = zero_private;
201:   for (f = 0; f < ce->Nf; ++f) {
202:     PetscDSGetExactSolution(prob, f, &ce->exactSol[f]);
203:     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);
204:   }
205:   return(0);
206: }

208: static PetscErrorCode PetscConvEstLinearRegression_Private(PetscConvEst ce, PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept)
209: {
210:   PetscScalar    H[4];
211:   PetscReal     *X, *Y, beta[2];
212:   PetscInt       i, j, k;

216:   *slope = *intercept = 0.0;
217:   PetscMalloc2(n*2, &X, n*2, &Y);
218:   for (k = 0; k < n; ++k) {
219:     /* X[n,2] = [1, x] */
220:     X[k*2+0] = 1.0;
221:     X[k*2+1] = x[k];
222:   }
223:   /* H = X^T X */
224:   for (i = 0; i < 2; ++i) {
225:     for (j = 0; j < 2; ++j) {
226:       H[i*2+j] = 0.0;
227:       for (k = 0; k < n; ++k) {
228:         H[i*2+j] += X[k*2+i] * X[k*2+j];
229:       }
230:     }
231:   }
232:   /* H = (X^T X)^{-1} */
233:   {
234:     PetscBLASInt two = 2, ipiv[2], info;
235:     PetscScalar  work[2];

237:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
238:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
239:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
240:     PetscFPTrapPop();
241:   }
242:     /* Y = H X^T */
243:   for (i = 0; i < 2; ++i) {
244:     for (k = 0; k < n; ++k) {
245:       Y[i*n+k] = 0.0;
246:       for (j = 0; j < 2; ++j) {
247:         Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j];
248:       }
249:     }
250:   }
251:   /* beta = Y error = [y-intercept, slope] */
252:   for (i = 0; i < 2; ++i) {
253:     beta[i] = 0.0;
254:     for (k = 0; k < n; ++k) {
255:       beta[i] += Y[i*n+k] * y[k];
256:     }
257:   }
258:   PetscFree2(X, Y);
259:   *intercept = beta[0];
260:   *slope     = beta[1];
261:   return(0);
262: }

264: /*@
265:   PetscConvEstGetConvRate - Returns an estimate of the convergence rate for the discretization

267:   Not collective

269:   Input Parameter:
270: . ce   - The PetscConvEst object

272:   Output Parameter:
273: . alpha - The convergence rate

275:   Note: The convergence rate alpha is defined by
276: $ || u_h - u_exact || < C h^alpha
277: where u_h is the discrete solution, and h is a measure of the discretization size.

279: We solve a series of problems on refined meshes, calculate an error based upon the exact solution in the DS,
280: and then fit the result to our model above using linear regression.

282:   Options database keys:
283: . -snes_convergence_estimate : Execute convergence estimation and print out the rate

285:   Level: intermediate

287: .keywords: PetscConvEst, convergence
288: .seealso: PetscConvEstSetSolver(), PetscConvEstCreate(), PetscConvEstGetConvRate()
289: @*/
290: PetscErrorCode PetscConvEstGetConvRate(PetscConvEst ce, PetscReal *alpha)
291: {
292:   DM            *dm;
293:   PetscDS        prob;
294:   PetscObject    disc;
295:   MPI_Comm       comm;
296:   const char    *uname, *dmname;
297:   void          *ctx;
298:   Vec            u;
299:   PetscReal      t = 0.0, *x, *y, slope, intercept;
300:   PetscInt      *dof, dim, Nr = ce->Nr, r;

304:   PetscObjectGetComm((PetscObject) ce, &comm);
305:   DMGetDimension(ce->idm, &dim);
306:   DMGetApplicationContext(ce->idm, &ctx);
307:   DMGetDS(ce->idm, &prob);
308:   DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
309:   PetscMalloc2((Nr+1), &dm, (Nr+1), &dof);
310:   dm[0]  = ce->idm;
311:   *alpha = 0.0;
312:   /* Loop over meshes */
313:   for (r = 0; r <= Nr; ++r) {
314:     if (r > 0) {
315:       DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
316:       DMSetCoarseDM(dm[r], dm[r-1]);
317:       DMSetDS(dm[r], prob);
318:       PetscObjectGetName((PetscObject) dm[r-1], &dmname);
319:       PetscObjectSetName((PetscObject) dm[r], dmname);
320:     }
321:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
322:     DMPlexGetHeightStratum(dm[r], 0, NULL, &dof[r]);
323:     /* Create solution */
324:     DMCreateGlobalVector(dm[r], &u);
325:     PetscDSGetDiscretization(prob, 0, &disc);
326:     PetscObjectGetName(disc, &uname);
327:     PetscObjectSetName((PetscObject) u, uname);
328:     /* Setup solver */
329:     SNESReset(ce->snes);
330:     SNESSetDM(ce->snes, dm[r]);
331:     DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);
332:     SNESSetFromOptions(ce->snes);
333:     /* Create initial guess */
334:     DMProjectFunction(dm[r], t, ce->initGuess, NULL, INSERT_VALUES, u);
335:     SNESSolve(ce->snes, NULL, u);
336:     DMComputeL2FieldDiff(dm[r], t, ce->exactSol, NULL, u, &ce->errors[r*ce->Nf]);
337:     /* Monitor */
338:     if (ce->monitor) {
339:       PetscReal *errors = &ce->errors[r*ce->Nf];
340:       PetscInt f;

342:       PetscPrintf(comm, "L_2 Error: [");
343:       for (f = 0; f < ce->Nf; ++f) {
344:         if (f > 0) {PetscPrintf(comm, ", ");}
345:         if (errors[f] < 1.0e-11) {PetscPrintf(comm, "< 1e-11");}
346:         else                     {PetscPrintf(comm, "%g", (double)errors[f]);}
347:       }
348:       PetscPrintf(comm, "]\n");
349:     }
350:     /* Cleanup */
351:     VecDestroy(&u);
352:   }
353:   for (r = 1; r <= Nr; ++r) {
354:     DMDestroy(&dm[r]);
355:   }
356:   /* Fit convergence rate */
357:   PetscMalloc2(Nr+1, &x, Nr+1, &y);
358:   for (r = 0; r <= Nr; ++r) {
359:     x[r] = PetscLog10Real(dof[r]);
360:     y[r] = PetscLog10Real(ce->errors[r*ce->Nf+0]);
361:   }
362:   PetscConvEstLinearRegression_Private(ce, Nr+1, x, y, &slope, &intercept);
363:   PetscFree2(x, y);
364:   /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
365:   *alpha = -slope * dim;
366:   PetscFree2(dm, dof);
367:   /* Restore solver */
368:   SNESReset(ce->snes);
369:   SNESSetDM(ce->snes, ce->idm);
370:   DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);
371:   SNESSetFromOptions(ce->snes);
372:   return(0);
373: }

375: /*@
376:   PetscConvEstRateView - Displays the convergence rate to a viewer

378:    Collective on SNES

380:    Parameter:
381: +  snes - iterative context obtained from SNESCreate()
382: .  alpha - the convergence rate
383: -  viewer - the viewer to display the reason

385:    Options Database Keys:
386: .  -snes_convergence_estimate - print the convergence rate

388:    Level: developer

390: .seealso: PetscConvEstGetRate()
391: @*/
392: PetscErrorCode PetscConvEstRateView(PetscConvEst ce, PetscReal alpha, PetscViewer viewer)
393: {
394:   PetscBool      isAscii;

398:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
399:   if (isAscii) {
400:     PetscViewerASCIIAddTab(viewer, ((PetscObject) ce)->tablevel);
401:     PetscViewerASCIIPrintf(viewer, "L_2 convergence rate: %g\n", (double) alpha);
402:     PetscViewerASCIISubtractTab(viewer, ((PetscObject) ce)->tablevel);
403:   }
404:   return(0);
405: }