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[] =
9: "@article{Chung1993,\n"
10: " title = {A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation: The Generalized-$\\alpha$ Method},\n"
11: " author = {J. Chung, G. M. Hubert},\n"
12: " journal = {ASME Journal of Applied Mechanics},\n"
13: " volume = {60},\n"
14: " number = {2},\n"
15: " pages = {371--375},\n"
16: " year = {1993},\n"
17: " issn = {0021-8936},\n"
18: " doi = {http://dx.doi.org/10.1115/1.2900803}\n}\n";
20: typedef struct {
21: PetscReal stage_time;
22: PetscReal shift_V;
23: PetscReal shift_A;
24: PetscReal scale_F;
25: Vec X0,Xa,X1;
26: Vec V0,Va,V1;
27: Vec A0,Aa,A1;
29: Vec vec_dot;
31: PetscReal Alpha_m;
32: PetscReal Alpha_f;
33: PetscReal Gamma;
34: PetscReal Beta;
35: PetscInt order;
37: Vec vec_sol_prev;
38: Vec vec_dot_prev;
39: Vec vec_lte_work[2];
41: TSStepStatus status;
42: } TS_Alpha;
44: static PetscErrorCode TSAlpha_StageTime(TS ts)
45: {
46: TS_Alpha *th = (TS_Alpha*)ts->data;
47: PetscReal t = ts->ptime;
48: PetscReal dt = ts->time_step;
49: PetscReal Alpha_m = th->Alpha_m;
50: PetscReal Alpha_f = th->Alpha_f;
51: PetscReal Gamma = th->Gamma;
52: PetscReal Beta = th->Beta;
54: th->stage_time = t + Alpha_f*dt;
55: th->shift_V = Gamma/(dt*Beta);
56: th->shift_A = Alpha_m/(Alpha_f*dt*dt*Beta);
57: th->scale_F = 1/Alpha_f;
58: return 0;
59: }
61: static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X)
62: {
63: TS_Alpha *th = (TS_Alpha*)ts->data;
64: Vec X1 = X, V1 = th->V1, A1 = th->A1;
65: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
66: Vec X0 = th->X0, V0 = th->V0, A0 = th->A0;
67: PetscReal dt = ts->time_step;
68: PetscReal Alpha_m = th->Alpha_m;
69: PetscReal Alpha_f = th->Alpha_f;
70: PetscReal Gamma = th->Gamma;
71: PetscReal Beta = th->Beta;
73: /* A1 = ... */
74: VecWAXPY(A1,-1.0,X0,X1);
75: VecAXPY (A1,-dt,V0);
76: VecAXPBY(A1,-(1-2*Beta)/(2*Beta),1/(dt*dt*Beta),A0);
77: /* V1 = ... */
78: VecWAXPY(V1,(1.0-Gamma)/Gamma,A0,A1);
79: VecAYPX (V1,dt*Gamma,V0);
80: /* Xa = X0 + Alpha_f*(X1-X0) */
81: VecWAXPY(Xa,-1.0,X0,X1);
82: VecAYPX (Xa,Alpha_f,X0);
83: /* Va = V0 + Alpha_f*(V1-V0) */
84: VecWAXPY(Va,-1.0,V0,V1);
85: VecAYPX (Va,Alpha_f,V0);
86: /* Aa = A0 + Alpha_m*(A1-A0) */
87: VecWAXPY(Aa,-1.0,A0,A1);
88: VecAYPX (Aa,Alpha_m,A0);
89: return 0;
90: }
92: static PetscErrorCode TSAlpha_SNESSolve(TS ts,Vec b,Vec x)
93: {
94: PetscInt nits,lits;
96: SNESSolve(ts->snes,b,x);
97: SNESGetIterationNumber(ts->snes,&nits);
98: SNESGetLinearSolveIterations(ts->snes,&lits);
99: ts->snes_its += nits; 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) {
207: if (ts->steprestart) {
208: TSAlpha_Restart(ts,&stageok);
209: if (!stageok) goto reject_step;
210: }
212: TSAlpha_StageTime(ts);
213: VecCopy(th->X0,th->X1);
214: TSPreStage(ts,th->stage_time);
215: TSAlpha_SNESSolve(ts,NULL,th->X1);
216: TSPostStage(ts,th->stage_time,0,&th->Xa);
217: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
218: if (!stageok) goto reject_step;
220: th->status = TS_STEP_PENDING;
221: VecCopy(th->X1,ts->vec_sol);
222: VecCopy(th->V1,ts->vec_dot);
223: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
224: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
225: if (!accept) {
226: VecCopy(th->X0,ts->vec_sol);
227: VecCopy(th->V0,ts->vec_dot);
228: ts->time_step = next_time_step;
229: goto reject_step;
230: }
232: ts->ptime += ts->time_step;
233: ts->time_step = next_time_step;
234: break;
236: reject_step:
237: ts->reject++; 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=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
241: }
243: }
244: return 0;
245: }
247: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
248: {
249: TS_Alpha *th = (TS_Alpha*)ts->data;
250: Vec X = th->X1; /* X = solution */
251: Vec V = th->V1; /* V = solution */
252: Vec Y = th->vec_lte_work[0]; /* Y = X + LTE */
253: Vec Z = th->vec_lte_work[1]; /* Z = V + LTE */
254: PetscReal enormX,enormV,enormXa,enormVa,enormXr,enormVr;
256: if (!th->vec_sol_prev) {*wlte = -1; return 0;}
257: if (!th->vec_dot_prev) {*wlte = -1; return 0;}
258: if (!th->vec_lte_work[0]) {*wlte = -1; return 0;}
259: if (!th->vec_lte_work[1]) {*wlte = -1; return 0;}
260: if (ts->steprestart) {
261: /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
262: VecAXPY(Y,1,X);
263: VecAXPY(Z,1,V);
264: } else {
265: /* Compute LTE using backward differences with non-constant time step */
266: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
267: PetscReal a = 1 + h_prev/h;
268: PetscScalar scal[3]; Vec vecX[3],vecV[3];
269: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
270: vecX[0] = th->X1; vecX[1] = th->X0; vecX[2] = th->vec_sol_prev;
271: vecV[0] = th->V1; vecV[1] = th->V0; vecV[2] = th->vec_dot_prev;
272: VecCopy(X,Y);
273: VecMAXPY(Y,3,scal,vecX);
274: VecCopy(V,Z);
275: VecMAXPY(Z,3,scal,vecV);
276: }
277: /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
278: TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX,&enormXa,&enormXr);
279: TSErrorWeightedNorm(ts,V,Z,wnormtype,&enormV,&enormVa,&enormVr);
280: if (wnormtype == NORM_2)
281: *wlte = PetscSqrtReal(PetscSqr(enormX)/2 + PetscSqr(enormV)/2);
282: else
283: *wlte = PetscMax(enormX,enormV);
284: if (order) *order = 2;
285: return 0;
286: }
288: static PetscErrorCode TSRollBack_Alpha(TS ts)
289: {
290: TS_Alpha *th = (TS_Alpha*)ts->data;
292: VecCopy(th->X0,ts->vec_sol);
293: VecCopy(th->V0,ts->vec_dot);
294: return 0;
295: }
297: /*
298: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
299: {
300: TS_Alpha *th = (TS_Alpha*)ts->data;
301: PetscReal dt = t - ts->ptime;
303: VecCopy(ts->vec_dot,V);
304: VecAXPY(V,dt*(1-th->Gamma),th->A0);
305: VecAXPY(V,dt*th->Gamma,th->A1);
306: VecCopy(ts->vec_sol,X);
307: VecAXPY(X,dt,V);
308: VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
309: VecAXPY(X,dt*dt*th->Beta,th->A1);
310: return 0;
311: }
312: */
314: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
315: {
316: TS_Alpha *th = (TS_Alpha*)ts->data;
317: PetscReal ta = th->stage_time;
318: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
320: TSAlpha_StageVecs(ts,X);
321: /* F = Function(ta,Xa,Va,Aa) */
322: TSComputeI2Function(ts,ta,Xa,Va,Aa,F);
323: VecScale(F,th->scale_F);
324: return 0;
325: }
327: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)
328: {
329: TS_Alpha *th = (TS_Alpha*)ts->data;
330: PetscReal ta = th->stage_time;
331: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
332: PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
334: /* J,P = Jacobian(ta,Xa,Va,Aa) */
335: TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P);
336: return 0;
337: }
339: static PetscErrorCode TSReset_Alpha(TS ts)
340: {
341: TS_Alpha *th = (TS_Alpha*)ts->data;
343: VecDestroy(&th->X0);
344: VecDestroy(&th->Xa);
345: VecDestroy(&th->X1);
346: VecDestroy(&th->V0);
347: VecDestroy(&th->Va);
348: VecDestroy(&th->V1);
349: VecDestroy(&th->A0);
350: VecDestroy(&th->Aa);
351: VecDestroy(&th->A1);
352: VecDestroy(&th->vec_sol_prev);
353: VecDestroy(&th->vec_dot_prev);
354: VecDestroy(&th->vec_lte_work[0]);
355: VecDestroy(&th->vec_lte_work[1]);
356: return 0;
357: }
359: static PetscErrorCode TSDestroy_Alpha(TS ts)
360: {
361: TSReset_Alpha(ts);
362: PetscFree(ts->data);
364: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL);
365: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL);
366: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL);
367: return 0;
368: }
370: static PetscErrorCode TSSetUp_Alpha(TS ts)
371: {
372: TS_Alpha *th = (TS_Alpha*)ts->data;
373: PetscBool match;
375: VecDuplicate(ts->vec_sol,&th->X0);
376: VecDuplicate(ts->vec_sol,&th->Xa);
377: VecDuplicate(ts->vec_sol,&th->X1);
378: VecDuplicate(ts->vec_sol,&th->V0);
379: VecDuplicate(ts->vec_sol,&th->Va);
380: VecDuplicate(ts->vec_sol,&th->V1);
381: VecDuplicate(ts->vec_sol,&th->A0);
382: VecDuplicate(ts->vec_sol,&th->Aa);
383: VecDuplicate(ts->vec_sol,&th->A1);
385: TSGetAdapt(ts,&ts->adapt);
386: TSAdaptCandidatesClear(ts->adapt);
387: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
388: if (!match) {
389: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
390: VecDuplicate(ts->vec_sol,&th->vec_dot_prev);
391: VecDuplicate(ts->vec_sol,&th->vec_lte_work[0]);
392: VecDuplicate(ts->vec_sol,&th->vec_lte_work[1]);
393: }
395: TSGetSNES(ts,&ts->snes);
396: return 0;
397: }
399: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
400: {
401: TS_Alpha *th = (TS_Alpha*)ts->data;
403: PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
404: {
405: PetscBool flg;
406: PetscReal radius = 1;
407: PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg);
408: if (flg) TSAlpha2SetRadius(ts,radius);
409: PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL);
410: PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL);
411: PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL);
412: PetscOptionsReal("-ts_alpha_beta","Algorithmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL);
413: TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta);
414: }
415: PetscOptionsTail();
416: return 0;
417: }
419: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
420: {
421: TS_Alpha *th = (TS_Alpha*)ts->data;
422: PetscBool iascii;
424: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
425: if (iascii) {
426: 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);
427: }
428: return 0;
429: }
431: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)
432: {
433: PetscReal alpha_m,alpha_f,gamma,beta;
436: alpha_m = (2-radius)/(1+radius);
437: alpha_f = 1/(1+radius);
438: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
439: beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta;
440: TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
441: return 0;
442: }
444: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
445: {
446: TS_Alpha *th = (TS_Alpha*)ts->data;
447: PetscReal tol = 100*PETSC_MACHINE_EPSILON;
448: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
450: th->Alpha_m = alpha_m;
451: th->Alpha_f = alpha_f;
452: th->Gamma = gamma;
453: th->Beta = beta;
454: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
455: return 0;
456: }
458: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
459: {
460: TS_Alpha *th = (TS_Alpha*)ts->data;
462: if (alpha_m) *alpha_m = th->Alpha_m;
463: if (alpha_f) *alpha_f = th->Alpha_f;
464: if (gamma) *gamma = th->Gamma;
465: if (beta) *beta = th->Beta;
466: return 0;
467: }
469: /*MC
470: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
471: for second-order systems
473: Level: beginner
475: References:
476: . * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
477: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
478: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
480: .seealso: TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams()
481: M*/
482: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
483: {
484: TS_Alpha *th;
486: ts->ops->reset = TSReset_Alpha;
487: ts->ops->destroy = TSDestroy_Alpha;
488: ts->ops->view = TSView_Alpha;
489: ts->ops->setup = TSSetUp_Alpha;
490: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
491: ts->ops->step = TSStep_Alpha;
492: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
493: ts->ops->rollback = TSRollBack_Alpha;
494: /*ts->ops->interpolate = TSInterpolate_Alpha;*/
495: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
496: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
497: ts->default_adapt_type = TSADAPTNONE;
499: ts->usessnes = PETSC_TRUE;
501: PetscNewLog(ts,&th);
502: ts->data = (void*)th;
504: th->Alpha_m = 0.5;
505: th->Alpha_f = 0.5;
506: th->Gamma = 0.5;
507: th->Beta = 0.25;
508: th->order = 2;
510: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha);
511: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha);
512: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha);
513: return 0;
514: }
516: /*@
517: TSAlpha2SetRadius - sets the desired spectral radius of the method
518: (i.e. high-frequency numerical damping)
520: Logically Collective on TS
522: The algorithmic parameters \alpha_m and \alpha_f of the
523: generalized-\alpha method can be computed in terms of a specified
524: spectral radius \rho in [0,1] for infinite time step in order to
525: control high-frequency numerical damping:
526: \alpha_m = (2-\rho)/(1+\rho)
527: \alpha_f = 1/(1+\rho)
529: Input Parameters:
530: + ts - timestepping context
531: - radius - the desired spectral radius
533: Options Database:
534: . -ts_alpha_radius <radius> - set the desired spectral radius
536: Level: intermediate
538: .seealso: TSAlpha2SetParams(), TSAlpha2GetParams()
539: @*/
540: PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius)
541: {
545: PetscTryMethod(ts,"TSAlpha2SetRadius_C",(TS,PetscReal),(ts,radius));
546: return 0;
547: }
549: /*@
550: TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2
552: Logically Collective on TS
554: Second-order accuracy can be obtained so long as:
555: \gamma = 1/2 + alpha_m - alpha_f
556: \beta = 1/4 (1 + alpha_m - alpha_f)^2
558: Unconditional stability requires:
559: \alpha_m >= \alpha_f >= 1/2
561: Input Parameters:
562: + ts - timestepping context
563: . \alpha_m - algorithmic parameter
564: . \alpha_f - algorithmic parameter
565: . \gamma - algorithmic parameter
566: - \beta - algorithmic parameter
568: Options Database:
569: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
570: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
571: . -ts_alpha_gamma <gamma> - set gamma
572: - -ts_alpha_beta <beta> - set beta
574: Note:
575: Use of this function is normally only required to hack TSALPHA2 to
576: use a modified integration scheme. Users should call
577: TSAlpha2SetRadius() to set the desired spectral radius of the methods
578: (i.e. high-frequency damping) in order so select optimal values for
579: these parameters.
581: Level: advanced
583: .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams()
584: @*/
585: PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
586: {
592: PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta));
593: return 0;
594: }
596: /*@
597: TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2
599: Not Collective
601: Input Parameter:
602: . ts - timestepping context
604: Output Parameters:
605: + \alpha_m - algorithmic parameter
606: . \alpha_f - algorithmic parameter
607: . \gamma - algorithmic parameter
608: - \beta - algorithmic parameter
610: Note:
611: Use of this function is normally only required to hack TSALPHA2 to
612: use a modified integration scheme. Users should call
613: TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral
614: radius of the method) in order so select optimal values for these
615: parameters.
617: Level: advanced
619: .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams()
620: @*/
621: PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
622: {
628: PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta));
629: return 0;
630: }