Actual source code: alpha2.c

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

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

 28:   Vec vec_dot;

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

 36:   Vec vec_sol_prev;
 37:   Vec vec_dot_prev;
 38:   Vec vec_lte_work[2];

 40:   TSStepStatus status;
 41: } TS_Alpha;

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

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

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

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

 91: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
 92: {
 93:   PetscInt nits, lits;

 95:   SNESSolve(ts->snes, b, x);
 96:   SNESGetIterationNumber(ts->snes, &nits);
 97:   SNESGetLinearSolveIterations(ts->snes, &lits);
 98:   ts->snes_its += nits;
 99:   ts->ksp_its += lits;
100:   return 0;
101: }

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

118:   VecDuplicate(X0, &X1);
119:   VecDuplicate(V0, &V1);

121:   /* Setup backward Euler with halved time step */
122:   TSAlpha2GetParams(ts, &alpha_m, &alpha_f, &gamma, &beta);
123:   TSAlpha2SetParams(ts, 1, 1, 1, 0.5);
124:   TSGetTimeStep(ts, &time_step);
125:   ts->time_step = time_step / 2;
126:   TSAlpha_StageTime(ts);
127:   th->stage_time = ts->ptime;
128:   VecZeroEntries(th->A0);

130:   /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
131:   th->stage_time += ts->time_step;
132:   VecCopy(X0, th->X0);
133:   VecCopy(V0, th->V0);
134:   TSPreStage(ts, th->stage_time);
135:   VecCopy(th->X0, X1);
136:   TSAlpha_SNESSolve(ts, NULL, X1);
137:   VecCopy(th->V1, V1);
138:   TSPostStage(ts, th->stage_time, 0, &X1);
139:   TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok);
140:   if (!stageok) goto finally;

142:   /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
143:   th->stage_time += ts->time_step;
144:   VecCopy(X1, th->X0);
145:   VecCopy(V1, th->V0);
146:   TSPreStage(ts, th->stage_time);
147:   VecCopy(th->X0, X2);
148:   TSAlpha_SNESSolve(ts, NULL, X2);
149:   VecCopy(th->V1, V2);
150:   TSPostStage(ts, th->stage_time, 0, &X2);
151:   TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok);
152:   if (!stageok) goto finally;

154:   /* Compute A0 ~ dV/dt at t0 with backward differences */
155:   VecZeroEntries(th->A0);
156:   VecAXPY(th->A0, -3 / ts->time_step, V0);
157:   VecAXPY(th->A0, +4 / ts->time_step, V1);
158:   VecAXPY(th->A0, -1 / ts->time_step, V2);

160:   /* Rough, lower-order estimate LTE of the initial step */
161:   if (th->vec_lte_work[0]) {
162:     VecZeroEntries(th->vec_lte_work[0]);
163:     VecAXPY(th->vec_lte_work[0], +2, X2);
164:     VecAXPY(th->vec_lte_work[0], -4, X1);
165:     VecAXPY(th->vec_lte_work[0], +2, X0);
166:   }
167:   if (th->vec_lte_work[1]) {
168:     VecZeroEntries(th->vec_lte_work[1]);
169:     VecAXPY(th->vec_lte_work[1], +2, V2);
170:     VecAXPY(th->vec_lte_work[1], -4, V1);
171:     VecAXPY(th->vec_lte_work[1], +2, V0);
172:   }

174: finally:
175:   /* Revert TSAlpha to the initial state (t0,X0,V0) */
176:   if (initok) *initok = stageok;
177:   TSSetTimeStep(ts, time_step);
178:   TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta);
179:   VecCopy(ts->vec_sol, th->X0);
180:   VecCopy(ts->vec_dot, th->V0);

182:   VecDestroy(&X1);
183:   VecDestroy(&V1);
184:   return 0;
185: }

187: static PetscErrorCode TSStep_Alpha(TS ts)
188: {
189:   TS_Alpha *th         = (TS_Alpha *)ts->data;
190:   PetscInt  rejections = 0;
191:   PetscBool stageok, accept = PETSC_TRUE;
192:   PetscReal next_time_step = ts->time_step;

194:   PetscCitationsRegister(citation, &cited);

196:   if (!ts->steprollback) {
197:     if (th->vec_sol_prev) VecCopy(th->X0, th->vec_sol_prev);
198:     if (th->vec_dot_prev) VecCopy(th->V0, th->vec_dot_prev);
199:     VecCopy(ts->vec_sol, th->X0);
200:     VecCopy(ts->vec_dot, th->V0);
201:     VecCopy(th->A1, th->A0);
202:   }

204:   th->status = TS_STEP_INCOMPLETE;
205:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
206:     if (ts->steprestart) {
207:       TSAlpha_Restart(ts, &stageok);
208:       if (!stageok) goto reject_step;
209:     }

211:     TSAlpha_StageTime(ts);
212:     VecCopy(th->X0, th->X1);
213:     TSPreStage(ts, th->stage_time);
214:     TSAlpha_SNESSolve(ts, NULL, th->X1);
215:     TSPostStage(ts, th->stage_time, 0, &th->Xa);
216:     TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok);
217:     if (!stageok) goto reject_step;

219:     th->status = TS_STEP_PENDING;
220:     VecCopy(th->X1, ts->vec_sol);
221:     VecCopy(th->V1, ts->vec_dot);
222:     TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
223:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
224:     if (!accept) {
225:       VecCopy(th->X0, ts->vec_sol);
226:       VecCopy(th->V0, ts->vec_dot);
227:       ts->time_step = next_time_step;
228:       goto reject_step;
229:     }

231:     ts->ptime += ts->time_step;
232:     ts->time_step = next_time_step;
233:     break;

235:   reject_step:
236:     ts->reject++;
237:     accept = PETSC_FALSE;
238:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
239:       ts->reason = TS_DIVERGED_STEP_REJECTED;
240:       PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
241:     }
242:   }
243:   return 0;
244: }

246: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
247: {
248:   TS_Alpha *th = (TS_Alpha *)ts->data;
249:   Vec       X  = th->X1;              /* X = solution */
250:   Vec       V  = th->V1;              /* V = solution */
251:   Vec       Y  = th->vec_lte_work[0]; /* Y = X + LTE  */
252:   Vec       Z  = th->vec_lte_work[1]; /* Z = V + LTE  */
253:   PetscReal enormX, enormV, enormXa, enormVa, enormXr, enormVr;

255:   if (!th->vec_sol_prev) {
256:     *wlte = -1;
257:     return 0;
258:   }
259:   if (!th->vec_dot_prev) {
260:     *wlte = -1;
261:     return 0;
262:   }
263:   if (!th->vec_lte_work[0]) {
264:     *wlte = -1;
265:     return 0;
266:   }
267:   if (!th->vec_lte_work[1]) {
268:     *wlte = -1;
269:     return 0;
270:   }
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];
280:     Vec         vecX[3], vecV[3];
281:     scal[0] = +1 / a;
282:     scal[1] = -1 / (a - 1);
283:     scal[2] = +1 / (a * (a - 1));
284:     vecX[0] = th->X1;
285:     vecX[1] = th->X0;
286:     vecX[2] = th->vec_sol_prev;
287:     vecV[0] = th->V1;
288:     vecV[1] = th->V0;
289:     vecV[2] = th->vec_dot_prev;
290:     VecCopy(X, Y);
291:     VecMAXPY(Y, 3, scal, vecX);
292:     VecCopy(V, Z);
293:     VecMAXPY(Z, 3, scal, vecV);
294:   }
295:   /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
296:   TSErrorWeightedNorm(ts, X, Y, wnormtype, &enormX, &enormXa, &enormXr);
297:   TSErrorWeightedNorm(ts, V, Z, wnormtype, &enormV, &enormVa, &enormVr);
298:   if (wnormtype == NORM_2) *wlte = PetscSqrtReal(PetscSqr(enormX) / 2 + PetscSqr(enormV) / 2);
299:   else *wlte = PetscMax(enormX, enormV);
300:   if (order) *order = 2;
301:   return 0;
302: }

304: static PetscErrorCode TSRollBack_Alpha(TS ts)
305: {
306:   TS_Alpha *th = (TS_Alpha *)ts->data;

308:   VecCopy(th->X0, ts->vec_sol);
309:   VecCopy(th->V0, ts->vec_dot);
310:   return 0;
311: }

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

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

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

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

343: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
344: {
345:   TS_Alpha *th = (TS_Alpha *)ts->data;
346:   PetscReal ta = th->stage_time;
347:   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
348:   PetscReal dVdX = th->shift_V, dAdX = th->shift_A;

350:   /* J,P = Jacobian(ta,Xa,Va,Aa) */
351:   TSComputeI2Jacobian(ts, ta, Xa, Va, Aa, dVdX, dAdX, J, P);
352:   return 0;
353: }

355: static PetscErrorCode TSReset_Alpha(TS ts)
356: {
357:   TS_Alpha *th = (TS_Alpha *)ts->data;

359:   VecDestroy(&th->X0);
360:   VecDestroy(&th->Xa);
361:   VecDestroy(&th->X1);
362:   VecDestroy(&th->V0);
363:   VecDestroy(&th->Va);
364:   VecDestroy(&th->V1);
365:   VecDestroy(&th->A0);
366:   VecDestroy(&th->Aa);
367:   VecDestroy(&th->A1);
368:   VecDestroy(&th->vec_sol_prev);
369:   VecDestroy(&th->vec_dot_prev);
370:   VecDestroy(&th->vec_lte_work[0]);
371:   VecDestroy(&th->vec_lte_work[1]);
372:   return 0;
373: }

375: static PetscErrorCode TSDestroy_Alpha(TS ts)
376: {
377:   TSReset_Alpha(ts);
378:   PetscFree(ts->data);

380:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", NULL);
381:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", NULL);
382:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", NULL);
383:   return 0;
384: }

386: static PetscErrorCode TSSetUp_Alpha(TS ts)
387: {
388:   TS_Alpha *th = (TS_Alpha *)ts->data;
389:   PetscBool match;

391:   VecDuplicate(ts->vec_sol, &th->X0);
392:   VecDuplicate(ts->vec_sol, &th->Xa);
393:   VecDuplicate(ts->vec_sol, &th->X1);
394:   VecDuplicate(ts->vec_sol, &th->V0);
395:   VecDuplicate(ts->vec_sol, &th->Va);
396:   VecDuplicate(ts->vec_sol, &th->V1);
397:   VecDuplicate(ts->vec_sol, &th->A0);
398:   VecDuplicate(ts->vec_sol, &th->Aa);
399:   VecDuplicate(ts->vec_sol, &th->A1);

401:   TSGetAdapt(ts, &ts->adapt);
402:   TSAdaptCandidatesClear(ts->adapt);
403:   PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match);
404:   if (!match) {
405:     VecDuplicate(ts->vec_sol, &th->vec_sol_prev);
406:     VecDuplicate(ts->vec_sol, &th->vec_dot_prev);
407:     VecDuplicate(ts->vec_sol, &th->vec_lte_work[0]);
408:     VecDuplicate(ts->vec_sol, &th->vec_lte_work[1]);
409:   }

411:   TSGetSNES(ts, &ts->snes);
412:   return 0;
413: }

415: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
416: {
417:   TS_Alpha *th = (TS_Alpha *)ts->data;

419:   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
420:   {
421:     PetscBool flg;
422:     PetscReal radius = 1;
423:     PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlpha2SetRadius", radius, &radius, &flg);
424:     if (flg) TSAlpha2SetRadius(ts, radius);
425:     PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlpha2SetParams", th->Alpha_m, &th->Alpha_m, NULL);
426:     PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlpha2SetParams", th->Alpha_f, &th->Alpha_f, NULL);
427:     PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlpha2SetParams", th->Gamma, &th->Gamma, NULL);
428:     PetscOptionsReal("-ts_alpha_beta", "Algorithmic parameter beta", "TSAlpha2SetParams", th->Beta, &th->Beta, NULL);
429:     TSAlpha2SetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma, th->Beta);
430:   }
431:   PetscOptionsHeadEnd();
432:   return 0;
433: }

435: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
436: {
437:   TS_Alpha *th = (TS_Alpha *)ts->data;
438:   PetscBool iascii;

440:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
441:   if (iascii) 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);
442:   return 0;
443: }

445: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
446: {
447:   PetscReal alpha_m, alpha_f, gamma, beta;

450:   alpha_m = (2 - radius) / (1 + radius);
451:   alpha_f = 1 / (1 + radius);
452:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
453:   beta    = (PetscReal)0.5 * (1 + alpha_m - alpha_f);
454:   beta *= beta;
455:   TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta);
456:   return 0;
457: }

459: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
460: {
461:   TS_Alpha *th  = (TS_Alpha *)ts->data;
462:   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
463:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

465:   th->Alpha_m = alpha_m;
466:   th->Alpha_f = alpha_f;
467:   th->Gamma   = gamma;
468:   th->Beta    = beta;
469:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
470:   return 0;
471: }

473: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
474: {
475:   TS_Alpha *th = (TS_Alpha *)ts->data;

477:   if (alpha_m) *alpha_m = th->Alpha_m;
478:   if (alpha_f) *alpha_f = th->Alpha_f;
479:   if (gamma) *gamma = th->Gamma;
480:   if (beta) *beta = th->Beta;
481:   return 0;
482: }

484: /*MC
485:       TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
486:                  for second-order systems

488:   Level: beginner

490:   References:
491: . * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
492:   Dynamics with Improved Numerical Dissipation: The Generalized-alpha
493:   Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.

495: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
496: M*/
497: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
498: {
499:   TS_Alpha *th;

501:   ts->ops->reset          = TSReset_Alpha;
502:   ts->ops->destroy        = TSDestroy_Alpha;
503:   ts->ops->view           = TSView_Alpha;
504:   ts->ops->setup          = TSSetUp_Alpha;
505:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
506:   ts->ops->step           = TSStep_Alpha;
507:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
508:   ts->ops->rollback       = TSRollBack_Alpha;
509:   /*ts->ops->interpolate  = TSInterpolate_Alpha;*/
510:   ts->ops->snesfunction  = SNESTSFormFunction_Alpha;
511:   ts->ops->snesjacobian  = SNESTSFormJacobian_Alpha;
512:   ts->default_adapt_type = TSADAPTNONE;

514:   ts->usessnes = PETSC_TRUE;

516:   PetscNew(&th);
517:   ts->data = (void *)th;

519:   th->Alpha_m = 0.5;
520:   th->Alpha_f = 0.5;
521:   th->Gamma   = 0.5;
522:   th->Beta    = 0.25;
523:   th->order   = 2;

525:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", TSAlpha2SetRadius_Alpha);
526:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", TSAlpha2SetParams_Alpha);
527:   PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", TSAlpha2GetParams_Alpha);
528:   return 0;
529: }

531: /*@
532:   TSAlpha2SetRadius - sets the desired spectral radius of the method for `TSALPHA2`
533:                       (i.e. high-frequency numerical damping)

535:   Logically Collective

537:   The algorithmic parameters \alpha_m and \alpha_f of the
538:   generalized-\alpha method can be computed in terms of a specified
539:   spectral radius \rho in [0,1] for infinite time step in order to
540:   control high-frequency numerical damping:
541:     \alpha_m = (2-\rho)/(1+\rho)
542:     \alpha_f = 1/(1+\rho)

544:   Input Parameters:
545: +  ts - timestepping context
546: -  radius - the desired spectral radius

548:   Options Database Key:
549: .  -ts_alpha_radius <radius> - set the desired spectral radius

551:   Level: intermediate

553: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
554: @*/
555: PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
556: {
560:   PetscTryMethod(ts, "TSAlpha2SetRadius_C", (TS, PetscReal), (ts, radius));
561:   return 0;
562: }

564: /*@
565:   TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`

567:   Logically Collective

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

573:   Unconditional stability requires:
574:     \alpha_m >= \alpha_f >= 1/2

576:   Input Parameters:
577: + ts - timestepping context
578: . \alpha_m - algorithmic parameter
579: . \alpha_f - algorithmic parameter
580: . \gamma   - algorithmic parameter
581: - \beta    - algorithmic parameter

583:   Options Database Keys:
584: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
585: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
586: . -ts_alpha_gamma   <gamma> - set gamma
587: - -ts_alpha_beta    <beta> - set beta

589:   Level: advanced

591:   Note:
592:   Use of this function is normally only required to hack `TSALPHA2` to
593:   use a modified integration scheme. Users should call
594:   `TSAlpha2SetRadius()` to set the desired spectral radius of the methods
595:   (i.e. high-frequency damping) in order so select optimal values for
596:   these parameters.

598: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
599: @*/
600: PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
601: {
607:   PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
608:   return 0;
609: }

611: /*@
612:   TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`

614:   Not Collective

616:   Input Parameter:
617: . ts - timestepping context

619:   Output Parameters:
620: + \alpha_m - algorithmic parameter
621: . \alpha_f - algorithmic parameter
622: . \gamma   - algorithmic parameter
623: - \beta    - algorithmic parameter

625:   Level: advanced

627:   Note:
628:   Use of this function is normally only required to hack `TSALPHA2` to
629:   use a modified integration scheme. Users should call
630:   `TSAlpha2SetRadius()` to set the high-frequency damping (i.e. spectral
631:   radius of the method) in order so select optimal values for these
632:   parameters.

634: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
635: @*/
636: PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
637: {
643:   PetscUseMethod(ts, "TSAlpha2GetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma, beta));
644:   return 0;
645: }