Actual source code: alpha2.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: /*
  2:   Code for timestepping with implicit generalized-\alpha method
  3:   for second order systems.
  4: */
  5: #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/

  7: static PetscBool  cited = PETSC_FALSE;
  8: static const char citation[] =
  9:   "@article{Chung1993,\n"
 10:   "  title   = {A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation: The Generalized-$\\alpha$ Method},\n"
 11:   "  author  = {J. Chung, G. M. Hubert},\n"
 12:   "  journal = {ASME Journal of Applied Mechanics},\n"
 13:   "  volume  = {60},\n"
 14:   "  number  = {2},\n"
 15:   "  pages   = {371--375},\n"
 16:   "  year    = {1993},\n"
 17:   "  issn    = {0021-8936},\n"
 18:   "  doi     = {http://dx.doi.org/10.1115/1.2900803}\n}\n";

 20: typedef struct {
 21:   PetscReal stage_time;
 22:   PetscReal shift_V;
 23:   PetscReal shift_A;
 24:   PetscReal scale_F;
 25:   Vec       X0,Xa,X1;
 26:   Vec       V0,Va,V1;
 27:   Vec       A0,Aa,A1;

 29:   Vec       vec_dot;

 31:   PetscReal Alpha_m;
 32:   PetscReal Alpha_f;
 33:   PetscReal Gamma;
 34:   PetscReal Beta;
 35:   PetscInt  order;

 37:   PetscBool adapt;
 38:   Vec       vec_sol_prev;
 39:   Vec       vec_dot_prev;
 40:   Vec       vec_lte_work[2];

 42:   TSStepStatus status;
 43: } TS_Alpha;

 47: static PetscErrorCode TSAlpha_StageTime(TS ts)
 48: {
 49:   TS_Alpha  *th = (TS_Alpha*)ts->data;
 50:   PetscReal t  = ts->ptime;
 51:   PetscReal dt = ts->time_step;
 52:   PetscReal Alpha_m = th->Alpha_m;
 53:   PetscReal Alpha_f = th->Alpha_f;
 54:   PetscReal Gamma   = th->Gamma;
 55:   PetscReal Beta    = th->Beta;

 58:   th->stage_time = t + Alpha_f*dt;
 59:   th->shift_V = Gamma/(dt*Beta);
 60:   th->shift_A = Alpha_m/(Alpha_f*dt*dt*Beta);
 61:   th->scale_F = 1/Alpha_f;
 62:   return(0);
 63: }

 67: static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X)
 68: {
 69:   TS_Alpha       *th = (TS_Alpha*)ts->data;
 70:   Vec            X1 = X,      V1 = th->V1, A1 = th->A1;
 71:   Vec            Xa = th->Xa, Va = th->Va, Aa = th->Aa;
 72:   Vec            X0 = th->X0, V0 = th->V0, A0 = th->A0;
 73:   PetscReal      dt = ts->time_step;
 74:   PetscReal      Alpha_m = th->Alpha_m;
 75:   PetscReal      Alpha_f = th->Alpha_f;
 76:   PetscReal      Gamma   = th->Gamma;
 77:   PetscReal      Beta    = th->Beta;

 81:   /* A1 = ... */
 82:   VecWAXPY(A1,-1.0,X0,X1);
 83:   VecAXPY (A1,-dt,V0);
 84:   VecAXPBY(A1,-(1-2*Beta)/(2*Beta),1/(dt*dt*Beta),A0);
 85:   /* V1 = ... */
 86:   VecWAXPY(V1,(1.0-Gamma)/Gamma,A0,A1);
 87:   VecAYPX (V1,dt*Gamma,V0);
 88:   /* Xa = X0 + Alpha_f*(X1-X0) */
 89:   VecWAXPY(Xa,-1.0,X0,X1);
 90:   VecAYPX (Xa,Alpha_f,X0);
 91:   /* Va = V0 + Alpha_f*(V1-V0) */
 92:   VecWAXPY(Va,-1.0,V0,V1);
 93:   VecAYPX (Va,Alpha_f,V0);
 94:   /* Aa = A0 + Alpha_m*(A1-A0) */
 95:   VecWAXPY(Aa,-1.0,A0,A1);
 96:   VecAYPX (Aa,Alpha_m,A0);
 97:   return(0);
 98: }

102: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
103: {
104:   PetscInt       nits,lits;

108:   SNESSolve(ts->snes,b,x);
109:   SNESGetIterationNumber(ts->snes,&nits);
110:   SNESGetLinearSolveIterations(ts->snes,&lits);
111:   ts->snes_its += nits; ts->ksp_its += lits;
112:   return(0);
113: }

115: /*
116:   Compute a consistent initial state for the generalized-alpha method.
117:   - Solve two successive backward Euler steps with halved time step.
118:   - Compute the initial second time derivative using backward differences.
119:   - If using adaptivity, estimate the LTE of the initial step.
120: */
123: static PetscErrorCode TSAlpha_Restart(TS ts,PetscBool *initok)
124: {
125:   TS_Alpha       *th = (TS_Alpha*)ts->data;
126:   PetscReal      time_step;
127:   PetscReal      alpha_m,alpha_f,gamma,beta;
128:   Vec            X0 = ts->vec_sol, X1, X2 = th->X1;
129:   Vec            V0 = ts->vec_dot, V1, V2 = th->V1;
130:   PetscBool      stageok;

134:   VecDuplicate(X0,&X1);
135:   VecDuplicate(V0,&V1);

137:   /* Setup backward Euler with halved time step */
138:   TSAlpha2GetParams(ts,&alpha_m,&alpha_f,&gamma,&beta);
139:   TSAlpha2SetParams(ts,1,1,1,0.5);
140:   TSGetTimeStep(ts,&time_step);
141:   ts->time_step = time_step/2;
142:   TSAlpha_StageTime(ts);
143:   th->stage_time = ts->ptime;
144:   VecZeroEntries(th->A0);

146:   /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
147:   th->stage_time += ts->time_step;
148:   VecCopy(X0,th->X0);
149:   VecCopy(V0,th->V0);
150:   TSPreStage(ts,th->stage_time);
151:   VecCopy(th->X0,X1);
152:   TS_SNESSolve(ts,NULL,X1);
153:   VecCopy(th->V1,V1);
154:   TSPostStage(ts,th->stage_time,0,&X1);
155:   TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
156:   if (!stageok) goto finally;

158:   /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
159:   th->stage_time += ts->time_step;
160:   VecCopy(X1,th->X0);
161:   VecCopy(V1,th->V0);
162:   TSPreStage(ts,th->stage_time);
163:   VecCopy(th->X0,X2);
164:   TS_SNESSolve(ts,NULL,X2);
165:   VecCopy(th->V1,V2);
166:   TSPostStage(ts,th->stage_time,0,&X2);
167:   TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
168:   if (!stageok) goto finally;

170:   /* Compute A0 ~ dV/dt at t0 with backward differences */
171:   VecZeroEntries(th->A0);
172:   VecAXPY(th->A0,-3/ts->time_step,V0);
173:   VecAXPY(th->A0,+4/ts->time_step,V1);
174:   VecAXPY(th->A0,-1/ts->time_step,V2);

176:   /* Rough, lower-order estimate LTE of the initial step */
177:   if (th->adapt) {
178:     VecZeroEntries(th->vec_lte_work[0]);
179:     VecAXPY(th->vec_lte_work[0],+2,X2);
180:     VecAXPY(th->vec_lte_work[0],-4,X1);
181:     VecAXPY(th->vec_lte_work[0],+2,X0);
182:   }
183:   if (th->adapt) {
184:     VecZeroEntries(th->vec_lte_work[1]);
185:     VecAXPY(th->vec_lte_work[1],+2,V2);
186:     VecAXPY(th->vec_lte_work[1],-4,V1);
187:     VecAXPY(th->vec_lte_work[1],+2,V0);
188:   }

190:  finally:
191:   /* Revert TSAlpha to the initial state (t0,X0,V0) */
192:   if (initok) *initok = stageok;
193:   TSSetTimeStep(ts,time_step);
194:   TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
195:   VecCopy(ts->vec_sol,th->X0);
196:   VecCopy(ts->vec_dot,th->V0);

198:   VecDestroy(&X1);
199:   VecDestroy(&V1);
200:   return(0);
201: }

205: static PetscErrorCode TSStep_Alpha(TS ts)
206: {
207:   TS_Alpha       *th = (TS_Alpha*)ts->data;
208:   PetscInt       rejections = 0;
209:   PetscBool      stageok,accept = PETSC_TRUE;
210:   PetscReal      next_time_step = ts->time_step;

214:   PetscCitationsRegister(citation,&cited);

216:   if (!ts->steprollback) {
217:     if (th->adapt) { VecCopy(th->X0,th->vec_sol_prev); }
218:     if (th->adapt) { VecCopy(th->V0,th->vec_dot_prev); }
219:     VecCopy(ts->vec_sol,th->X0);
220:     VecCopy(ts->vec_dot,th->V0);
221:     VecCopy(th->A1,th->A0);
222:   }

224:   th->status = TS_STEP_INCOMPLETE;
225:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {

227:     if (ts->steprestart) {
228:       TSAlpha_Restart(ts,&stageok);
229:       if (!stageok) goto reject_step;
230:     }

232:     TSAlpha_StageTime(ts);
233:     VecCopy(th->X0,th->X1);
234:     TSPreStage(ts,th->stage_time);
235:     TS_SNESSolve(ts,NULL,th->X1);
236:     TSPostStage(ts,th->stage_time,0,&th->Xa);
237:     TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
238:     if (!stageok) goto reject_step;

240:     th->status = TS_STEP_PENDING;
241:     VecCopy(th->X1,ts->vec_sol);
242:     VecCopy(th->V1,ts->vec_dot);
243:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
244:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
245:     if (!accept) {
246:       VecCopy(th->X0,ts->vec_sol);
247:       VecCopy(th->V0,ts->vec_dot);
248:       ts->time_step = next_time_step;
249:       goto reject_step;
250:     }

252:     ts->ptime += ts->time_step;
253:     ts->time_step = next_time_step;
254:     break;

256:   reject_step:
257:     ts->reject++; accept = PETSC_FALSE;
258:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
259:       ts->reason = TS_DIVERGED_STEP_REJECTED;
260:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
261:     }

263:   }
264:   return(0);
265: }

269: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
270: {
271:   TS_Alpha       *th = (TS_Alpha*)ts->data;
272:   Vec            X = th->X1;              /* X = solution */
273:   Vec            V = th->V1;              /* V = solution */
274:   Vec            Y = th->vec_lte_work[0]; /* Y = X + LTE  */
275:   Vec            Z = th->vec_lte_work[1]; /* Z = V + LTE  */
276:   PetscReal      enormX,enormV;

280:   if (ts->steprestart) {
281:     /* th->vec_{sol|dot}_prev is set to the LTE in TSAlpha_Restart() */
282:     VecAXPY(Y,1,X);
283:     VecAXPY(Z,1,V);
284:   } else {
285:     /* Compute LTE using backward differences with non-constant time step */
286:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
287:     PetscReal   a = 1 + h_prev/h;
288:     PetscScalar scal[3]; Vec vecX[3],vecV[3];
289:     scal[0] = +1/a;   scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
290:     vecX[0] = th->X1; vecX[1] = th->X0;   vecX[2] = th->vec_sol_prev;
291:     vecV[0] = th->V1; vecV[1] = th->V0;   vecV[2] = th->vec_dot_prev;
292:     VecCopy(X,Y);
293:     VecMAXPY(Y,3,scal,vecX);
294:     VecCopy(V,Z);
295:     VecMAXPY(Z,3,scal,vecV);
296:   }
297:   /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
298:   TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX);
299:   TSErrorWeightedNorm(ts,V,Z,wnormtype,&enormV);
300:   if (wnormtype == NORM_2)
301:     *wlte = PetscSqrtReal(PetscSqr(enormX)/2 + PetscSqr(enormV)/2);
302:   else
303:     *wlte = PetscMax(enormX,enormV);
304:   if (order) *order = 2;
305:   return(0);
306: }

310: static PetscErrorCode TSRollBack_Alpha(TS ts)
311: {
312:   TS_Alpha       *th = (TS_Alpha*)ts->data;

316:   VecCopy(th->X0,ts->vec_sol);
317:   VecCopy(th->V0,ts->vec_dot);
318:   return(0);
319: }

321: /*
324: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
325: {
326:   TS_Alpha       *th = (TS_Alpha*)ts->data;
327:   PetscReal      dt  = t - ts->ptime;

331:   VecCopy(ts->vec_dot,V);
332:   VecAXPY(V,dt*(1-th->Gamma),th->A0);
333:   VecAXPY(V,dt*th->Gamma,th->A1);
334:   VecCopy(ts->vec_sol,X);
335:   VecAXPY(X,dt,V);
336:   VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
337:   VecAXPY(X,dt*dt*th->Beta,th->A1);
338:   return(0);
339: }
340: */

344: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
345: {
346:   TS_Alpha       *th = (TS_Alpha*)ts->data;
347:   PetscReal      ta = th->stage_time;
348:   Vec            Xa = th->Xa, Va = th->Va, Aa = th->Aa;

352:   TSAlpha_StageVecs(ts,X);
353:   /* F = Function(ta,Xa,Va,Aa) */
354:   TSComputeI2Function(ts,ta,Xa,Va,Aa,F);
355:   VecScale(F,th->scale_F);
356:   return(0);
357: }

361: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)
362: {
363:   TS_Alpha       *th = (TS_Alpha*)ts->data;
364:   PetscReal      ta = th->stage_time;
365:   Vec            Xa = th->Xa, Va = th->Va, Aa = th->Aa;
366:   PetscReal      dVdX = th->shift_V, dAdX = th->shift_A;

370:   /* J,P = Jacobian(ta,Xa,Va,Aa) */
371:   TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P);
372:   return(0);
373: }

377: static PetscErrorCode TSReset_Alpha(TS ts)
378: {
379:   TS_Alpha       *th = (TS_Alpha*)ts->data;

383:   VecDestroy(&th->X0);
384:   VecDestroy(&th->Xa);
385:   VecDestroy(&th->X1);
386:   VecDestroy(&th->V0);
387:   VecDestroy(&th->Va);
388:   VecDestroy(&th->V1);
389:   VecDestroy(&th->A0);
390:   VecDestroy(&th->Aa);
391:   VecDestroy(&th->A1);
392:   VecDestroy(&th->vec_sol_prev);
393:   VecDestroy(&th->vec_dot_prev);
394:   VecDestroy(&th->vec_lte_work[0]);
395:   VecDestroy(&th->vec_lte_work[1]);
396:   return(0);
397: }

401: static PetscErrorCode TSDestroy_Alpha(TS ts)
402: {

406:   TSReset_Alpha(ts);
407:   PetscFree(ts->data);

409:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2UseAdapt_C",NULL);
410:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL);
411:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL);
412:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL);
413:   return(0);
414: }

418: static PetscErrorCode TSSetUp_Alpha(TS ts)
419: {
420:   TS_Alpha       *th = (TS_Alpha*)ts->data;

424:   VecDuplicate(ts->vec_sol,&th->X0);
425:   VecDuplicate(ts->vec_sol,&th->Xa);
426:   VecDuplicate(ts->vec_sol,&th->X1);
427:   VecDuplicate(ts->vec_sol,&th->V0);
428:   VecDuplicate(ts->vec_sol,&th->Va);
429:   VecDuplicate(ts->vec_sol,&th->V1);
430:   VecDuplicate(ts->vec_sol,&th->A0);
431:   VecDuplicate(ts->vec_sol,&th->Aa);
432:   VecDuplicate(ts->vec_sol,&th->A1);

434:   TSGetAdapt(ts,&ts->adapt);
435:   TSAdaptCandidatesClear(ts->adapt);
436:   if (!th->adapt) {
437:     TSAdaptSetType(ts->adapt,TSADAPTNONE);
438:   } else {
439:     VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
440:     VecDuplicate(ts->vec_sol,&th->vec_dot_prev);
441:     VecDuplicate(ts->vec_sol,&th->vec_lte_work[0]);
442:     VecDuplicate(ts->vec_sol,&th->vec_lte_work[1]);
443:     if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
444:       ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
445:   }

447:   TSGetSNES(ts,&ts->snes);
448:   return(0);
449: }

453: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
454: {
455:   TS_Alpha       *th = (TS_Alpha*)ts->data;

459:   PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
460:   {
461:     PetscBool flg;
462:     PetscReal radius = 1;
463:     PetscBool adapt  = th->adapt;
464:     PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg);
465:     if (flg) {TSAlpha2SetRadius(ts,radius);}
466:     PetscOptionsReal("-ts_alpha_alpha_m","Algoritmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL);
467:     PetscOptionsReal("-ts_alpha_alpha_f","Algoritmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL);
468:     PetscOptionsReal("-ts_alpha_gamma","Algoritmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL);
469:     PetscOptionsReal("-ts_alpha_beta","Algoritmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL);
470:     TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta);
471:     PetscOptionsBool("-ts_alpha_adapt","Use time-step adaptivity with the Alpha method","TSAlpha2UseAdapt",adapt,&adapt,&flg);
472:     if (flg) {TSAlpha2UseAdapt(ts,adapt);}
473:   }
474:   PetscOptionsTail();
475:   return(0);
476: }

480: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
481: {
482:   TS_Alpha       *th = (TS_Alpha*)ts->data;
483:   PetscBool      iascii;

487:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
488:   if (iascii) {
489:     PetscViewerASCIIPrintf(viewer,"  Alpha_m=%g, Alpha_f=%g, Gamma=%g, Beta=%g\n",(double)th->Alpha_m,(double)th->Alpha_f,(double)th->Gamma,(double)th->Beta);
490:   }
491:   if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
492:   if (ts->snes)  {SNESView(ts->snes,viewer);}
493:   return(0);
494: }

498: static PetscErrorCode TSAlpha2UseAdapt_Alpha(TS ts,PetscBool use)
499: {
500:   TS_Alpha *th = (TS_Alpha*)ts->data;

503:   if (use == th->adapt) return(0);
504:   if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
505:   th->adapt = use;
506:   return(0);
507: }

511: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)
512: {
513:   PetscReal      alpha_m,alpha_f,gamma,beta;

517:   if (radius < 0 || radius > 1) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
518:   alpha_m = (2-radius)/(1+radius);
519:   alpha_f = 1/(1+radius);
520:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
521:   beta    = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta;
522:   TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
523:   return(0);
524: }

528: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
529: {
530:   TS_Alpha  *th = (TS_Alpha*)ts->data;
531:   PetscReal tol = 100*PETSC_MACHINE_EPSILON;
532:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

535:   th->Alpha_m = alpha_m;
536:   th->Alpha_f = alpha_f;
537:   th->Gamma   = gamma;
538:   th->Beta    = beta;
539:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
540:   return(0);
541: }

545: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
546: {
547:   TS_Alpha *th = (TS_Alpha*)ts->data;

550:   if (alpha_m) *alpha_m = th->Alpha_m;
551:   if (alpha_f) *alpha_f = th->Alpha_f;
552:   if (gamma)   *gamma   = th->Gamma;
553:   if (beta)    *beta    = th->Beta;
554:   return(0);
555: }

557: /*MC
558:       TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
559:                  for second-order systems

561:   Level: beginner

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

568: .seealso:  TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams()
569: M*/
572: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
573: {
574:   TS_Alpha       *th;

578:   ts->ops->reset          = TSReset_Alpha;
579:   ts->ops->destroy        = TSDestroy_Alpha;
580:   ts->ops->view           = TSView_Alpha;
581:   ts->ops->setup          = TSSetUp_Alpha;
582:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
583:   ts->ops->step           = TSStep_Alpha;
584:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
585:   ts->ops->rollback       = TSRollBack_Alpha;
586:   /*ts->ops->interpolate  = TSInterpolate_Alpha;*/
587:   ts->ops->snesfunction   = SNESTSFormFunction_Alpha;
588:   ts->ops->snesjacobian   = SNESTSFormJacobian_Alpha;

590:   PetscNewLog(ts,&th);
591:   ts->data = (void*)th;

593:   th->Alpha_m = 0.5;
594:   th->Alpha_f = 0.5;
595:   th->Gamma   = 0.5;
596:   th->Beta    = 0.25;
597:   th->order   = 2;

599:   th->adapt = PETSC_FALSE;

601:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2UseAdapt_C",TSAlpha2UseAdapt_Alpha);
602:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha);
603:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha);
604:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha);
605:   return(0);
606: }

610: /*@
611:   TSAlpha2UseAdapt - Use time-step adaptivity with the Alpha method

613:   Logically Collective on TS

615:   Input Parameter:
616: +  ts - timestepping context
617: -  use - flag to use adaptivity

619:   Options Database:
620: .  -ts_alpha_adapt

622:   Level: intermediate

624: .seealso: TSAdapt, TSADAPTBASIC
625: @*/
626: PetscErrorCode TSAlpha2UseAdapt(TS ts,PetscBool use)
627: {

633:   PetscTryMethod(ts,"TSAlpha2UseAdapt_C",(TS,PetscBool),(ts,use));
634:   return(0);
635: }

639: /*@
640:   TSAlpha2SetRadius - sets the desired spectral radius of the method
641:                       (i.e. high-frequency numerical damping)

643:   Logically Collective on TS

645:   The algorithmic parameters \alpha_m and \alpha_f of the
646:   generalized-\alpha method can be computed in terms of a specified
647:   spectral radius \rho in [0,1] for infinite time step in order to
648:   control high-frequency numerical damping:
649:     \alpha_m = (2-\rho)/(1+\rho)
650:     \alpha_f = 1/(1+\rho)

652:   Input Parameter:
653: +  ts - timestepping context
654: -  radius - the desired spectral radius

656:   Options Database:
657: .  -ts_alpha_radius <radius>

659:   Level: intermediate

661: .seealso: TSAlpha2SetParams(), TSAlpha2GetParams()
662: @*/
663: PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius)
664: {

670:   if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
671:   PetscTryMethod(ts,"TSAlpha2SetRadius_C",(TS,PetscReal),(ts,radius));
672:   return(0);
673: }

677: /*@
678:   TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2

680:   Logically Collective on TS

682:   Second-order accuracy can be obtained so long as:
683:     \gamma = 1/2 + alpha_m - alpha_f
684:     \beta  = 1/4 (1 + alpha_m - alpha_f)^2

686:   Unconditional stability requires:
687:     \alpha_m >= \alpha_f >= 1/2


690:   Input Parameter:
691: + ts - timestepping context
692: . \alpha_m - algorithmic paramenter
693: . \alpha_f - algorithmic paramenter
694: . \gamma   - algorithmic paramenter
695: - \beta    - algorithmic paramenter

697:   Options Database:
698: + -ts_alpha_alpha_m <alpha_m>
699: . -ts_alpha_alpha_f <alpha_f>
700: . -ts_alpha_gamma   <gamma>
701: - -ts_alpha_beta    <beta>

703:   Note:
704:   Use of this function is normally only required to hack TSALPHA2 to
705:   use a modified integration scheme. Users should call
706:   TSAlpha2SetRadius() to set the desired spectral radius of the methods
707:   (i.e. high-frequency damping) in order so select optimal values for
708:   these parameters.

710:   Level: advanced

712: .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams()
713: @*/
714: PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
715: {

724:   PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta));
725:   return(0);
726: }

730: /*@
731:   TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2

733:   Not Collective

735:   Input Parameter:
736: . ts - timestepping context

738:   Output Parameters:
739: + \alpha_m - algorithmic parameter
740: . \alpha_f - algorithmic parameter
741: . \gamma   - algorithmic parameter
742: - \beta    - algorithmic parameter

744:   Note:
745:   Use of this function is normally only required to hack TSALPHA2 to
746:   use a modified integration scheme. Users should call
747:   TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral
748:   radius of the method) in order so select optimal values for these
749:   parameters.

751:   Level: advanced

753: .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams()
754: @*/
755: PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
756: {

765:   PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta));
766:   return(0);
767: }