Actual source code: alpha2.c
petsc-3.7.5 2017-01-01
1: /*
2: Code for timestepping with implicit generalized-\alpha method
3: for second order systems.
4: */
5: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
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: PetscBool adapt;
38: Vec vec_sol_prev;
39: Vec vec_dot_prev;
40: Vec vec_lte_work[2];
42: TSStepStatus status;
43: } TS_Alpha;
47: static PetscErrorCode TSAlpha_StageTime(TS ts)
48: {
49: TS_Alpha *th = (TS_Alpha*)ts->data;
50: PetscReal t = ts->ptime;
51: PetscReal dt = ts->time_step;
52: PetscReal Alpha_m = th->Alpha_m;
53: PetscReal Alpha_f = th->Alpha_f;
54: PetscReal Gamma = th->Gamma;
55: PetscReal Beta = th->Beta;
58: th->stage_time = t + Alpha_f*dt;
59: th->shift_V = Gamma/(dt*Beta);
60: th->shift_A = Alpha_m/(Alpha_f*dt*dt*Beta);
61: th->scale_F = 1/Alpha_f;
62: return(0);
63: }
67: static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X)
68: {
69: TS_Alpha *th = (TS_Alpha*)ts->data;
70: Vec X1 = X, V1 = th->V1, A1 = th->A1;
71: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
72: Vec X0 = th->X0, V0 = th->V0, A0 = th->A0;
73: PetscReal dt = ts->time_step;
74: PetscReal Alpha_m = th->Alpha_m;
75: PetscReal Alpha_f = th->Alpha_f;
76: PetscReal Gamma = th->Gamma;
77: PetscReal Beta = th->Beta;
81: /* A1 = ... */
82: VecWAXPY(A1,-1.0,X0,X1);
83: VecAXPY (A1,-dt,V0);
84: VecAXPBY(A1,-(1-2*Beta)/(2*Beta),1/(dt*dt*Beta),A0);
85: /* V1 = ... */
86: VecWAXPY(V1,(1.0-Gamma)/Gamma,A0,A1);
87: VecAYPX (V1,dt*Gamma,V0);
88: /* Xa = X0 + Alpha_f*(X1-X0) */
89: VecWAXPY(Xa,-1.0,X0,X1);
90: VecAYPX (Xa,Alpha_f,X0);
91: /* Va = V0 + Alpha_f*(V1-V0) */
92: VecWAXPY(Va,-1.0,V0,V1);
93: VecAYPX (Va,Alpha_f,V0);
94: /* Aa = A0 + Alpha_m*(A1-A0) */
95: VecWAXPY(Aa,-1.0,A0,A1);
96: VecAYPX (Aa,Alpha_m,A0);
97: return(0);
98: }
102: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
103: {
104: PetscInt nits,lits;
108: SNESSolve(ts->snes,b,x);
109: SNESGetIterationNumber(ts->snes,&nits);
110: SNESGetLinearSolveIterations(ts->snes,&lits);
111: ts->snes_its += nits; ts->ksp_its += lits;
112: return(0);
113: }
115: /*
116: Compute a consistent initial state for the generalized-alpha method.
117: - Solve two successive backward Euler steps with halved time step.
118: - Compute the initial second time derivative using backward differences.
119: - If using adaptivity, estimate the LTE of the initial step.
120: */
123: static PetscErrorCode TSAlpha_Restart(TS ts,PetscBool *initok)
124: {
125: TS_Alpha *th = (TS_Alpha*)ts->data;
126: PetscReal time_step;
127: PetscReal alpha_m,alpha_f,gamma,beta;
128: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
129: Vec V0 = ts->vec_dot, V1, V2 = th->V1;
130: PetscBool stageok;
134: VecDuplicate(X0,&X1);
135: VecDuplicate(V0,&V1);
137: /* Setup backward Euler with halved time step */
138: TSAlpha2GetParams(ts,&alpha_m,&alpha_f,&gamma,&beta);
139: TSAlpha2SetParams(ts,1,1,1,0.5);
140: TSGetTimeStep(ts,&time_step);
141: ts->time_step = time_step/2;
142: TSAlpha_StageTime(ts);
143: th->stage_time = ts->ptime;
144: VecZeroEntries(th->A0);
146: /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
147: th->stage_time += ts->time_step;
148: VecCopy(X0,th->X0);
149: VecCopy(V0,th->V0);
150: TSPreStage(ts,th->stage_time);
151: VecCopy(th->X0,X1);
152: TS_SNESSolve(ts,NULL,X1);
153: VecCopy(th->V1,V1);
154: TSPostStage(ts,th->stage_time,0,&X1);
155: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
156: if (!stageok) goto finally;
158: /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
159: th->stage_time += ts->time_step;
160: VecCopy(X1,th->X0);
161: VecCopy(V1,th->V0);
162: TSPreStage(ts,th->stage_time);
163: VecCopy(th->X0,X2);
164: TS_SNESSolve(ts,NULL,X2);
165: VecCopy(th->V1,V2);
166: TSPostStage(ts,th->stage_time,0,&X2);
167: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
168: if (!stageok) goto finally;
170: /* Compute A0 ~ dV/dt at t0 with backward differences */
171: VecZeroEntries(th->A0);
172: VecAXPY(th->A0,-3/ts->time_step,V0);
173: VecAXPY(th->A0,+4/ts->time_step,V1);
174: VecAXPY(th->A0,-1/ts->time_step,V2);
176: /* Rough, lower-order estimate LTE of the initial step */
177: if (th->adapt) {
178: VecZeroEntries(th->vec_lte_work[0]);
179: VecAXPY(th->vec_lte_work[0],+2,X2);
180: VecAXPY(th->vec_lte_work[0],-4,X1);
181: VecAXPY(th->vec_lte_work[0],+2,X0);
182: }
183: if (th->adapt) {
184: VecZeroEntries(th->vec_lte_work[1]);
185: VecAXPY(th->vec_lte_work[1],+2,V2);
186: VecAXPY(th->vec_lte_work[1],-4,V1);
187: VecAXPY(th->vec_lte_work[1],+2,V0);
188: }
190: finally:
191: /* Revert TSAlpha to the initial state (t0,X0,V0) */
192: if (initok) *initok = stageok;
193: TSSetTimeStep(ts,time_step);
194: TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
195: VecCopy(ts->vec_sol,th->X0);
196: VecCopy(ts->vec_dot,th->V0);
198: VecDestroy(&X1);
199: VecDestroy(&V1);
200: return(0);
201: }
205: static PetscErrorCode TSStep_Alpha(TS ts)
206: {
207: TS_Alpha *th = (TS_Alpha*)ts->data;
208: PetscInt rejections = 0;
209: PetscBool stageok,accept = PETSC_TRUE;
210: PetscReal next_time_step = ts->time_step;
214: PetscCitationsRegister(citation,&cited);
216: if (!ts->steprollback) {
217: if (th->adapt) { VecCopy(th->X0,th->vec_sol_prev); }
218: if (th->adapt) { VecCopy(th->V0,th->vec_dot_prev); }
219: VecCopy(ts->vec_sol,th->X0);
220: VecCopy(ts->vec_dot,th->V0);
221: VecCopy(th->A1,th->A0);
222: }
224: th->status = TS_STEP_INCOMPLETE;
225: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
227: if (ts->steprestart) {
228: TSAlpha_Restart(ts,&stageok);
229: if (!stageok) goto reject_step;
230: }
232: TSAlpha_StageTime(ts);
233: VecCopy(th->X0,th->X1);
234: TSPreStage(ts,th->stage_time);
235: TS_SNESSolve(ts,NULL,th->X1);
236: TSPostStage(ts,th->stage_time,0,&th->Xa);
237: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
238: if (!stageok) goto reject_step;
240: th->status = TS_STEP_PENDING;
241: VecCopy(th->X1,ts->vec_sol);
242: VecCopy(th->V1,ts->vec_dot);
243: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
244: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
245: if (!accept) {
246: VecCopy(th->X0,ts->vec_sol);
247: VecCopy(th->V0,ts->vec_dot);
248: ts->time_step = next_time_step;
249: goto reject_step;
250: }
252: ts->ptime += ts->time_step;
253: ts->time_step = next_time_step;
254: break;
256: reject_step:
257: ts->reject++; accept = PETSC_FALSE;
258: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
259: ts->reason = TS_DIVERGED_STEP_REJECTED;
260: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
261: }
263: }
264: return(0);
265: }
269: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
270: {
271: TS_Alpha *th = (TS_Alpha*)ts->data;
272: Vec X = th->X1; /* X = solution */
273: Vec V = th->V1; /* V = solution */
274: Vec Y = th->vec_lte_work[0]; /* Y = X + LTE */
275: Vec Z = th->vec_lte_work[1]; /* Z = V + LTE */
276: PetscReal enormX,enormV;
280: if (ts->steprestart) {
281: /* th->vec_{sol|dot}_prev is set to the LTE in TSAlpha_Restart() */
282: VecAXPY(Y,1,X);
283: VecAXPY(Z,1,V);
284: } else {
285: /* Compute LTE using backward differences with non-constant time step */
286: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
287: PetscReal a = 1 + h_prev/h;
288: PetscScalar scal[3]; Vec vecX[3],vecV[3];
289: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
290: vecX[0] = th->X1; vecX[1] = th->X0; vecX[2] = th->vec_sol_prev;
291: vecV[0] = th->V1; vecV[1] = th->V0; vecV[2] = th->vec_dot_prev;
292: VecCopy(X,Y);
293: VecMAXPY(Y,3,scal,vecX);
294: VecCopy(V,Z);
295: VecMAXPY(Z,3,scal,vecV);
296: }
297: /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
298: TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX);
299: TSErrorWeightedNorm(ts,V,Z,wnormtype,&enormV);
300: if (wnormtype == NORM_2)
301: *wlte = PetscSqrtReal(PetscSqr(enormX)/2 + PetscSqr(enormV)/2);
302: else
303: *wlte = PetscMax(enormX,enormV);
304: if (order) *order = 2;
305: return(0);
306: }
310: static PetscErrorCode TSRollBack_Alpha(TS ts)
311: {
312: TS_Alpha *th = (TS_Alpha*)ts->data;
316: VecCopy(th->X0,ts->vec_sol);
317: VecCopy(th->V0,ts->vec_dot);
318: return(0);
319: }
321: /*
324: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
325: {
326: TS_Alpha *th = (TS_Alpha*)ts->data;
327: PetscReal dt = t - ts->ptime;
331: VecCopy(ts->vec_dot,V);
332: VecAXPY(V,dt*(1-th->Gamma),th->A0);
333: VecAXPY(V,dt*th->Gamma,th->A1);
334: VecCopy(ts->vec_sol,X);
335: VecAXPY(X,dt,V);
336: VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
337: VecAXPY(X,dt*dt*th->Beta,th->A1);
338: return(0);
339: }
340: */
344: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
345: {
346: TS_Alpha *th = (TS_Alpha*)ts->data;
347: PetscReal ta = th->stage_time;
348: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
352: TSAlpha_StageVecs(ts,X);
353: /* F = Function(ta,Xa,Va,Aa) */
354: TSComputeI2Function(ts,ta,Xa,Va,Aa,F);
355: VecScale(F,th->scale_F);
356: return(0);
357: }
361: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)
362: {
363: TS_Alpha *th = (TS_Alpha*)ts->data;
364: PetscReal ta = th->stage_time;
365: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
366: PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
370: /* J,P = Jacobian(ta,Xa,Va,Aa) */
371: TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P);
372: return(0);
373: }
377: static PetscErrorCode TSReset_Alpha(TS ts)
378: {
379: TS_Alpha *th = (TS_Alpha*)ts->data;
383: VecDestroy(&th->X0);
384: VecDestroy(&th->Xa);
385: VecDestroy(&th->X1);
386: VecDestroy(&th->V0);
387: VecDestroy(&th->Va);
388: VecDestroy(&th->V1);
389: VecDestroy(&th->A0);
390: VecDestroy(&th->Aa);
391: VecDestroy(&th->A1);
392: VecDestroy(&th->vec_sol_prev);
393: VecDestroy(&th->vec_dot_prev);
394: VecDestroy(&th->vec_lte_work[0]);
395: VecDestroy(&th->vec_lte_work[1]);
396: return(0);
397: }
401: static PetscErrorCode TSDestroy_Alpha(TS ts)
402: {
406: TSReset_Alpha(ts);
407: PetscFree(ts->data);
409: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2UseAdapt_C",NULL);
410: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL);
411: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL);
412: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL);
413: return(0);
414: }
418: static PetscErrorCode TSSetUp_Alpha(TS ts)
419: {
420: TS_Alpha *th = (TS_Alpha*)ts->data;
424: VecDuplicate(ts->vec_sol,&th->X0);
425: VecDuplicate(ts->vec_sol,&th->Xa);
426: VecDuplicate(ts->vec_sol,&th->X1);
427: VecDuplicate(ts->vec_sol,&th->V0);
428: VecDuplicate(ts->vec_sol,&th->Va);
429: VecDuplicate(ts->vec_sol,&th->V1);
430: VecDuplicate(ts->vec_sol,&th->A0);
431: VecDuplicate(ts->vec_sol,&th->Aa);
432: VecDuplicate(ts->vec_sol,&th->A1);
434: TSGetAdapt(ts,&ts->adapt);
435: TSAdaptCandidatesClear(ts->adapt);
436: if (!th->adapt) {
437: TSAdaptSetType(ts->adapt,TSADAPTNONE);
438: } else {
439: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
440: VecDuplicate(ts->vec_sol,&th->vec_dot_prev);
441: VecDuplicate(ts->vec_sol,&th->vec_lte_work[0]);
442: VecDuplicate(ts->vec_sol,&th->vec_lte_work[1]);
443: if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
444: ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
445: }
447: TSGetSNES(ts,&ts->snes);
448: return(0);
449: }
453: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
454: {
455: TS_Alpha *th = (TS_Alpha*)ts->data;
459: PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
460: {
461: PetscBool flg;
462: PetscReal radius = 1;
463: PetscBool adapt = th->adapt;
464: PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg);
465: if (flg) {TSAlpha2SetRadius(ts,radius);}
466: PetscOptionsReal("-ts_alpha_alpha_m","Algoritmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL);
467: PetscOptionsReal("-ts_alpha_alpha_f","Algoritmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL);
468: PetscOptionsReal("-ts_alpha_gamma","Algoritmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL);
469: PetscOptionsReal("-ts_alpha_beta","Algoritmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL);
470: TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta);
471: PetscOptionsBool("-ts_alpha_adapt","Use time-step adaptivity with the Alpha method","TSAlpha2UseAdapt",adapt,&adapt,&flg);
472: if (flg) {TSAlpha2UseAdapt(ts,adapt);}
473: }
474: PetscOptionsTail();
475: return(0);
476: }
480: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
481: {
482: TS_Alpha *th = (TS_Alpha*)ts->data;
483: PetscBool iascii;
487: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
488: if (iascii) {
489: 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);
490: }
491: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
492: if (ts->snes) {SNESView(ts->snes,viewer);}
493: return(0);
494: }
498: static PetscErrorCode TSAlpha2UseAdapt_Alpha(TS ts,PetscBool use)
499: {
500: TS_Alpha *th = (TS_Alpha*)ts->data;
503: if (use == th->adapt) return(0);
504: if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
505: th->adapt = use;
506: return(0);
507: }
511: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)
512: {
513: PetscReal alpha_m,alpha_f,gamma,beta;
517: if (radius < 0 || radius > 1) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
518: alpha_m = (2-radius)/(1+radius);
519: alpha_f = 1/(1+radius);
520: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
521: beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta;
522: TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
523: return(0);
524: }
528: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
529: {
530: TS_Alpha *th = (TS_Alpha*)ts->data;
531: PetscReal tol = 100*PETSC_MACHINE_EPSILON;
532: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
535: th->Alpha_m = alpha_m;
536: th->Alpha_f = alpha_f;
537: th->Gamma = gamma;
538: th->Beta = beta;
539: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
540: return(0);
541: }
545: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
546: {
547: TS_Alpha *th = (TS_Alpha*)ts->data;
550: if (alpha_m) *alpha_m = th->Alpha_m;
551: if (alpha_f) *alpha_f = th->Alpha_f;
552: if (gamma) *gamma = th->Gamma;
553: if (beta) *beta = th->Beta;
554: return(0);
555: }
557: /*MC
558: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
559: for second-order systems
561: Level: beginner
563: References:
564: J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
565: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
566: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
568: .seealso: TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams()
569: M*/
572: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
573: {
574: TS_Alpha *th;
578: ts->ops->reset = TSReset_Alpha;
579: ts->ops->destroy = TSDestroy_Alpha;
580: ts->ops->view = TSView_Alpha;
581: ts->ops->setup = TSSetUp_Alpha;
582: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
583: ts->ops->step = TSStep_Alpha;
584: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
585: ts->ops->rollback = TSRollBack_Alpha;
586: /*ts->ops->interpolate = TSInterpolate_Alpha;*/
587: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
588: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
590: PetscNewLog(ts,&th);
591: ts->data = (void*)th;
593: th->Alpha_m = 0.5;
594: th->Alpha_f = 0.5;
595: th->Gamma = 0.5;
596: th->Beta = 0.25;
597: th->order = 2;
599: th->adapt = PETSC_FALSE;
601: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2UseAdapt_C",TSAlpha2UseAdapt_Alpha);
602: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha);
603: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha);
604: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha);
605: return(0);
606: }
610: /*@
611: TSAlpha2UseAdapt - Use time-step adaptivity with the Alpha method
613: Logically Collective on TS
615: Input Parameter:
616: + ts - timestepping context
617: - use - flag to use adaptivity
619: Options Database:
620: . -ts_alpha_adapt
622: Level: intermediate
624: .seealso: TSAdapt, TSADAPTBASIC
625: @*/
626: PetscErrorCode TSAlpha2UseAdapt(TS ts,PetscBool use)
627: {
633: PetscTryMethod(ts,"TSAlpha2UseAdapt_C",(TS,PetscBool),(ts,use));
634: return(0);
635: }
639: /*@
640: TSAlpha2SetRadius - sets the desired spectral radius of the method
641: (i.e. high-frequency numerical damping)
643: Logically Collective on TS
645: The algorithmic parameters \alpha_m and \alpha_f of the
646: generalized-\alpha method can be computed in terms of a specified
647: spectral radius \rho in [0,1] for infinite time step in order to
648: control high-frequency numerical damping:
649: \alpha_m = (2-\rho)/(1+\rho)
650: \alpha_f = 1/(1+\rho)
652: Input Parameter:
653: + ts - timestepping context
654: - radius - the desired spectral radius
656: Options Database:
657: . -ts_alpha_radius <radius>
659: Level: intermediate
661: .seealso: TSAlpha2SetParams(), TSAlpha2GetParams()
662: @*/
663: PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius)
664: {
670: if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
671: PetscTryMethod(ts,"TSAlpha2SetRadius_C",(TS,PetscReal),(ts,radius));
672: return(0);
673: }
677: /*@
678: TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2
680: Logically Collective on TS
682: Second-order accuracy can be obtained so long as:
683: \gamma = 1/2 + alpha_m - alpha_f
684: \beta = 1/4 (1 + alpha_m - alpha_f)^2
686: Unconditional stability requires:
687: \alpha_m >= \alpha_f >= 1/2
690: Input Parameter:
691: + ts - timestepping context
692: . \alpha_m - algorithmic paramenter
693: . \alpha_f - algorithmic paramenter
694: . \gamma - algorithmic paramenter
695: - \beta - algorithmic paramenter
697: Options Database:
698: + -ts_alpha_alpha_m <alpha_m>
699: . -ts_alpha_alpha_f <alpha_f>
700: . -ts_alpha_gamma <gamma>
701: - -ts_alpha_beta <beta>
703: Note:
704: Use of this function is normally only required to hack TSALPHA2 to
705: use a modified integration scheme. Users should call
706: TSAlpha2SetRadius() to set the desired spectral radius of the methods
707: (i.e. high-frequency damping) in order so select optimal values for
708: these parameters.
710: Level: advanced
712: .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams()
713: @*/
714: PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
715: {
724: PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta));
725: return(0);
726: }
730: /*@
731: TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2
733: Not Collective
735: Input Parameter:
736: . ts - timestepping context
738: Output Parameters:
739: + \alpha_m - algorithmic parameter
740: . \alpha_f - algorithmic parameter
741: . \gamma - algorithmic parameter
742: - \beta - algorithmic parameter
744: Note:
745: Use of this function is normally only required to hack TSALPHA2 to
746: use a modified integration scheme. Users should call
747: TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral
748: radius of the method) in order so select optimal values for these
749: parameters.
751: Level: advanced
753: .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams()
754: @*/
755: PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
756: {
765: PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta));
766: return(0);
767: }