Actual source code: alpha2.c

petsc-3.14.6 2021-03-30
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>

  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:   Vec       vec_sol_prev;
 38:   Vec       vec_dot_prev;
 39:   Vec       vec_lte_work[2];

 41:   TSStepStatus status;
 42: } TS_Alpha;

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

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

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

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

 95: static PetscErrorCode TSAlpha_SNESSolve(TS ts,Vec b,Vec x)
 96: {
 97:   PetscInt       nits,lits;

101:   SNESSolve(ts->snes,b,x);
102:   SNESGetIterationNumber(ts->snes,&nits);
103:   SNESGetLinearSolveIterations(ts->snes,&lits);
104:   ts->snes_its += nits; ts->ksp_its += lits;
105:   return(0);
106: }

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

125:   VecDuplicate(X0,&X1);
126:   VecDuplicate(V0,&V1);

128:   /* Setup backward Euler with halved time step */
129:   TSAlpha2GetParams(ts,&alpha_m,&alpha_f,&gamma,&beta);
130:   TSAlpha2SetParams(ts,1,1,1,0.5);
131:   TSGetTimeStep(ts,&time_step);
132:   ts->time_step = time_step/2;
133:   TSAlpha_StageTime(ts);
134:   th->stage_time = ts->ptime;
135:   VecZeroEntries(th->A0);

137:   /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
138:   th->stage_time += ts->time_step;
139:   VecCopy(X0,th->X0);
140:   VecCopy(V0,th->V0);
141:   TSPreStage(ts,th->stage_time);
142:   VecCopy(th->X0,X1);
143:   TSAlpha_SNESSolve(ts,NULL,X1);
144:   VecCopy(th->V1,V1);
145:   TSPostStage(ts,th->stage_time,0,&X1);
146:   TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
147:   if (!stageok) goto finally;

149:   /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
150:   th->stage_time += ts->time_step;
151:   VecCopy(X1,th->X0);
152:   VecCopy(V1,th->V0);
153:   TSPreStage(ts,th->stage_time);
154:   VecCopy(th->X0,X2);
155:   TSAlpha_SNESSolve(ts,NULL,X2);
156:   VecCopy(th->V1,V2);
157:   TSPostStage(ts,th->stage_time,0,&X2);
158:   TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
159:   if (!stageok) goto finally;

161:   /* Compute A0 ~ dV/dt at t0 with backward differences */
162:   VecZeroEntries(th->A0);
163:   VecAXPY(th->A0,-3/ts->time_step,V0);
164:   VecAXPY(th->A0,+4/ts->time_step,V1);
165:   VecAXPY(th->A0,-1/ts->time_step,V2);

167:   /* Rough, lower-order estimate LTE of the initial step */
168:   if (th->vec_lte_work[0]) {
169:     VecZeroEntries(th->vec_lte_work[0]);
170:     VecAXPY(th->vec_lte_work[0],+2,X2);
171:     VecAXPY(th->vec_lte_work[0],-4,X1);
172:     VecAXPY(th->vec_lte_work[0],+2,X0);
173:   }
174:   if (th->vec_lte_work[1]) {
175:     VecZeroEntries(th->vec_lte_work[1]);
176:     VecAXPY(th->vec_lte_work[1],+2,V2);
177:     VecAXPY(th->vec_lte_work[1],-4,V1);
178:     VecAXPY(th->vec_lte_work[1],+2,V0);
179:   }

181:  finally:
182:   /* Revert TSAlpha to the initial state (t0,X0,V0) */
183:   if (initok) *initok = stageok;
184:   TSSetTimeStep(ts,time_step);
185:   TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
186:   VecCopy(ts->vec_sol,th->X0);
187:   VecCopy(ts->vec_dot,th->V0);

189:   VecDestroy(&X1);
190:   VecDestroy(&V1);
191:   return(0);
192: }

194: static PetscErrorCode TSStep_Alpha(TS ts)
195: {
196:   TS_Alpha       *th = (TS_Alpha*)ts->data;
197:   PetscInt       rejections = 0;
198:   PetscBool      stageok,accept = PETSC_TRUE;
199:   PetscReal      next_time_step = ts->time_step;

203:   PetscCitationsRegister(citation,&cited);

205:   if (!ts->steprollback) {
206:     if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
207:     if (th->vec_dot_prev) { VecCopy(th->V0,th->vec_dot_prev); }
208:     VecCopy(ts->vec_sol,th->X0);
209:     VecCopy(ts->vec_dot,th->V0);
210:     VecCopy(th->A1,th->A0);
211:   }

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

216:     if (ts->steprestart) {
217:       TSAlpha_Restart(ts,&stageok);
218:       if (!stageok) goto reject_step;
219:     }

221:     TSAlpha_StageTime(ts);
222:     VecCopy(th->X0,th->X1);
223:     TSPreStage(ts,th->stage_time);
224:     TSAlpha_SNESSolve(ts,NULL,th->X1);
225:     TSPostStage(ts,th->stage_time,0,&th->Xa);
226:     TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
227:     if (!stageok) goto reject_step;

229:     th->status = TS_STEP_PENDING;
230:     VecCopy(th->X1,ts->vec_sol);
231:     VecCopy(th->V1,ts->vec_dot);
232:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
233:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
234:     if (!accept) {
235:       VecCopy(th->X0,ts->vec_sol);
236:       VecCopy(th->V0,ts->vec_dot);
237:       ts->time_step = next_time_step;
238:       goto reject_step;
239:     }

241:     ts->ptime += ts->time_step;
242:     ts->time_step = next_time_step;
243:     break;

245:   reject_step:
246:     ts->reject++; accept = PETSC_FALSE;
247:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
248:       ts->reason = TS_DIVERGED_STEP_REJECTED;
249:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
250:     }

252:   }
253:   return(0);
254: }

256: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
257: {
258:   TS_Alpha       *th = (TS_Alpha*)ts->data;
259:   Vec            X = th->X1;              /* X = solution */
260:   Vec            V = th->V1;              /* V = solution */
261:   Vec            Y = th->vec_lte_work[0]; /* Y = X + LTE  */
262:   Vec            Z = th->vec_lte_work[1]; /* Z = V + LTE  */
263:   PetscReal      enormX,enormV,enormXa,enormVa,enormXr,enormVr;

267:   if (!th->vec_sol_prev) {*wlte = -1; return(0);}
268:   if (!th->vec_dot_prev) {*wlte = -1; return(0);}
269:   if (!th->vec_lte_work[0]) {*wlte = -1; return(0);}
270:   if (!th->vec_lte_work[1]) {*wlte = -1; return(0);}
271:   if (ts->steprestart) {
272:     /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
273:     VecAXPY(Y,1,X);
274:     VecAXPY(Z,1,V);
275:   } else {
276:     /* Compute LTE using backward differences with non-constant time step */
277:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
278:     PetscReal   a = 1 + h_prev/h;
279:     PetscScalar scal[3]; Vec vecX[3],vecV[3];
280:     scal[0] = +1/a;   scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
281:     vecX[0] = th->X1; vecX[1] = th->X0;   vecX[2] = th->vec_sol_prev;
282:     vecV[0] = th->V1; vecV[1] = th->V0;   vecV[2] = th->vec_dot_prev;
283:     VecCopy(X,Y);
284:     VecMAXPY(Y,3,scal,vecX);
285:     VecCopy(V,Z);
286:     VecMAXPY(Z,3,scal,vecV);
287:   }
288:   /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
289:   TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX,&enormXa,&enormXr);
290:   TSErrorWeightedNorm(ts,V,Z,wnormtype,&enormV,&enormVa,&enormVr);
291:   if (wnormtype == NORM_2)
292:     *wlte = PetscSqrtReal(PetscSqr(enormX)/2 + PetscSqr(enormV)/2);
293:   else
294:     *wlte = PetscMax(enormX,enormV);
295:   if (order) *order = 2;
296:   return(0);
297: }

299: static PetscErrorCode TSRollBack_Alpha(TS ts)
300: {
301:   TS_Alpha       *th = (TS_Alpha*)ts->data;

305:   VecCopy(th->X0,ts->vec_sol);
306:   VecCopy(th->V0,ts->vec_dot);
307:   return(0);
308: }

310: /*
311: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
312: {
313:   TS_Alpha       *th = (TS_Alpha*)ts->data;
314:   PetscReal      dt  = t - ts->ptime;

318:   VecCopy(ts->vec_dot,V);
319:   VecAXPY(V,dt*(1-th->Gamma),th->A0);
320:   VecAXPY(V,dt*th->Gamma,th->A1);
321:   VecCopy(ts->vec_sol,X);
322:   VecAXPY(X,dt,V);
323:   VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
324:   VecAXPY(X,dt*dt*th->Beta,th->A1);
325:   return(0);
326: }
327: */

329: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
330: {
331:   TS_Alpha       *th = (TS_Alpha*)ts->data;
332:   PetscReal      ta = th->stage_time;
333:   Vec            Xa = th->Xa, Va = th->Va, Aa = th->Aa;

337:   TSAlpha_StageVecs(ts,X);
338:   /* F = Function(ta,Xa,Va,Aa) */
339:   TSComputeI2Function(ts,ta,Xa,Va,Aa,F);
340:   VecScale(F,th->scale_F);
341:   return(0);
342: }

344: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,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;
349:   PetscReal      dVdX = th->shift_V, dAdX = th->shift_A;

353:   /* J,P = Jacobian(ta,Xa,Va,Aa) */
354:   TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P);
355:   return(0);
356: }

358: static PetscErrorCode TSReset_Alpha(TS ts)
359: {
360:   TS_Alpha       *th = (TS_Alpha*)ts->data;

364:   VecDestroy(&th->X0);
365:   VecDestroy(&th->Xa);
366:   VecDestroy(&th->X1);
367:   VecDestroy(&th->V0);
368:   VecDestroy(&th->Va);
369:   VecDestroy(&th->V1);
370:   VecDestroy(&th->A0);
371:   VecDestroy(&th->Aa);
372:   VecDestroy(&th->A1);
373:   VecDestroy(&th->vec_sol_prev);
374:   VecDestroy(&th->vec_dot_prev);
375:   VecDestroy(&th->vec_lte_work[0]);
376:   VecDestroy(&th->vec_lte_work[1]);
377:   return(0);
378: }

380: static PetscErrorCode TSDestroy_Alpha(TS ts)
381: {

385:   TSReset_Alpha(ts);
386:   PetscFree(ts->data);

388:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL);
389:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL);
390:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL);
391:   return(0);
392: }

394: static PetscErrorCode TSSetUp_Alpha(TS ts)
395: {
396:   TS_Alpha       *th = (TS_Alpha*)ts->data;
397:   PetscBool      match;

401:   VecDuplicate(ts->vec_sol,&th->X0);
402:   VecDuplicate(ts->vec_sol,&th->Xa);
403:   VecDuplicate(ts->vec_sol,&th->X1);
404:   VecDuplicate(ts->vec_sol,&th->V0);
405:   VecDuplicate(ts->vec_sol,&th->Va);
406:   VecDuplicate(ts->vec_sol,&th->V1);
407:   VecDuplicate(ts->vec_sol,&th->A0);
408:   VecDuplicate(ts->vec_sol,&th->Aa);
409:   VecDuplicate(ts->vec_sol,&th->A1);

411:   TSGetAdapt(ts,&ts->adapt);
412:   TSAdaptCandidatesClear(ts->adapt);
413:   PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
414:   if (!match) {
415:     VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
416:     VecDuplicate(ts->vec_sol,&th->vec_dot_prev);
417:     VecDuplicate(ts->vec_sol,&th->vec_lte_work[0]);
418:     VecDuplicate(ts->vec_sol,&th->vec_lte_work[1]);
419:   }

421:   TSGetSNES(ts,&ts->snes);
422:   return(0);
423: }

425: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
426: {
427:   TS_Alpha       *th = (TS_Alpha*)ts->data;

431:   PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
432:   {
433:     PetscBool flg;
434:     PetscReal radius = 1;
435:     PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg);
436:     if (flg) {TSAlpha2SetRadius(ts,radius);}
437:     PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL);
438:     PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL);
439:     PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL);
440:     PetscOptionsReal("-ts_alpha_beta","Algorithmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL);
441:     TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta);
442:   }
443:   PetscOptionsTail();
444:   return(0);
445: }

447: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
448: {
449:   TS_Alpha       *th = (TS_Alpha*)ts->data;
450:   PetscBool      iascii;

454:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
455:   if (iascii) {
456:     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);
457:   }
458:   return(0);
459: }

461: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)
462: {
463:   PetscReal      alpha_m,alpha_f,gamma,beta;

467:   if (radius < 0 || radius > 1) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
468:   alpha_m = (2-radius)/(1+radius);
469:   alpha_f = 1/(1+radius);
470:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
471:   beta    = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta;
472:   TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
473:   return(0);
474: }

476: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
477: {
478:   TS_Alpha  *th = (TS_Alpha*)ts->data;
479:   PetscReal tol = 100*PETSC_MACHINE_EPSILON;
480:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

483:   th->Alpha_m = alpha_m;
484:   th->Alpha_f = alpha_f;
485:   th->Gamma   = gamma;
486:   th->Beta    = beta;
487:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
488:   return(0);
489: }

491: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
492: {
493:   TS_Alpha *th = (TS_Alpha*)ts->data;

496:   if (alpha_m) *alpha_m = th->Alpha_m;
497:   if (alpha_f) *alpha_f = th->Alpha_f;
498:   if (gamma)   *gamma   = th->Gamma;
499:   if (beta)    *beta    = th->Beta;
500:   return(0);
501: }

503: /*MC
504:       TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
505:                  for second-order systems

507:   Level: beginner

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

514: .seealso:  TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams()
515: M*/
516: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
517: {
518:   TS_Alpha       *th;

522:   ts->ops->reset          = TSReset_Alpha;
523:   ts->ops->destroy        = TSDestroy_Alpha;
524:   ts->ops->view           = TSView_Alpha;
525:   ts->ops->setup          = TSSetUp_Alpha;
526:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
527:   ts->ops->step           = TSStep_Alpha;
528:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
529:   ts->ops->rollback       = TSRollBack_Alpha;
530:   /*ts->ops->interpolate  = TSInterpolate_Alpha;*/
531:   ts->ops->snesfunction   = SNESTSFormFunction_Alpha;
532:   ts->ops->snesjacobian   = SNESTSFormJacobian_Alpha;
533:   ts->default_adapt_type  = TSADAPTNONE;

535:   ts->usessnes = PETSC_TRUE;

537:   PetscNewLog(ts,&th);
538:   ts->data = (void*)th;

540:   th->Alpha_m = 0.5;
541:   th->Alpha_f = 0.5;
542:   th->Gamma   = 0.5;
543:   th->Beta    = 0.25;
544:   th->order   = 2;

546:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha);
547:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha);
548:   PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha);
549:   return(0);
550: }

552: /*@
553:   TSAlpha2SetRadius - sets the desired spectral radius of the method
554:                       (i.e. high-frequency numerical damping)

556:   Logically Collective on TS

558:   The algorithmic parameters \alpha_m and \alpha_f of the
559:   generalized-\alpha method can be computed in terms of a specified
560:   spectral radius \rho in [0,1] for infinite time step in order to
561:   control high-frequency numerical damping:
562:     \alpha_m = (2-\rho)/(1+\rho)
563:     \alpha_f = 1/(1+\rho)

565:   Input Parameter:
566: +  ts - timestepping context
567: -  radius - the desired spectral radius

569:   Options Database:
570: .  -ts_alpha_radius <radius>

572:   Level: intermediate

574: .seealso: TSAlpha2SetParams(), TSAlpha2GetParams()
575: @*/
576: PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius)
577: {

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

588: /*@
589:   TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2

591:   Logically Collective on TS

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

597:   Unconditional stability requires:
598:     \alpha_m >= \alpha_f >= 1/2


601:   Input Parameter:
602: + ts - timestepping context
603: . \alpha_m - algorithmic paramenter
604: . \alpha_f - algorithmic paramenter
605: . \gamma   - algorithmic paramenter
606: - \beta    - algorithmic paramenter

608:   Options Database:
609: + -ts_alpha_alpha_m <alpha_m>
610: . -ts_alpha_alpha_f <alpha_f>
611: . -ts_alpha_gamma   <gamma>
612: - -ts_alpha_beta    <beta>

614:   Note:
615:   Use of this function is normally only required to hack TSALPHA2 to
616:   use a modified integration scheme. Users should call
617:   TSAlpha2SetRadius() to set the desired spectral radius of the methods
618:   (i.e. high-frequency damping) in order so select optimal values for
619:   these parameters.

621:   Level: advanced

623: .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams()
624: @*/
625: PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
626: {

635:   PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta));
636:   return(0);
637: }

639: /*@
640:   TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2

642:   Not Collective

644:   Input Parameter:
645: . ts - timestepping context

647:   Output Parameters:
648: + \alpha_m - algorithmic parameter
649: . \alpha_f - algorithmic parameter
650: . \gamma   - algorithmic parameter
651: - \beta    - algorithmic parameter

653:   Note:
654:   Use of this function is normally only required to hack TSALPHA2 to
655:   use a modified integration scheme. Users should call
656:   TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral
657:   radius of the method) in order so select optimal values for these
658:   parameters.

660:   Level: advanced

662: .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams()
663: @*/
664: PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
665: {

674:   PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta));
675:   return(0);
676: }