Actual source code: alpha.c

petsc-3.3-p7 2013-05-11
  1: /*
  2:   Code for timestepping with implicit generalized-\alpha method
  3:   for first order systems.
  4: */
  5: #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/

  7: typedef PetscErrorCode (*TSAlphaAdaptFunction)(TS,PetscReal,Vec,Vec,PetscReal*,PetscBool*,void*);

  9: typedef struct {
 10:   Vec X0,Xa,X1;
 11:   Vec V0,Va,V1;
 12:   Vec R,E;
 13:   PetscReal Alpha_m;
 14:   PetscReal Alpha_f;
 15:   PetscReal Gamma;
 16:   PetscReal stage_time;
 17:   PetscReal shift;

 19:   TSAlphaAdaptFunction adapt;
 20:   void *adaptctx;
 21:   PetscReal rtol;
 22:   PetscReal atol;
 23:   PetscReal rho;
 24:   PetscReal scale_min;
 25:   PetscReal scale_max;
 26:   PetscReal dt_min;
 27:   PetscReal dt_max;
 28: } TS_Alpha;

 32: static PetscErrorCode TSStep_Alpha(TS ts)
 33: {
 34:   TS_Alpha            *th    = (TS_Alpha*)ts->data;
 35:   PetscInt            its,lits,reject;
 36:   PetscReal           next_time_step;
 37:   SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
 38:   PetscErrorCode      ierr;

 41:   if (ts->steps == 0) {
 42:     VecSet(th->V0,0.0);
 43:   } else {
 44:     VecCopy(th->V1,th->V0);
 45:   }
 46:   VecCopy(ts->vec_sol,th->X0);
 47:   next_time_step = ts->time_step;
 48:   for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
 49:     ts->time_step = next_time_step;
 50:     th->stage_time = ts->ptime + th->Alpha_f*ts->time_step;
 51:     th->shift = th->Alpha_m/(th->Alpha_f*th->Gamma*ts->time_step);
 52:     TSPreStep(ts);
 53:     TSPreStage(ts,th->stage_time);
 54:     /* predictor */
 55:     VecCopy(th->X0,th->X1);
 56:     /* solve R(X,V) = 0 */
 57:     SNESSolve(ts->snes,PETSC_NULL,th->X1);
 58:     /* V1 = (1-1/Gamma)*V0 + 1/(Gamma*dT)*(X1-X0) */
 59:     VecWAXPY(th->V1,-1,th->X0,th->X1);
 60:     VecAXPBY(th->V1,1-1/th->Gamma,1/(th->Gamma*ts->time_step),th->V0);
 61:     /* nonlinear solve convergence */
 62:     SNESGetConvergedReason(ts->snes,&snesreason);
 63:     if (snesreason < 0 && !th->adapt) break;
 64:     SNESGetIterationNumber(ts->snes,&its);
 65:     SNESGetLinearSolveIterations(ts->snes,&lits);
 66:     ts->snes_its += its; ts->ksp_its += lits;
 67:     PetscInfo3(ts,"step=%D, nonlinear solve iterations=%D, linear solve iterations=%D\n",ts->steps,its,lits);
 68:     /* time step adaptativity */
 69:     if (!th->adapt) break;
 70:     else {
 71:       PetscReal t1 = ts->ptime + ts->time_step;
 72:       PetscBool stepok = (reject==0) ? PETSC_TRUE : PETSC_FALSE;
 73:       th->adapt(ts,t1,th->X1,th->V1,&next_time_step,&stepok,th->adaptctx);
 74:       PetscInfo5(ts,"Step %D (t=%G,dt=%G) %s, next dt=%G\n",ts->steps,ts->ptime,ts->time_step,stepok?"accepted":"rejected",next_time_step);
 75:       if (stepok) break;
 76:     }
 77:   }
 78:   if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
 79:     ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
 80:     PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
 81:     return(0);
 82:   }
 83:   if (reject >= ts->max_reject) {
 84:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 85:     PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
 86:     return(0);
 87:   }
 88:   VecCopy(th->X1,ts->vec_sol);
 89:   ts->ptime += ts->time_step;
 90:   ts->time_step = next_time_step;
 91:   ts->steps++;
 92:   return(0);
 93: }

 97: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X)
 98: {
 99:   TS_Alpha       *th = (TS_Alpha*)ts->data;
100:   PetscReal      dt = t - ts->ptime;

104:   VecCopy(ts->vec_sol,X);
105:   VecAXPY(X,th->Gamma*dt,th->V1);
106:   VecAXPY(X,(1-th->Gamma)*dt,th->V0);
107:   return(0);
108: }

110: /*------------------------------------------------------------*/
113: static PetscErrorCode TSReset_Alpha(TS ts)
114: {
115:   TS_Alpha       *th = (TS_Alpha*)ts->data;
116:   PetscErrorCode  ierr;

119:   VecDestroy(&th->X0);
120:   VecDestroy(&th->Xa);
121:   VecDestroy(&th->X1);
122:   VecDestroy(&th->V0);
123:   VecDestroy(&th->Va);
124:   VecDestroy(&th->V1);
125:   VecDestroy(&th->E);
126:   return(0);
127: }

131: static PetscErrorCode TSDestroy_Alpha(TS ts)
132: {
133:   PetscErrorCode  ierr;

136:   TSReset_Alpha(ts);
137:   PetscFree(ts->data);

139:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetRadius_C","",PETSC_NULL);
140:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetParams_C","",PETSC_NULL);
141:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaGetParams_C","",PETSC_NULL);
142:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetAdapt_C","",PETSC_NULL);
143:   return(0);
144: }

148: static PetscErrorCode SNESTSFormFunction_Alpha(SNES snes,Vec x,Vec y,TS ts)
149: {
150:   TS_Alpha       *th = (TS_Alpha*)ts->data;
151:   Vec            X0 = th->X0, V0 = th->V0;
152:   Vec            X1 = x, V1 = th->V1, R = y;

156:   /* V1 = (1-1/Gamma)*V0 + 1/(Gamma*dT)*(X1-X0) */
157:   VecWAXPY(V1,-1,X0,X1);
158:   VecAXPBY(V1,1-1/th->Gamma,1/(th->Gamma*ts->time_step),V0);
159:   /* Xa = X0 + Alpha_f*(X1-X0) */
160:   VecWAXPY(th->Xa,-1,X0,X1);
161:   VecAYPX(th->Xa,th->Alpha_f,X0);
162:   /* Va = V0 + Alpha_m*(V1-V0) */
163:   VecWAXPY(th->Va,-1,V0,V1);
164:   VecAYPX(th->Va,th->Alpha_m,V0);
165:   /* F = Function(ta,Xa,Va) */
166:   TSComputeIFunction(ts,th->stage_time,th->Xa,th->Va,R,PETSC_FALSE);
167:   VecScale(R,1/th->Alpha_f);
168:   return(0);
169: }

173: static PetscErrorCode SNESTSFormJacobian_Alpha(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts)
174: {
175:   TS_Alpha       *th = (TS_Alpha*)ts->data;

179:   /* A,B = Jacobian(ta,Xa,Va) */
180:   TSComputeIJacobian(ts,th->stage_time,th->Xa,th->Va,th->shift,A,B,str,PETSC_FALSE);
181:   return(0);
182: }

186: static PetscErrorCode TSSetUp_Alpha(TS ts)
187: {
188:   TS_Alpha *th = (TS_Alpha*)ts->data;

192:   VecDuplicate(ts->vec_sol,&th->X0);
193:   VecDuplicate(ts->vec_sol,&th->Xa);
194:   VecDuplicate(ts->vec_sol,&th->X1);
195:   VecDuplicate(ts->vec_sol,&th->V0);
196:   VecDuplicate(ts->vec_sol,&th->Va);
197:   VecDuplicate(ts->vec_sol,&th->V1);
198:   return(0);
199: }

203: static PetscErrorCode TSSetFromOptions_Alpha(TS ts)
204: {
205:   TS_Alpha *th = (TS_Alpha*)ts->data;

209:   PetscOptionsHead("Alpha ODE solver options");
210:   {
211:     PetscBool flag, adapt = PETSC_FALSE;
212:     PetscReal radius = 1.0;
213:     PetscOptionsReal("-ts_alpha_radius","spectral radius","TSAlphaSetRadius",radius,&radius,&flag);
214:     if (flag) { TSAlphaSetRadius(ts,radius); }
215:     PetscOptionsReal("-ts_alpha_alpha_m","algoritmic parameter alpha_m","TSAlphaSetParams",th->Alpha_m,&th->Alpha_m,PETSC_NULL);
216:     PetscOptionsReal("-ts_alpha_alpha_f","algoritmic parameter alpha_f","TSAlphaSetParams",th->Alpha_f,&th->Alpha_f,PETSC_NULL);
217:     PetscOptionsReal("-ts_alpha_gamma","algoritmic parameter gamma","TSAlphaSetParams",th->Gamma,&th->Gamma,PETSC_NULL);
218:     TSAlphaSetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma);

220:     PetscOptionsBool("-ts_alpha_adapt","default time step adaptativity","TSAlphaSetAdapt",adapt,&adapt,&flag);
221:     if (flag) { TSAlphaSetAdapt(ts,adapt?TSAlphaAdaptDefault:PETSC_NULL,PETSC_NULL);  }
222:     PetscOptionsReal("-ts_alpha_adapt_rtol","relative tolerance for dt adaptativity","",th->rtol,&th->rtol,PETSC_NULL);
223:     PetscOptionsReal("-ts_alpha_adapt_atol","absolute tolerance for dt adaptativity","",th->atol,&th->atol,PETSC_NULL);
224:     PetscOptionsReal("-ts_alpha_adapt_min","minimum dt scale","",th->scale_min,&th->scale_min,PETSC_NULL);
225:     PetscOptionsReal("-ts_alpha_adapt_max","maximum dt scale","",th->scale_max,&th->scale_max,PETSC_NULL);
226:     PetscOptionsReal("-ts_alpha_adapt_dt_min","minimum dt","",th->dt_min,&th->dt_min,PETSC_NULL);
227:     PetscOptionsReal("-ts_alpha_adapt_dt_max","maximum dt","",th->dt_max,&th->dt_max,PETSC_NULL);
228:     SNESSetFromOptions(ts->snes);
229:   }
230:   PetscOptionsTail();
231:   return(0);
232: }

236: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
237: {
238:   TS_Alpha       *th = (TS_Alpha*)ts->data;
239:   PetscBool       iascii;
240:   PetscErrorCode  ierr;

243:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
244:   if (iascii) {
245:     PetscViewerASCIIPrintf(viewer,"  Alpha_m=%G, Alpha_f=%G, Gamma=%G\n",th->Alpha_m,th->Alpha_f,th->Gamma);
246:   }
247:   SNESView(ts->snes,viewer);
248:   return(0);
249: }

251: /*------------------------------------------------------------*/

253: EXTERN_C_BEGIN

257: PetscErrorCode  TSAlphaSetRadius_Alpha(TS ts,PetscReal radius)
258: {
259:   TS_Alpha *th = (TS_Alpha*)ts->data;

262:   if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %G not in range [0,1]",radius);
263:   th->Alpha_m = 0.5*(3-radius)/(1+radius);
264:   th->Alpha_f = 1/(1+radius);
265:   th->Gamma   = 0.5 + th->Alpha_m - th->Alpha_f;
266:   return(0);
267: }

271: PetscErrorCode  TSAlphaSetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
272: {
273:   TS_Alpha *th = (TS_Alpha*)ts->data;

276:   th->Alpha_m = alpha_m;
277:   th->Alpha_f = alpha_f;
278:   th->Gamma   = gamma;
279:   return(0);
280: }

284: PetscErrorCode  TSAlphaGetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
285: {
286:   TS_Alpha *th = (TS_Alpha*)ts->data;

289:   if (alpha_m) *alpha_m = th->Alpha_m;
290:   if (alpha_f) *alpha_f = th->Alpha_f;
291:   if (gamma)   *gamma   = th->Gamma;
292:   return(0);
293: }

297: PetscErrorCode  TSAlphaSetAdapt_Alpha(TS ts,TSAlphaAdaptFunction adapt,void *ctx)
298: {
299:   TS_Alpha *th = (TS_Alpha*)ts->data;

302:   th->adapt    = adapt;
303:   th->adaptctx = ctx;
304:   return(0);
305: }

307: EXTERN_C_END

309: /* ------------------------------------------------------------ */
310: /*MC
311:       TSALPHA - DAE solver using the implicit Generalized-Alpha method

313:   Level: beginner

315:   References:
316:   K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
317:   method for integrating the filtered Navier-Stokes equations with a
318:   stabilized finite element method", Computer Methods in Applied
319:   Mechanics and Engineering, 190, 305-319, 2000.
320:   DOI: 10.1016/S0045-7825(00)00203-6.

322:   J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
323:   Dynamics with Improved Numerical Dissipation: The Generalized-alpha
324:   Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.

326: .seealso:  TSCreate(), TS, TSSetType()

328: M*/
329: EXTERN_C_BEGIN
332: PetscErrorCode  TSCreate_Alpha(TS ts)
333: {
334:   TS_Alpha       *th;

338:   ts->ops->reset          = TSReset_Alpha;
339:   ts->ops->destroy        = TSDestroy_Alpha;
340:   ts->ops->view           = TSView_Alpha;
341:   ts->ops->setup          = TSSetUp_Alpha;
342:   ts->ops->step           = TSStep_Alpha;
343:   ts->ops->interpolate    = TSInterpolate_Alpha;
344:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
345:   ts->ops->snesfunction   = SNESTSFormFunction_Alpha;
346:   ts->ops->snesjacobian   = SNESTSFormJacobian_Alpha;

348:   PetscNewLog(ts,TS_Alpha,&th);
349:   ts->data = (void*)th;

351:   th->Alpha_m = 0.5;
352:   th->Alpha_f = 0.5;
353:   th->Gamma   = 0.5;

355:   th->rtol      = 1e-3;
356:   th->atol      = 1e-3;
357:   th->rho       = 0.9;
358:   th->scale_min = 0.1;
359:   th->scale_max = 5.0;
360:   th->dt_min    = 0.0;
361:   th->dt_max    = PETSC_MAX_REAL;

363:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetAdapt_C","TSAlphaSetAdapt_Alpha",TSAlphaSetAdapt_Alpha);
364:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetRadius_C","TSAlphaSetRadius_Alpha",TSAlphaSetRadius_Alpha);
365:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetParams_C","TSAlphaSetParams_Alpha",TSAlphaSetParams_Alpha);
366:   PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaGetParams_C","TSAlphaGetParams_Alpha",TSAlphaGetParams_Alpha);

368:   return(0);
369: }
370: EXTERN_C_END

374: /*@C
375:   TSAlphaSetAdapt - sets the time step adaptativity and acceptance test routine

377:   This function allows to accept/reject a step and select the
378:   next time step to use.

380:   Not Collective

382:   Input Parameter:
383: +  ts - timestepping context
384: .  adapt - user-defined adapt routine
385: -  ctx  - [optional] user-defined context for private data for the
386:          adapt routine (may be PETSC_NULL)

388:    Calling sequence of adapt:
389: $    adapt (TS ts,PetscReal t,Vec X,Vec Xdot,
390: $            PetscReal *next_dt,PetscBool *accepted,void *ctx);

392:   Level: intermediate

394: @*/
395: PetscErrorCode  TSAlphaSetAdapt(TS ts,TSAlphaAdaptFunction adapt,void *ctx)
396: {
400:   PetscTryMethod(ts,"TSAlphaSetAdapt_C",(TS,TSAlphaAdaptFunction,void*),(ts,adapt,ctx));
401:   return(0);
402: }

406: PetscErrorCode  TSAlphaAdaptDefault(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
407: {
408:   TS_Alpha            *th;
409:   SNESConvergedReason snesreason;
410:   PetscReal           dt,normX,normE,Emax,scale;
411:   PetscErrorCode      ierr;

415: #if PETSC_USE_DEBUG
416:   {
417:     PetscBool match;
418:     PetscObjectTypeCompare((PetscObject)ts,TSALPHA,&match);
419:     if (!match) SETERRQ(((PetscObject)ts)->comm,1,"Only for TSALPHA");
420:   }
421: #endif
422:   th = (TS_Alpha*)ts->data;

424:   SNESGetConvergedReason(ts->snes,&snesreason);
425:   if (snesreason < 0) {
426:     *ok = PETSC_FALSE;
427:     *nextdt *= th->scale_min;
428:     goto finally;
429:   }

431:   /* first-order aproximation to the local error */
432:   /* E = (X0 + dt*Xdot) - X */
433:   TSGetTimeStep(ts,&dt);
434:   if (!th->E) {VecDuplicate(th->X0,&th->E);}
435:   VecWAXPY(th->E,dt,Xdot,th->X0);
436:   VecAXPY(th->E,-1,X);
437:   VecNorm(th->E,NORM_2,&normE);
438:   /* compute maximum allowable error */
439:   VecNorm(X,NORM_2,&normX);
440:   if (normX == 0) {VecNorm(th->X0,NORM_2,&normX);}
441:   Emax =  th->rtol * normX + th->atol;
442:   /* compute next time step */
443:   if (normE > 0) {
444:     scale = th->rho * PetscRealPart(PetscSqrtScalar((PetscScalar)(Emax/normE)));
445:     scale = PetscMax(scale,th->scale_min);
446:     scale = PetscMin(scale,th->scale_max);
447:     if (!(*ok))
448:       scale = PetscMin(1.0,scale);
449:     *nextdt *= scale;
450:   }
451:   /* accept or reject step */
452:   if (normE <= Emax)
453:     *ok = PETSC_TRUE;
454:   else
455:     *ok = PETSC_FALSE;

457:   finally:
458:   *nextdt = PetscMax(*nextdt,th->dt_min);
459:   *nextdt = PetscMin(*nextdt,th->dt_max);
460:   return(0);
461: }

465: /*@
466:   TSAlphaSetRadius - sets the desired spectral radius of the method
467:                      (i.e. high-frequency numerical damping)

469:   Logically Collective on TS

471:   The algorithmic parameters \alpha_m and \alpha_f of the
472:   generalized-\alpha method can be computed in terms of a specified
473:   spectral radius \rho in [0,1] for infinite time step in order to
474:   control high-frequency numerical damping:
475:     alpha_m = 0.5*(3-\rho)/(1+\rho)
476:     alpha_f = 1/(1+\rho)

478:   Input Parameter:
479: +  ts - timestepping context
480: -  radius - the desired spectral radius

482:   Options Database:
483: .  -ts_alpha_radius <radius>

485:   Level: intermediate

487: .seealso: TSAlphaSetParams(), TSAlphaGetParams()
488: @*/
489: PetscErrorCode  TSAlphaSetRadius(TS ts,PetscReal radius)
490: {

495:   PetscTryMethod(ts,"TSAlphaSetRadius_C",(TS,PetscReal),(ts,radius));
496:   return(0);
497: }

501: /*@
502:   TSAlphaSetParams - sets the algorithmic parameters for TSALPHA

504:   Not Collective

506:   Second-order accuracy can be obtained so long as:
507:     \gamma = 0.5 + alpha_m - alpha_f

509:   Unconditional stability requires:
510:     \alpha_m >= \alpha_f >= 0.5

512:   Backward Euler method is recovered when:
513:     \alpha_m = \alpha_f = gamma = 1


516:   Input Parameter:
517: +  ts - timestepping context
518: .  \alpha_m - algorithmic paramenter
519: .  \alpha_f - algorithmic paramenter
520: -  \gamma   - algorithmic paramenter

522:    Options Database:
523: +  -ts_alpha_alpha_m <alpha_m>
524: .  -ts_alpha_alpha_f <alpha_f>
525: -  -ts_alpha_gamma <gamma>

527:   Note:
528:   Use of this function is normally only required to hack TSALPHA to
529:   use a modified integration scheme. Users should call
530:   TSAlphaSetRadius() to set the desired spectral radius of the methods
531:   (i.e. high-frequency damping) in order so select optimal values for
532:   these parameters.

534:   Level: advanced

536: .seealso: TSAlphaSetRadius(), TSAlphaGetParams()
537: @*/
538: PetscErrorCode  TSAlphaSetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
539: {

544:   PetscTryMethod(ts,"TSAlphaSetParams_C",(TS,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma));
545:   return(0);
546: }

550: /*@
551:   TSAlphaGetParams - gets the algorithmic parameters for TSALPHA

553:   Not Collective

555:   Input Parameter:
556: +  ts - timestepping context
557: .  \alpha_m - algorithmic parameter
558: .  \alpha_f - algorithmic parameter
559: -  \gamma   - algorithmic parameter

561:   Note:
562:   Use of this function is normally only required to hack TSALPHA to
563:   use a modified integration scheme. Users should call
564:   TSAlphaSetRadius() to set the high-frequency damping (i.e. spectral
565:   radius of the method) in order so select optimal values for these
566:   parameters.

568:   Level: advanced

570: .seealso: TSAlphaSetRadius(), TSAlphaSetParams()
571: @*/
572: PetscErrorCode  TSAlphaGetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
573: {

581:   PetscUseMethod(ts,"TSAlphaGetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma));
582:   return(0);
583: }