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 TSResizeRegister_Alpha(TS ts, PetscBool reg)
 38: {
 39:   TS_Alpha *th = (TS_Alpha *)ts->data;

 41:   PetscFunctionBegin;
 42:   if (reg) {
 43:     PetscCall(TSResizeRegisterVec(ts, "ts:theta:sol_prev", th->vec_sol_prev));
 44:     PetscCall(TSResizeRegisterVec(ts, "ts:theta:X0", th->X0));
 45:   } else {
 46:     PetscCall(TSResizeRetrieveVec(ts, "ts:theta:sol_prev", &th->vec_sol_prev));
 47:     PetscCall(PetscObjectReference((PetscObject)th->vec_sol_prev));
 48:     PetscCall(TSResizeRetrieveVec(ts, "ts:theta:X0", &th->X0));
 49:     PetscCall(PetscObjectReference((PetscObject)th->X0));
 50:   }
 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

 54: static PetscErrorCode TSAlpha_StageTime(TS ts)
 55: {
 56:   TS_Alpha *th      = (TS_Alpha *)ts->data;
 57:   PetscReal t       = ts->ptime;
 58:   PetscReal dt      = ts->time_step;
 59:   PetscReal Alpha_m = th->Alpha_m;
 60:   PetscReal Alpha_f = th->Alpha_f;
 61:   PetscReal Gamma   = th->Gamma;

 63:   PetscFunctionBegin;
 64:   th->stage_time = t + Alpha_f * dt;
 65:   th->shift_V    = Alpha_m / (Alpha_f * Gamma * dt);
 66:   th->scale_F    = 1 / Alpha_f;
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

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

 81:   PetscFunctionBegin;
 82:   /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
 83:   PetscCall(VecWAXPY(V1, -1.0, X0, X1));
 84:   PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
 85:   /* Xa = X0 + Alpha_f*(X1-X0) */
 86:   PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
 87:   PetscCall(VecAYPX(Xa, Alpha_f, X0));
 88:   /* Va = V0 + Alpha_m*(V1-V0) */
 89:   PetscCall(VecWAXPY(Va, -1.0, V0, V1));
 90:   PetscCall(VecAYPX(Va, Alpha_m, V0));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

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

 98:   PetscFunctionBegin;
 99:   PetscCall(SNESSolve(ts->snes, b, x));
100:   PetscCall(SNESGetIterationNumber(ts->snes, &nits));
101:   PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
102:   ts->snes_its += nits;
103:   ts->ksp_its += lits;
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

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

121:   PetscFunctionBegin;
122:   PetscCall(VecDuplicate(X0, &X1));

124:   /* Setup backward Euler with halved time step */
125:   PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
126:   PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
127:   PetscCall(TSGetTimeStep(ts, &time_step));
128:   ts->time_step = time_step / 2;
129:   PetscCall(TSAlpha_StageTime(ts));
130:   th->stage_time = ts->ptime;
131:   PetscCall(VecZeroEntries(th->V0));

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

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

153:   /* Compute V0 ~ dX/dt at t0 with backward differences */
154:   PetscCall(VecZeroEntries(th->V0));
155:   PetscCall(VecAXPY(th->V0, -3 / ts->time_step, X0));
156:   PetscCall(VecAXPY(th->V0, +4 / ts->time_step, X1));
157:   PetscCall(VecAXPY(th->V0, -1 / ts->time_step, X2));

159:   /* Rough, lower-order estimate LTE of the initial step */
160:   if (th->vec_lte_work) {
161:     PetscCall(VecZeroEntries(th->vec_lte_work));
162:     PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
163:     PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
164:     PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
165:   }

167: finally:
168:   /* Revert TSAlpha to the initial state (t0,X0) */
169:   if (initok) *initok = stageok;
170:   PetscCall(TSSetTimeStep(ts, time_step));
171:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
172:   PetscCall(VecCopy(ts->vec_sol, th->X0));

174:   PetscCall(VecDestroy(&X1));
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: static PetscErrorCode TSStep_Alpha(TS ts)
179: {
180:   TS_Alpha *th         = (TS_Alpha *)ts->data;
181:   PetscInt  rejections = 0;
182:   PetscBool stageok, accept = PETSC_TRUE;
183:   PetscReal next_time_step = ts->time_step;

185:   PetscFunctionBegin;
186:   PetscCall(PetscCitationsRegister(citation, &cited));

188:   if (!ts->steprollback) {
189:     if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
190:     PetscCall(VecCopy(ts->vec_sol, th->X0));
191:     PetscCall(VecCopy(th->V1, th->V0));
192:   }

194:   th->status = TS_STEP_INCOMPLETE;
195:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
196:     if (ts->steprestart) {
197:       PetscCall(TSAlpha_Restart(ts, &stageok));
198:       if (!stageok) goto reject_step;
199:     }

201:     PetscCall(TSAlpha_StageTime(ts));
202:     PetscCall(VecCopy(th->X0, th->X1));
203:     PetscCall(TSPreStage(ts, th->stage_time));
204:     PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
205:     PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
206:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
207:     if (!stageok) goto reject_step;

209:     th->status = TS_STEP_PENDING;
210:     PetscCall(VecCopy(th->X1, ts->vec_sol));
211:     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
212:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
213:     if (!accept) {
214:       PetscCall(VecCopy(th->X0, ts->vec_sol));
215:       ts->time_step = next_time_step;
216:       goto reject_step;
217:     }

219:     ts->ptime += ts->time_step;
220:     ts->time_step = next_time_step;
221:     break;

223:   reject_step:
224:     ts->reject++;
225:     accept = PETSC_FALSE;
226:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
227:       ts->reason = TS_DIVERGED_STEP_REJECTED;
228:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
229:     }
230:   }
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
235: {
236:   TS_Alpha *th = (TS_Alpha *)ts->data;
237:   Vec       X  = th->X1;           /* X = solution */
238:   Vec       Y  = th->vec_lte_work; /* Y = X + LTE  */
239:   PetscReal wltea, wlter;

241:   PetscFunctionBegin;
242:   if (!th->vec_sol_prev) {
243:     *wlte = -1;
244:     PetscFunctionReturn(PETSC_SUCCESS);
245:   }
246:   if (!th->vec_lte_work) {
247:     *wlte = -1;
248:     PetscFunctionReturn(PETSC_SUCCESS);
249:   }
250:   if (ts->steprestart) {
251:     /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
252:     PetscCall(VecAXPY(Y, 1, X));
253:   } else {
254:     /* Compute LTE using backward differences with non-constant time step */
255:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
256:     PetscReal   a = 1 + h_prev / h;
257:     PetscScalar scal[3];
258:     Vec         vecs[3];
259:     scal[0] = +1 / a;
260:     scal[1] = -1 / (a - 1);
261:     scal[2] = +1 / (a * (a - 1));
262:     vecs[0] = th->X1;
263:     vecs[1] = th->X0;
264:     vecs[2] = th->vec_sol_prev;
265:     PetscCall(VecCopy(X, Y));
266:     PetscCall(VecMAXPY(Y, 3, scal, vecs));
267:   }
268:   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
269:   if (order) *order = 2;
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
274: {
275:   TS_Alpha *th    = (TS_Alpha *)ts->data;
276:   PetscReal dt    = t - ts->ptime;
277:   PetscReal Gamma = th->Gamma;

279:   PetscFunctionBegin;
280:   PetscCall(VecWAXPY(th->V1, -1.0, th->X0, ts->vec_sol));
281:   PetscCall(VecAXPBY(th->V1, 1 - 1 / Gamma, 1 / (Gamma * ts->time_step), th->V0));
282:   PetscCall(VecCopy(ts->vec_sol, X));
283:   /* X = X + Gamma*dT*V1 */
284:   PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
285:   /* X = X + (1-Gamma)*dT*V0 */
286:   PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
287:   PetscFunctionReturn(PETSC_SUCCESS);
288: }

290: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
291: {
292:   TS_Alpha *th = (TS_Alpha *)ts->data;
293:   PetscReal ta = th->stage_time;
294:   Vec       Xa = th->Xa, Va = th->Va;

296:   PetscFunctionBegin;
297:   PetscCall(TSAlpha_StageVecs(ts, X));
298:   /* F = Function(ta,Xa,Va) */
299:   PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
300:   PetscCall(VecScale(F, th->scale_F));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
305: {
306:   TS_Alpha *th = (TS_Alpha *)ts->data;
307:   PetscReal ta = th->stage_time;
308:   Vec       Xa = th->Xa, Va = th->Va;
309:   PetscReal dVdX = th->shift_V;

311:   PetscFunctionBegin;
312:   /* J,P = Jacobian(ta,Xa,Va) */
313:   PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: static PetscErrorCode TSReset_Alpha(TS ts)
318: {
319:   TS_Alpha *th = (TS_Alpha *)ts->data;

321:   PetscFunctionBegin;
322:   PetscCall(VecDestroy(&th->X0));
323:   PetscCall(VecDestroy(&th->Xa));
324:   PetscCall(VecDestroy(&th->X1));
325:   PetscCall(VecDestroy(&th->V0));
326:   PetscCall(VecDestroy(&th->Va));
327:   PetscCall(VecDestroy(&th->V1));
328:   PetscCall(VecDestroy(&th->vec_sol_prev));
329:   PetscCall(VecDestroy(&th->vec_lte_work));
330:   PetscFunctionReturn(PETSC_SUCCESS);
331: }

333: static PetscErrorCode TSDestroy_Alpha(TS ts)
334: {
335:   PetscFunctionBegin;
336:   PetscCall(TSReset_Alpha(ts));
337:   PetscCall(PetscFree(ts->data));

339:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
340:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
341:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
342:   PetscFunctionReturn(PETSC_SUCCESS);
343: }

345: static PetscErrorCode TSSetUp_Alpha(TS ts)
346: {
347:   TS_Alpha *th = (TS_Alpha *)ts->data;
348:   PetscBool match;

350:   PetscFunctionBegin;
351:   if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
352:   PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
353:   PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
354:   PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
355:   PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
356:   PetscCall(VecDuplicate(ts->vec_sol, &th->V1));

358:   PetscCall(TSGetAdapt(ts, &ts->adapt));
359:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
360:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
361:   if (!match) {
362:     if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
363:     if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
364:   }

366:   PetscCall(TSGetSNES(ts, &ts->snes));
367:   PetscFunctionReturn(PETSC_SUCCESS);
368: }

370: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
371: {
372:   TS_Alpha *th = (TS_Alpha *)ts->data;

374:   PetscFunctionBegin;
375:   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
376:   {
377:     PetscBool flg;
378:     PetscReal radius = 1;
379:     PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
380:     if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
381:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
382:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
383:     PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
384:     PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
385:   }
386:   PetscOptionsHeadEnd();
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
391: {
392:   TS_Alpha *th = (TS_Alpha *)ts->data;
393:   PetscBool iascii;

395:   PetscFunctionBegin;
396:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
397:   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));
398:   PetscFunctionReturn(PETSC_SUCCESS);
399: }

401: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
402: {
403:   PetscReal alpha_m, alpha_f, gamma;

405:   PetscFunctionBegin;
406:   PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
407:   alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
408:   alpha_f = 1 / (1 + radius);
409:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
410:   PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
411:   PetscFunctionReturn(PETSC_SUCCESS);
412: }

414: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
415: {
416:   TS_Alpha *th  = (TS_Alpha *)ts->data;
417:   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
418:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

420:   PetscFunctionBegin;
421:   th->Alpha_m = alpha_m;
422:   th->Alpha_f = alpha_f;
423:   th->Gamma   = gamma;
424:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

428: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
429: {
430:   TS_Alpha *th = (TS_Alpha *)ts->data;

432:   PetscFunctionBegin;
433:   if (alpha_m) *alpha_m = th->Alpha_m;
434:   if (alpha_f) *alpha_f = th->Alpha_f;
435:   if (gamma) *gamma = th->Gamma;
436:   PetscFunctionReturn(PETSC_SUCCESS);
437: }

439: /*MC
440:   TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems

442:   Level: beginner

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

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

464:   ts->usessnes = PETSC_TRUE;

466:   PetscCall(PetscNew(&th));
467:   ts->data = (void *)th;

469:   th->Alpha_m = 0.5;
470:   th->Alpha_f = 0.5;
471:   th->Gamma   = 0.5;
472:   th->order   = 2;

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

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

484:   Logically Collective

486:   Input Parameters:
487: + ts     - timestepping context
488: - radius - the desired spectral radius

490:   Options Database Key:
491: . -ts_alpha_radius <radius> - set alpha radius

493:   Level: intermediate

495:   Notes:
496:   The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
497:   be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
498:   in order to control high-frequency numerical damping\:

500:   $$
501:   \begin{align*}
502:   \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
503:   \alpha_f = 1/(1+\rho)
504:   \end{align*}
505:   $$

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

519: /*@
520:   TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`

522:   Logically Collective

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

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

535:   Level: advanced

537:   Note:
538:   Second-order accuracy can be obtained so long as\:  $\gamma = 0.5 + \alpha_m - \alpha_f$

540:   Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$

542:   Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$

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

549: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
550: @*/
551: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
552: {
553:   PetscFunctionBegin;
558:   PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: /*@
563:   TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`

565:   Not Collective

567:   Input Parameter:
568: . ts - timestepping context

570:   Output Parameters:
571: + alpha_m - algorithmic parameter
572: . alpha_f - algorithmic parameter
573: - gamma   - algorithmic parameter

575:   Level: advanced

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

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