Actual source code: alpha1.c

  1: /*
  2:   Code for timestepping with implicit generalized-\alpha method
  3:   for first order systems.
  4: */
  5: #include <petsc/private/tsimpl.h>

  7: static PetscBool  cited      = PETSC_FALSE;
  8: static const char citation[] = "@article{Jansen2000,\n"
  9:                                "  title   = {A generalized-$\\alpha$ method for integrating the filtered {N}avier--{S}tokes equations with a stabilized finite element method},\n"
 10:                                "  author  = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},\n"
 11:                                "  journal = {Computer Methods in Applied Mechanics and Engineering},\n"
 12:                                "  volume  = {190},\n"
 13:                                "  number  = {3--4},\n"
 14:                                "  pages   = {305--319},\n"
 15:                                "  year    = {2000},\n"
 16:                                "  issn    = {0045-7825},\n"
 17:                                "  doi     = {http://dx.doi.org/10.1016/S0045-7825(00)00203-6}\n}\n";

 19: typedef struct {
 20:   PetscReal stage_time;
 21:   PetscReal shift_V;
 22:   PetscReal scale_F;
 23:   Vec       X0, Xa, X1;
 24:   Vec       V0, Va, V1;

 26:   PetscReal Alpha_m;
 27:   PetscReal Alpha_f;
 28:   PetscReal Gamma;
 29:   PetscInt  order;

 31:   Vec vec_sol_prev;
 32:   Vec vec_lte_work;

 34:   TSStepStatus status;
 35: } TS_Alpha;

 37: static PetscErrorCode TSAlpha_StageTime(TS ts)
 38: {
 39:   TS_Alpha *th      = (TS_Alpha *)ts->data;
 40:   PetscReal t       = ts->ptime;
 41:   PetscReal dt      = ts->time_step;
 42:   PetscReal Alpha_m = th->Alpha_m;
 43:   PetscReal Alpha_f = th->Alpha_f;
 44:   PetscReal Gamma   = th->Gamma;

 46:   PetscFunctionBegin;
 47:   th->stage_time = t + Alpha_f * dt;
 48:   th->shift_V    = Alpha_m / (Alpha_f * Gamma * dt);
 49:   th->scale_F    = 1 / Alpha_f;
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
 54: {
 55:   TS_Alpha *th = (TS_Alpha *)ts->data;
 56:   Vec       X1 = X, V1 = th->V1;
 57:   Vec       Xa = th->Xa, Va = th->Va;
 58:   Vec       X0 = th->X0, V0 = th->V0;
 59:   PetscReal dt      = ts->time_step;
 60:   PetscReal Alpha_m = th->Alpha_m;
 61:   PetscReal Alpha_f = th->Alpha_f;
 62:   PetscReal Gamma   = th->Gamma;

 64:   PetscFunctionBegin;
 65:   /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
 66:   PetscCall(VecWAXPY(V1, -1.0, X0, X1));
 67:   PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
 68:   /* Xa = X0 + Alpha_f*(X1-X0) */
 69:   PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
 70:   PetscCall(VecAYPX(Xa, Alpha_f, X0));
 71:   /* Va = V0 + Alpha_m*(V1-V0) */
 72:   PetscCall(VecWAXPY(Va, -1.0, V0, V1));
 73:   PetscCall(VecAYPX(Va, Alpha_m, V0));
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
 78: {
 79:   PetscInt nits, lits;

 81:   PetscFunctionBegin;
 82:   PetscCall(SNESSolve(ts->snes, b, x));
 83:   PetscCall(SNESGetIterationNumber(ts->snes, &nits));
 84:   PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
 85:   ts->snes_its += nits;
 86:   ts->ksp_its += lits;
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

 90: /*
 91:   Compute a consistent initial state for the generalized-alpha method.
 92:   - Solve two successive backward Euler steps with halved time step.
 93:   - Compute the initial time derivative using backward differences.
 94:   - If using adaptivity, estimate the LTE of the initial step.
 95: */
 96: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
 97: {
 98:   TS_Alpha *th = (TS_Alpha *)ts->data;
 99:   PetscReal time_step;
100:   PetscReal alpha_m, alpha_f, gamma;
101:   Vec       X0 = ts->vec_sol, X1, X2 = th->X1;
102:   PetscBool stageok;

104:   PetscFunctionBegin;
105:   PetscCall(VecDuplicate(X0, &X1));

107:   /* Setup backward Euler with halved time step */
108:   PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
109:   PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
110:   PetscCall(TSGetTimeStep(ts, &time_step));
111:   ts->time_step = time_step / 2;
112:   PetscCall(TSAlpha_StageTime(ts));
113:   th->stage_time = ts->ptime;
114:   PetscCall(VecZeroEntries(th->V0));

116:   /* First BE step, (t0,X0) -> (t1,X1) */
117:   th->stage_time += ts->time_step;
118:   PetscCall(VecCopy(X0, th->X0));
119:   PetscCall(TSPreStage(ts, th->stage_time));
120:   PetscCall(VecCopy(th->X0, X1));
121:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
122:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
123:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
124:   if (!stageok) goto finally;

126:   /* Second BE step, (t1,X1) -> (t2,X2) */
127:   th->stage_time += ts->time_step;
128:   PetscCall(VecCopy(X1, th->X0));
129:   PetscCall(TSPreStage(ts, th->stage_time));
130:   PetscCall(VecCopy(th->X0, X2));
131:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
132:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
133:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
134:   if (!stageok) goto finally;

136:   /* Compute V0 ~ dX/dt at t0 with backward differences */
137:   PetscCall(VecZeroEntries(th->V0));
138:   PetscCall(VecAXPY(th->V0, -3 / ts->time_step, X0));
139:   PetscCall(VecAXPY(th->V0, +4 / ts->time_step, X1));
140:   PetscCall(VecAXPY(th->V0, -1 / ts->time_step, X2));

142:   /* Rough, lower-order estimate LTE of the initial step */
143:   if (th->vec_lte_work) {
144:     PetscCall(VecZeroEntries(th->vec_lte_work));
145:     PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
146:     PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
147:     PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
148:   }

150: finally:
151:   /* Revert TSAlpha to the initial state (t0,X0) */
152:   if (initok) *initok = stageok;
153:   PetscCall(TSSetTimeStep(ts, time_step));
154:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
155:   PetscCall(VecCopy(ts->vec_sol, th->X0));

157:   PetscCall(VecDestroy(&X1));
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: static PetscErrorCode TSStep_Alpha(TS ts)
162: {
163:   TS_Alpha *th         = (TS_Alpha *)ts->data;
164:   PetscInt  rejections = 0;
165:   PetscBool stageok, accept = PETSC_TRUE;
166:   PetscReal next_time_step = ts->time_step;

168:   PetscFunctionBegin;
169:   PetscCall(PetscCitationsRegister(citation, &cited));

171:   if (!ts->steprollback) {
172:     if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
173:     PetscCall(VecCopy(ts->vec_sol, th->X0));
174:     PetscCall(VecCopy(th->V1, th->V0));
175:   }

177:   th->status = TS_STEP_INCOMPLETE;
178:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
179:     if (ts->steprestart) {
180:       PetscCall(TSAlpha_Restart(ts, &stageok));
181:       if (!stageok) goto reject_step;
182:     }

184:     PetscCall(TSAlpha_StageTime(ts));
185:     PetscCall(VecCopy(th->X0, th->X1));
186:     PetscCall(TSPreStage(ts, th->stage_time));
187:     PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
188:     PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
189:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
190:     if (!stageok) goto reject_step;

192:     th->status = TS_STEP_PENDING;
193:     PetscCall(VecCopy(th->X1, ts->vec_sol));
194:     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
195:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
196:     if (!accept) {
197:       PetscCall(VecCopy(th->X0, ts->vec_sol));
198:       ts->time_step = next_time_step;
199:       goto reject_step;
200:     }

202:     ts->ptime += ts->time_step;
203:     ts->time_step = next_time_step;
204:     break;

206:   reject_step:
207:     ts->reject++;
208:     accept = PETSC_FALSE;
209:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
210:       ts->reason = TS_DIVERGED_STEP_REJECTED;
211:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
212:     }
213:   }
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
218: {
219:   TS_Alpha *th = (TS_Alpha *)ts->data;
220:   Vec       X  = th->X1;           /* X = solution */
221:   Vec       Y  = th->vec_lte_work; /* Y = X + LTE  */
222:   PetscReal wltea, wlter;

224:   PetscFunctionBegin;
225:   if (!th->vec_sol_prev) {
226:     *wlte = -1;
227:     PetscFunctionReturn(PETSC_SUCCESS);
228:   }
229:   if (!th->vec_lte_work) {
230:     *wlte = -1;
231:     PetscFunctionReturn(PETSC_SUCCESS);
232:   }
233:   if (ts->steprestart) {
234:     /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
235:     PetscCall(VecAXPY(Y, 1, X));
236:   } else {
237:     /* Compute LTE using backward differences with non-constant time step */
238:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
239:     PetscReal   a = 1 + h_prev / h;
240:     PetscScalar scal[3];
241:     Vec         vecs[3];
242:     scal[0] = +1 / a;
243:     scal[1] = -1 / (a - 1);
244:     scal[2] = +1 / (a * (a - 1));
245:     vecs[0] = th->X1;
246:     vecs[1] = th->X0;
247:     vecs[2] = th->vec_sol_prev;
248:     PetscCall(VecCopy(X, Y));
249:     PetscCall(VecMAXPY(Y, 3, scal, vecs));
250:   }
251:   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
252:   if (order) *order = 2;
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: static PetscErrorCode TSRollBack_Alpha(TS ts)
257: {
258:   TS_Alpha *th = (TS_Alpha *)ts->data;

260:   PetscFunctionBegin;
261:   PetscCall(VecCopy(th->X0, ts->vec_sol));
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

265: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
266: {
267:   TS_Alpha *th = (TS_Alpha *)ts->data;
268:   PetscReal dt = t - ts->ptime;

270:   PetscFunctionBegin;
271:   PetscCall(VecCopy(ts->vec_sol, X));
272:   PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
273:   PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
278: {
279:   TS_Alpha *th = (TS_Alpha *)ts->data;
280:   PetscReal ta = th->stage_time;
281:   Vec       Xa = th->Xa, Va = th->Va;

283:   PetscFunctionBegin;
284:   PetscCall(TSAlpha_StageVecs(ts, X));
285:   /* F = Function(ta,Xa,Va) */
286:   PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
287:   PetscCall(VecScale(F, th->scale_F));
288:   PetscFunctionReturn(PETSC_SUCCESS);
289: }

291: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
292: {
293:   TS_Alpha *th = (TS_Alpha *)ts->data;
294:   PetscReal ta = th->stage_time;
295:   Vec       Xa = th->Xa, Va = th->Va;
296:   PetscReal dVdX = th->shift_V;

298:   PetscFunctionBegin;
299:   /* J,P = Jacobian(ta,Xa,Va) */
300:   PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

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

308:   PetscFunctionBegin;
309:   PetscCall(VecDestroy(&th->X0));
310:   PetscCall(VecDestroy(&th->Xa));
311:   PetscCall(VecDestroy(&th->X1));
312:   PetscCall(VecDestroy(&th->V0));
313:   PetscCall(VecDestroy(&th->Va));
314:   PetscCall(VecDestroy(&th->V1));
315:   PetscCall(VecDestroy(&th->vec_sol_prev));
316:   PetscCall(VecDestroy(&th->vec_lte_work));
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: static PetscErrorCode TSDestroy_Alpha(TS ts)
321: {
322:   PetscFunctionBegin;
323:   PetscCall(TSReset_Alpha(ts));
324:   PetscCall(PetscFree(ts->data));

326:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
327:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
328:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: static PetscErrorCode TSSetUp_Alpha(TS ts)
333: {
334:   TS_Alpha *th = (TS_Alpha *)ts->data;
335:   PetscBool match;

337:   PetscFunctionBegin;
338:   PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
339:   PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
340:   PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
341:   PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
342:   PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
343:   PetscCall(VecDuplicate(ts->vec_sol, &th->V1));

345:   PetscCall(TSGetAdapt(ts, &ts->adapt));
346:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
347:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
348:   if (!match) {
349:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
350:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
351:   }

353:   PetscCall(TSGetSNES(ts, &ts->snes));
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
358: {
359:   TS_Alpha *th = (TS_Alpha *)ts->data;

361:   PetscFunctionBegin;
362:   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
363:   {
364:     PetscBool flg;
365:     PetscReal radius = 1;
366:     PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
367:     if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
368:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
369:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
370:     PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
371:     PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
372:   }
373:   PetscOptionsHeadEnd();
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }

377: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
378: {
379:   TS_Alpha *th = (TS_Alpha *)ts->data;
380:   PetscBool iascii;

382:   PetscFunctionBegin;
383:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
384:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  Alpha_m=%g, Alpha_f=%g, Gamma=%g\n", (double)th->Alpha_m, (double)th->Alpha_f, (double)th->Gamma));
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
389: {
390:   PetscReal alpha_m, alpha_f, gamma;

392:   PetscFunctionBegin;
393:   PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
394:   alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
395:   alpha_f = 1 / (1 + radius);
396:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
397:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
398:   PetscFunctionReturn(PETSC_SUCCESS);
399: }

401: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
402: {
403:   TS_Alpha *th  = (TS_Alpha *)ts->data;
404:   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
405:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

407:   PetscFunctionBegin;
408:   th->Alpha_m = alpha_m;
409:   th->Alpha_f = alpha_f;
410:   th->Gamma   = gamma;
411:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
416: {
417:   TS_Alpha *th = (TS_Alpha *)ts->data;

419:   PetscFunctionBegin;
420:   if (alpha_m) *alpha_m = th->Alpha_m;
421:   if (alpha_f) *alpha_f = th->Alpha_f;
422:   if (gamma) *gamma = th->Gamma;
423:   PetscFunctionReturn(PETSC_SUCCESS);
424: }

426: /*MC
427:       TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method
428:                 for first-order systems

430:   Level: beginner

432:   References:
433: + * - K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
434:   method for integrating the filtered Navier-Stokes equations with a
435:   stabilized finite element method", Computer Methods in Applied
436:   Mechanics and Engineering, 190, 305-319, 2000.
437:   DOI: 10.1016/S0045-7825(00)00203-6.
438: - * -  J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
439:   Dynamics with Improved Numerical Dissipation: The Generalized-alpha
440:   Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.

442: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
443: M*/
444: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
445: {
446:   TS_Alpha *th;

448:   PetscFunctionBegin;
449:   ts->ops->reset          = TSReset_Alpha;
450:   ts->ops->destroy        = TSDestroy_Alpha;
451:   ts->ops->view           = TSView_Alpha;
452:   ts->ops->setup          = TSSetUp_Alpha;
453:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
454:   ts->ops->step           = TSStep_Alpha;
455:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
456:   ts->ops->rollback       = TSRollBack_Alpha;
457:   ts->ops->interpolate    = TSInterpolate_Alpha;
458:   ts->ops->snesfunction   = SNESTSFormFunction_Alpha;
459:   ts->ops->snesjacobian   = SNESTSFormJacobian_Alpha;
460:   ts->default_adapt_type  = TSADAPTNONE;

462:   ts->usessnes = PETSC_TRUE;

464:   PetscCall(PetscNew(&th));
465:   ts->data = (void *)th;

467:   th->Alpha_m = 0.5;
468:   th->Alpha_f = 0.5;
469:   th->Gamma   = 0.5;
470:   th->order   = 2;

472:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
473:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
474:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

478: /*@
479:   TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
480:                      (i.e. high-frequency numerical damping)

482:   Logically Collective

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

491:   Input Parameters:
492: +  ts - timestepping context
493: -  radius - the desired spectral radius

495:   Options Database Key:
496: .  -ts_alpha_radius <radius> - set alpha radius

498:   Level: intermediate

500: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
501: @*/
502: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
503: {
504:   PetscFunctionBegin;
507:   PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
508:   PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

512: /*@
513:   TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`

515:   Logically Collective

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

520:   Unconditional stability requires:
521:     \alpha_m >= \alpha_f >= 0.5

523:   Backward Euler method is recovered with:
524:     \alpha_m = \alpha_f = gamma = 1

526:   Input Parameters:
527: +  ts - timestepping context
528: .  alpha_m - algorithmic parameter
529: .  alpha_f - algorithmic parameter
530: -  gamma   - algorithmic parameter

532:    Options Database Keys:
533: +  -ts_alpha_alpha_m <alpha_m> - set alpha_m
534: .  -ts_alpha_alpha_f <alpha_f> - set alpha_f
535: -  -ts_alpha_gamma   <gamma> - set gamma

537:   Level: advanced

539:   Note:
540:   Use of this function is normally only required to hack `TSALPHA` to
541:   use a modified integration scheme. Users should call
542:   `TSAlphaSetRadius()` to set the desired spectral radius of the methods
543:   (i.e. high-frequency damping) in order so select optimal values for
544:   these parameters.

546: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
547: @*/
548: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
549: {
550:   PetscFunctionBegin;
555:   PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

559: /*@
560:   TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`

562:   Not Collective

564:   Input Parameter:
565: .  ts - timestepping context

567:   Output Parameters:
568: +  alpha_m - algorithmic parameter
569: .  alpha_f - algorithmic parameter
570: -  gamma   - algorithmic parameter

572:   Level: advanced

574:   Note:
575:   Use of this function is normally only required to hack `TSALPHA` to
576:   use a modified integration scheme. Users should call
577:   `TSAlphaSetRadius()` to set the high-frequency damping (i.e. spectral
578:   radius of the method) in order so select optimal values for these
579:   parameters.

581: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
582: @*/
583: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
584: {
585:   PetscFunctionBegin;
590:   PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }