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;

 42:   TSAlpha2PredictorFn *predictor;
 43:   void                *predictor_ctx;
 44: } TS_Alpha;

 46: /*@C
 47:   TSAlpha2SetPredictor - sets the callback for computing a predictor (i.e., initial guess
 48:   for the nonlinear solver).

 50:   Input Parameters:
 51: + ts        - timestepping context
 52: . predictor - callback to set the predictor in each step
 53: - ctx       - the application context, which may be set to `NULL` if not used

 55:   Level: intermediate

 57:   Notes:

 59:   If this function is never called, a same-state-vector predictor will be used, i.e.,
 60:   the initial guess will be the converged solution from the previous time step, without regard
 61:   for the previous velocity or acceleration.

 63: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2PredictorFn`
 64: @*/
 65: PetscErrorCode TSAlpha2SetPredictor(TS ts, TSAlpha2PredictorFn *predictor, void *ctx)
 66: {
 67:   TS_Alpha *th = (TS_Alpha *)ts->data;

 69:   PetscFunctionBegin;
 70:   th->predictor     = predictor;
 71:   th->predictor_ctx = ctx;
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: static PetscErrorCode TSAlpha_ApplyPredictor(TS ts, Vec X1)
 76: {
 77:   /* Apply a custom predictor if set, or default to same-displacement. */
 78:   TS_Alpha *th = (TS_Alpha *)ts->data;

 80:   PetscFunctionBegin;
 81:   if (th->predictor) PetscCall(th->predictor(ts, th->X0, th->V0, th->A0, X1, th->predictor_ctx));
 82:   else PetscCall(VecCopy(th->X0, X1));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: static PetscErrorCode TSAlpha_StageTime(TS ts)
 87: {
 88:   TS_Alpha *th      = (TS_Alpha *)ts->data;
 89:   PetscReal t       = ts->ptime;
 90:   PetscReal dt      = ts->time_step;
 91:   PetscReal Alpha_m = th->Alpha_m;
 92:   PetscReal Alpha_f = th->Alpha_f;
 93:   PetscReal Gamma   = th->Gamma;
 94:   PetscReal Beta    = th->Beta;

 96:   PetscFunctionBegin;
 97:   th->stage_time = t + Alpha_f * dt;
 98:   th->shift_V    = Gamma / (dt * Beta);
 99:   th->shift_A    = Alpha_m / (Alpha_f * dt * dt * Beta);
100:   th->scale_F    = 1 / Alpha_f;
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
105: {
106:   TS_Alpha *th = (TS_Alpha *)ts->data;
107:   Vec       X1 = X, V1 = th->V1, A1 = th->A1;
108:   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
109:   Vec       X0 = th->X0, V0 = th->V0, A0 = th->A0;
110:   PetscReal dt      = ts->time_step;
111:   PetscReal Alpha_m = th->Alpha_m;
112:   PetscReal Alpha_f = th->Alpha_f;
113:   PetscReal Gamma   = th->Gamma;
114:   PetscReal Beta    = th->Beta;

116:   PetscFunctionBegin;
117:   /* A1 = ... */
118:   PetscCall(VecWAXPY(A1, -1.0, X0, X1));
119:   PetscCall(VecAXPY(A1, -dt, V0));
120:   PetscCall(VecAXPBY(A1, -(1 - 2 * Beta) / (2 * Beta), 1 / (dt * dt * Beta), A0));
121:   /* V1 = ... */
122:   PetscCall(VecWAXPY(V1, (1.0 - Gamma) / Gamma, A0, A1));
123:   PetscCall(VecAYPX(V1, dt * Gamma, V0));
124:   /* Xa = X0 + Alpha_f*(X1-X0) */
125:   PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
126:   PetscCall(VecAYPX(Xa, Alpha_f, X0));
127:   /* Va = V0 + Alpha_f*(V1-V0) */
128:   PetscCall(VecWAXPY(Va, -1.0, V0, V1));
129:   PetscCall(VecAYPX(Va, Alpha_f, V0));
130:   /* Aa = A0 + Alpha_m*(A1-A0) */
131:   PetscCall(VecWAXPY(Aa, -1.0, A0, A1));
132:   PetscCall(VecAYPX(Aa, Alpha_m, A0));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
137: {
138:   PetscInt nits, lits;

140:   PetscFunctionBegin;
141:   PetscCall(SNESSolve(ts->snes, b, x));
142:   PetscCall(SNESGetIterationNumber(ts->snes, &nits));
143:   PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
144:   ts->snes_its += nits;
145:   ts->ksp_its += lits;
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }

149: /*
150:   Compute a consistent initial state for the generalized-alpha method.
151:   - Solve two successive backward Euler steps with halved time step.
152:   - Compute the initial second time derivative using backward differences.
153:   - If using adaptivity, estimate the LTE of the initial step.
154: */
155: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
156: {
157:   TS_Alpha *th = (TS_Alpha *)ts->data;
158:   PetscReal time_step;
159:   PetscReal alpha_m, alpha_f, gamma, beta;
160:   Vec       X0 = ts->vec_sol, X1, X2 = th->X1;
161:   Vec       V0 = ts->vec_dot, V1, V2 = th->V1;
162:   PetscBool stageok;

164:   PetscFunctionBegin;
165:   PetscCall(VecDuplicate(X0, &X1));
166:   PetscCall(VecDuplicate(V0, &V1));

168:   /* Setup backward Euler with halved time step */
169:   PetscCall(TSAlpha2GetParams(ts, &alpha_m, &alpha_f, &gamma, &beta));
170:   PetscCall(TSAlpha2SetParams(ts, 1, 1, 1, 0.5));
171:   PetscCall(TSGetTimeStep(ts, &time_step));
172:   ts->time_step = time_step / 2;
173:   PetscCall(TSAlpha_StageTime(ts));
174:   th->stage_time = ts->ptime;
175:   PetscCall(VecZeroEntries(th->A0));

177:   /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
178:   th->stage_time += ts->time_step;
179:   PetscCall(VecCopy(X0, th->X0));
180:   PetscCall(VecCopy(V0, th->V0));
181:   PetscCall(TSPreStage(ts, th->stage_time));
182:   PetscCall(TSAlpha_ApplyPredictor(ts, X1));
183:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
184:   PetscCall(VecCopy(th->V1, V1));
185:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
186:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
187:   if (!stageok) goto finally;

189:   /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
190:   th->stage_time += ts->time_step;
191:   PetscCall(VecCopy(X1, th->X0));
192:   PetscCall(VecCopy(V1, th->V0));
193:   PetscCall(TSPreStage(ts, th->stage_time));
194:   PetscCall(TSAlpha_ApplyPredictor(ts, X2));
195:   PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
196:   PetscCall(VecCopy(th->V1, V2));
197:   PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
198:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
199:   if (!stageok) goto finally;

201:   /* Compute A0 ~ dV/dt at t0 with backward differences */
202:   PetscCall(VecZeroEntries(th->A0));
203:   PetscCall(VecAXPY(th->A0, -3 / ts->time_step, V0));
204:   PetscCall(VecAXPY(th->A0, +4 / ts->time_step, V1));
205:   PetscCall(VecAXPY(th->A0, -1 / ts->time_step, V2));

207:   /* Rough, lower-order estimate LTE of the initial step */
208:   if (th->vec_lte_work[0]) {
209:     PetscCall(VecZeroEntries(th->vec_lte_work[0]));
210:     PetscCall(VecAXPY(th->vec_lte_work[0], +2, X2));
211:     PetscCall(VecAXPY(th->vec_lte_work[0], -4, X1));
212:     PetscCall(VecAXPY(th->vec_lte_work[0], +2, X0));
213:   }
214:   if (th->vec_lte_work[1]) {
215:     PetscCall(VecZeroEntries(th->vec_lte_work[1]));
216:     PetscCall(VecAXPY(th->vec_lte_work[1], +2, V2));
217:     PetscCall(VecAXPY(th->vec_lte_work[1], -4, V1));
218:     PetscCall(VecAXPY(th->vec_lte_work[1], +2, V0));
219:   }

221: finally:
222:   /* Revert TSAlpha to the initial state (t0,X0,V0), but retain
223:      potential time step reduction by factor after failed solve. */
224:   if (initok) *initok = stageok;
225:   PetscCall(TSSetTimeStep(ts, 2 * ts->time_step));
226:   PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
227:   PetscCall(VecCopy(ts->vec_sol, th->X0));
228:   PetscCall(VecCopy(ts->vec_dot, th->V0));

230:   PetscCall(VecDestroy(&X1));
231:   PetscCall(VecDestroy(&V1));
232:   PetscFunctionReturn(PETSC_SUCCESS);
233: }

235: static PetscErrorCode TSStep_Alpha(TS ts)
236: {
237:   TS_Alpha *th         = (TS_Alpha *)ts->data;
238:   PetscInt  rejections = 0;
239:   PetscBool stageok, accept = PETSC_TRUE;
240:   PetscReal next_time_step = ts->time_step;

242:   PetscFunctionBegin;
243:   PetscCall(PetscCitationsRegister(citation, &cited));

245:   if (!ts->steprollback) {
246:     if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
247:     if (th->vec_dot_prev) PetscCall(VecCopy(th->V0, th->vec_dot_prev));
248:     PetscCall(VecCopy(ts->vec_sol, th->X0));
249:     PetscCall(VecCopy(ts->vec_dot, th->V0));
250:     PetscCall(VecCopy(th->A1, th->A0));
251:   }

253:   th->status = TS_STEP_INCOMPLETE;
254:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
255:     if (ts->steprestart) {
256:       PetscCall(TSAlpha_Restart(ts, &stageok));
257:       if (!stageok) goto reject_step;
258:     }

260:     PetscCall(TSAlpha_StageTime(ts));
261:     PetscCall(TSAlpha_ApplyPredictor(ts, th->X1));
262:     PetscCall(TSPreStage(ts, th->stage_time));
263:     PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
264:     PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
265:     PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
266:     if (!stageok) goto reject_step;

268:     th->status = TS_STEP_PENDING;
269:     PetscCall(VecCopy(th->X1, ts->vec_sol));
270:     PetscCall(VecCopy(th->V1, ts->vec_dot));
271:     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
272:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
273:     if (!accept) {
274:       PetscCall(VecCopy(th->X0, ts->vec_sol));
275:       PetscCall(VecCopy(th->V0, ts->vec_dot));
276:       ts->time_step = next_time_step;
277:       goto reject_step;
278:     }

280:     ts->ptime += ts->time_step;
281:     ts->time_step = next_time_step;
282:     break;

284:   reject_step:
285:     ts->reject++;
286:     accept = PETSC_FALSE;
287:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
288:       ts->reason = TS_DIVERGED_STEP_REJECTED;
289:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
290:     }
291:   }
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
296: {
297:   TS_Alpha *th = (TS_Alpha *)ts->data;
298:   Vec       X  = th->X1;              /* X = solution */
299:   Vec       V  = th->V1;              /* V = solution */
300:   Vec       Y  = th->vec_lte_work[0]; /* Y = X + LTE  */
301:   Vec       Z  = th->vec_lte_work[1]; /* Z = V + LTE  */
302:   PetscReal enormX, enormV, enormXa, enormVa, enormXr, enormVr;

304:   PetscFunctionBegin;
305:   if (!th->vec_sol_prev) {
306:     *wlte = -1;
307:     PetscFunctionReturn(PETSC_SUCCESS);
308:   }
309:   if (!th->vec_dot_prev) {
310:     *wlte = -1;
311:     PetscFunctionReturn(PETSC_SUCCESS);
312:   }
313:   if (!th->vec_lte_work[0]) {
314:     *wlte = -1;
315:     PetscFunctionReturn(PETSC_SUCCESS);
316:   }
317:   if (!th->vec_lte_work[1]) {
318:     *wlte = -1;
319:     PetscFunctionReturn(PETSC_SUCCESS);
320:   }
321:   if (ts->steprestart) {
322:     /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
323:     PetscCall(VecAXPY(Y, 1, X));
324:     PetscCall(VecAXPY(Z, 1, V));
325:   } else {
326:     /* Compute LTE using backward differences with non-constant time step */
327:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
328:     PetscReal   a = 1 + h_prev / h;
329:     PetscScalar scal[3];
330:     Vec         vecX[3], vecV[3];
331:     scal[0] = +1 / a;
332:     scal[1] = -1 / (a - 1);
333:     scal[2] = +1 / (a * (a - 1));
334:     vecX[0] = th->X1;
335:     vecX[1] = th->X0;
336:     vecX[2] = th->vec_sol_prev;
337:     vecV[0] = th->V1;
338:     vecV[1] = th->V0;
339:     vecV[2] = th->vec_dot_prev;
340:     PetscCall(VecCopy(X, Y));
341:     PetscCall(VecMAXPY(Y, 3, scal, vecX));
342:     PetscCall(VecCopy(V, Z));
343:     PetscCall(VecMAXPY(Z, 3, scal, vecV));
344:   }
345:   /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
346:   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, &enormX, &enormXa, &enormXr));
347:   PetscCall(TSErrorWeightedNorm(ts, V, Z, wnormtype, &enormV, &enormVa, &enormVr));
348:   if (wnormtype == NORM_2) *wlte = PetscSqrtReal(PetscSqr(enormX) / 2 + PetscSqr(enormV) / 2);
349:   else *wlte = PetscMax(enormX, enormV);
350:   if (order) *order = 2;
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

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

358:   PetscFunctionBegin;
359:   PetscCall(VecCopy(th->X0, ts->vec_sol));
360:   PetscCall(VecCopy(th->V0, ts->vec_dot));
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

364: /*
365: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
366: {
367:   TS_Alpha       *th = (TS_Alpha*)ts->data;
368:   PetscReal      dt  = t - ts->ptime;

370:   PetscFunctionBegin;
371:   PetscCall(VecCopy(ts->vec_dot,V));
372:   PetscCall(VecAXPY(V,dt*(1-th->Gamma),th->A0));
373:   PetscCall(VecAXPY(V,dt*th->Gamma,th->A1));
374:   PetscCall(VecCopy(ts->vec_sol,X));
375:   PetscCall(VecAXPY(X,dt,V));
376:   PetscCall(VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0));
377:   PetscCall(VecAXPY(X,dt*dt*th->Beta,th->A1));
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }
380: */

382: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
383: {
384:   TS_Alpha *th = (TS_Alpha *)ts->data;
385:   PetscReal ta = th->stage_time;
386:   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;

388:   PetscFunctionBegin;
389:   PetscCall(TSAlpha_StageVecs(ts, X));
390:   /* F = Function(ta,Xa,Va,Aa) */
391:   PetscCall(TSComputeI2Function(ts, ta, Xa, Va, Aa, F));
392:   PetscCall(VecScale(F, th->scale_F));
393:   PetscFunctionReturn(PETSC_SUCCESS);
394: }

396: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
397: {
398:   TS_Alpha *th = (TS_Alpha *)ts->data;
399:   PetscReal ta = th->stage_time;
400:   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
401:   PetscReal dVdX = th->shift_V, dAdX = th->shift_A;

403:   PetscFunctionBegin;
404:   /* J,P = Jacobian(ta,Xa,Va,Aa) */
405:   PetscCall(TSComputeI2Jacobian(ts, ta, Xa, Va, Aa, dVdX, dAdX, J, P));
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }

409: static PetscErrorCode TSReset_Alpha(TS ts)
410: {
411:   TS_Alpha *th = (TS_Alpha *)ts->data;

413:   PetscFunctionBegin;
414:   PetscCall(VecDestroy(&th->X0));
415:   PetscCall(VecDestroy(&th->Xa));
416:   PetscCall(VecDestroy(&th->X1));
417:   PetscCall(VecDestroy(&th->V0));
418:   PetscCall(VecDestroy(&th->Va));
419:   PetscCall(VecDestroy(&th->V1));
420:   PetscCall(VecDestroy(&th->A0));
421:   PetscCall(VecDestroy(&th->Aa));
422:   PetscCall(VecDestroy(&th->A1));
423:   PetscCall(VecDestroy(&th->vec_sol_prev));
424:   PetscCall(VecDestroy(&th->vec_dot_prev));
425:   PetscCall(VecDestroy(&th->vec_lte_work[0]));
426:   PetscCall(VecDestroy(&th->vec_lte_work[1]));
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: static PetscErrorCode TSDestroy_Alpha(TS ts)
431: {
432:   PetscFunctionBegin;
433:   PetscCall(TSReset_Alpha(ts));
434:   PetscCall(PetscFree(ts->data));

436:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", NULL));
437:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", NULL));
438:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", NULL));
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: static PetscErrorCode TSSetUp_Alpha(TS ts)
443: {
444:   TS_Alpha *th = (TS_Alpha *)ts->data;
445:   PetscBool match;

447:   PetscFunctionBegin;
448:   PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
449:   PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
450:   PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
451:   PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
452:   PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
453:   PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
454:   PetscCall(VecDuplicate(ts->vec_sol, &th->A0));
455:   PetscCall(VecDuplicate(ts->vec_sol, &th->Aa));
456:   PetscCall(VecDuplicate(ts->vec_sol, &th->A1));

458:   PetscCall(TSGetAdapt(ts, &ts->adapt));
459:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
460:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
461:   if (!match) {
462:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
463:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_dot_prev));
464:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[0]));
465:     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[1]));
466:   }

468:   PetscCall(TSGetSNES(ts, &ts->snes));
469:   PetscFunctionReturn(PETSC_SUCCESS);
470: }

472: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
473: {
474:   TS_Alpha *th = (TS_Alpha *)ts->data;

476:   PetscFunctionBegin;
477:   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
478:   {
479:     PetscBool flg;
480:     PetscReal radius = 1;
481:     PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlpha2SetRadius", radius, &radius, &flg));
482:     if (flg) PetscCall(TSAlpha2SetRadius(ts, radius));
483:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlpha2SetParams", th->Alpha_m, &th->Alpha_m, NULL));
484:     PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlpha2SetParams", th->Alpha_f, &th->Alpha_f, NULL));
485:     PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlpha2SetParams", th->Gamma, &th->Gamma, NULL));
486:     PetscCall(PetscOptionsReal("-ts_alpha_beta", "Algorithmic parameter beta", "TSAlpha2SetParams", th->Beta, &th->Beta, NULL));
487:     PetscCall(TSAlpha2SetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma, th->Beta));
488:   }
489:   PetscOptionsHeadEnd();
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
494: {
495:   TS_Alpha *th = (TS_Alpha *)ts->data;
496:   PetscBool iascii;

498:   PetscFunctionBegin;
499:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
500:   if (iascii) PetscCall(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));
501:   PetscFunctionReturn(PETSC_SUCCESS);
502: }

504: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
505: {
506:   PetscReal alpha_m, alpha_f, gamma, beta;

508:   PetscFunctionBegin;
509:   PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
510:   alpha_m = (2 - radius) / (1 + radius);
511:   alpha_f = 1 / (1 + radius);
512:   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
513:   beta    = (PetscReal)0.5 * (1 + alpha_m - alpha_f);
514:   beta *= beta;
515:   PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

519: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
520: {
521:   TS_Alpha *th  = (TS_Alpha *)ts->data;
522:   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
523:   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;

525:   PetscFunctionBegin;
526:   th->Alpha_m = alpha_m;
527:   th->Alpha_f = alpha_f;
528:   th->Gamma   = gamma;
529:   th->Beta    = beta;
530:   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
535: {
536:   TS_Alpha *th = (TS_Alpha *)ts->data;

538:   PetscFunctionBegin;
539:   if (alpha_m) *alpha_m = th->Alpha_m;
540:   if (alpha_f) *alpha_f = th->Alpha_f;
541:   if (gamma) *gamma = th->Gamma;
542:   if (beta) *beta = th->Beta;
543:   PetscFunctionReturn(PETSC_SUCCESS);
544: }

546: /*MC
547:   TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method for second-order systems {cite}`chung1993`

549:   Level: beginner

551: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
552: M*/
553: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
554: {
555:   TS_Alpha *th;

557:   PetscFunctionBegin;
558:   ts->ops->reset          = TSReset_Alpha;
559:   ts->ops->destroy        = TSDestroy_Alpha;
560:   ts->ops->view           = TSView_Alpha;
561:   ts->ops->setup          = TSSetUp_Alpha;
562:   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
563:   ts->ops->step           = TSStep_Alpha;
564:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
565:   ts->ops->rollback       = TSRollBack_Alpha;
566:   /*ts->ops->interpolate  = TSInterpolate_Alpha;*/
567:   ts->ops->snesfunction  = SNESTSFormFunction_Alpha;
568:   ts->ops->snesjacobian  = SNESTSFormJacobian_Alpha;
569:   ts->default_adapt_type = TSADAPTNONE;

571:   ts->usessnes = PETSC_TRUE;

573:   PetscCall(PetscNew(&th));
574:   ts->data = (void *)th;

576:   th->Alpha_m = 0.5;
577:   th->Alpha_f = 0.5;
578:   th->Gamma   = 0.5;
579:   th->Beta    = 0.25;
580:   th->order   = 2;

582:   th->predictor     = NULL;
583:   th->predictor_ctx = NULL;

585:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", TSAlpha2SetRadius_Alpha));
586:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", TSAlpha2SetParams_Alpha));
587:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", TSAlpha2GetParams_Alpha));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

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

595:   Logically Collective

597:   Input Parameters:
598: + ts     - timestepping context
599: - radius - the desired spectral radius

601:   Options Database Key:
602: . -ts_alpha_radius <radius> - set the desired spectral radius

604:   Level: intermediate

606:   Notes:

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

612:   $$
613:   \begin{align*}
614:   \alpha_m = (2-\rho)/(1+\rho) \\
615:   \alpha_f = 1/(1+\rho)
616:   \end{align*}
617:   $$

619: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
620: @*/
621: PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
622: {
623:   PetscFunctionBegin;
626:   PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
627:   PetscTryMethod(ts, "TSAlpha2SetRadius_C", (TS, PetscReal), (ts, radius));
628:   PetscFunctionReturn(PETSC_SUCCESS);
629: }

631: /*@
632:   TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`

634:   Logically Collective

636:   Input Parameters:
637: + ts      - timestepping context
638: . alpha_m - algorithmic parameter
639: . alpha_f - algorithmic parameter
640: . gamma   - algorithmic parameter
641: - beta    - algorithmic parameter

643:   Options Database Keys:
644: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
645: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
646: . -ts_alpha_gamma   <gamma>   - set gamma
647: - -ts_alpha_beta    <beta>    - set beta

649:   Level: advanced

651:   Notes:
652:   Second-order accuracy can be obtained so long as\:

654:   $$
655:   \begin{align*}
656:   \gamma = 1/2 + \alpha_m - \alpha_f \\
657:   \beta  = 1/4 (1 + \alpha_m - \alpha_f)^2.
658:   \end{align*}
659:   $$

661:   Unconditional stability requires\:
662:   $$
663:   \alpha_m >= \alpha_f >= 1/2.
664:   $$

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

671: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
672: @*/
673: PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
674: {
675:   PetscFunctionBegin;
681:   PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
682:   PetscFunctionReturn(PETSC_SUCCESS);
683: }

685: /*@
686:   TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`

688:   Not Collective

690:   Input Parameter:
691: . ts - timestepping context

693:   Output Parameters:
694: + alpha_m - algorithmic parameter
695: . alpha_f - algorithmic parameter
696: . gamma   - algorithmic parameter
697: - beta    - algorithmic parameter

699:   Level: advanced

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

706: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
707: @*/
708: PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
709: {
710:   PetscFunctionBegin;
712:   if (alpha_m) PetscAssertPointer(alpha_m, 2);
713:   if (alpha_f) PetscAssertPointer(alpha_f, 3);
714:   if (gamma) PetscAssertPointer(gamma, 4);
715:   if (beta) PetscAssertPointer(beta, 5);
716:   PetscUseMethod(ts, "TSAlpha2GetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma, beta));
717:   PetscFunctionReturn(PETSC_SUCCESS);
718: }