Actual source code: alpha2.c
petsc-3.14.6 2021-03-30
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;
55: th->stage_time = t + Alpha_f*dt;
56: th->shift_V = Gamma/(dt*Beta);
57: th->shift_A = Alpha_m/(Alpha_f*dt*dt*Beta);
58: th->scale_F = 1/Alpha_f;
59: return(0);
60: }
62: static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X)
63: {
64: TS_Alpha *th = (TS_Alpha*)ts->data;
65: Vec X1 = X, V1 = th->V1, A1 = th->A1;
66: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
67: Vec X0 = th->X0, V0 = th->V0, A0 = th->A0;
68: PetscReal dt = ts->time_step;
69: PetscReal Alpha_m = th->Alpha_m;
70: PetscReal Alpha_f = th->Alpha_f;
71: PetscReal Gamma = th->Gamma;
72: PetscReal Beta = th->Beta;
76: /* A1 = ... */
77: VecWAXPY(A1,-1.0,X0,X1);
78: VecAXPY (A1,-dt,V0);
79: VecAXPBY(A1,-(1-2*Beta)/(2*Beta),1/(dt*dt*Beta),A0);
80: /* V1 = ... */
81: VecWAXPY(V1,(1.0-Gamma)/Gamma,A0,A1);
82: VecAYPX (V1,dt*Gamma,V0);
83: /* Xa = X0 + Alpha_f*(X1-X0) */
84: VecWAXPY(Xa,-1.0,X0,X1);
85: VecAYPX (Xa,Alpha_f,X0);
86: /* Va = V0 + Alpha_f*(V1-V0) */
87: VecWAXPY(Va,-1.0,V0,V1);
88: VecAYPX (Va,Alpha_f,V0);
89: /* Aa = A0 + Alpha_m*(A1-A0) */
90: VecWAXPY(Aa,-1.0,A0,A1);
91: VecAYPX (Aa,Alpha_m,A0);
92: return(0);
93: }
95: static PetscErrorCode TSAlpha_SNESSolve(TS ts,Vec b,Vec x)
96: {
97: PetscInt nits,lits;
101: SNESSolve(ts->snes,b,x);
102: SNESGetIterationNumber(ts->snes,&nits);
103: SNESGetLinearSolveIterations(ts->snes,&lits);
104: ts->snes_its += nits; ts->ksp_its += lits;
105: return(0);
106: }
108: /*
109: Compute a consistent initial state for the generalized-alpha method.
110: - Solve two successive backward Euler steps with halved time step.
111: - Compute the initial second time derivative using backward differences.
112: - If using adaptivity, estimate the LTE of the initial step.
113: */
114: static PetscErrorCode TSAlpha_Restart(TS ts,PetscBool *initok)
115: {
116: TS_Alpha *th = (TS_Alpha*)ts->data;
117: PetscReal time_step;
118: PetscReal alpha_m,alpha_f,gamma,beta;
119: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
120: Vec V0 = ts->vec_dot, V1, V2 = th->V1;
121: PetscBool stageok;
125: VecDuplicate(X0,&X1);
126: VecDuplicate(V0,&V1);
128: /* Setup backward Euler with halved time step */
129: TSAlpha2GetParams(ts,&alpha_m,&alpha_f,&gamma,&beta);
130: TSAlpha2SetParams(ts,1,1,1,0.5);
131: TSGetTimeStep(ts,&time_step);
132: ts->time_step = time_step/2;
133: TSAlpha_StageTime(ts);
134: th->stage_time = ts->ptime;
135: VecZeroEntries(th->A0);
137: /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
138: th->stage_time += ts->time_step;
139: VecCopy(X0,th->X0);
140: VecCopy(V0,th->V0);
141: TSPreStage(ts,th->stage_time);
142: VecCopy(th->X0,X1);
143: TSAlpha_SNESSolve(ts,NULL,X1);
144: VecCopy(th->V1,V1);
145: TSPostStage(ts,th->stage_time,0,&X1);
146: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
147: if (!stageok) goto finally;
149: /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
150: th->stage_time += ts->time_step;
151: VecCopy(X1,th->X0);
152: VecCopy(V1,th->V0);
153: TSPreStage(ts,th->stage_time);
154: VecCopy(th->X0,X2);
155: TSAlpha_SNESSolve(ts,NULL,X2);
156: VecCopy(th->V1,V2);
157: TSPostStage(ts,th->stage_time,0,&X2);
158: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
159: if (!stageok) goto finally;
161: /* Compute A0 ~ dV/dt at t0 with backward differences */
162: VecZeroEntries(th->A0);
163: VecAXPY(th->A0,-3/ts->time_step,V0);
164: VecAXPY(th->A0,+4/ts->time_step,V1);
165: VecAXPY(th->A0,-1/ts->time_step,V2);
167: /* Rough, lower-order estimate LTE of the initial step */
168: if (th->vec_lte_work[0]) {
169: VecZeroEntries(th->vec_lte_work[0]);
170: VecAXPY(th->vec_lte_work[0],+2,X2);
171: VecAXPY(th->vec_lte_work[0],-4,X1);
172: VecAXPY(th->vec_lte_work[0],+2,X0);
173: }
174: if (th->vec_lte_work[1]) {
175: VecZeroEntries(th->vec_lte_work[1]);
176: VecAXPY(th->vec_lte_work[1],+2,V2);
177: VecAXPY(th->vec_lte_work[1],-4,V1);
178: VecAXPY(th->vec_lte_work[1],+2,V0);
179: }
181: finally:
182: /* Revert TSAlpha to the initial state (t0,X0,V0) */
183: if (initok) *initok = stageok;
184: TSSetTimeStep(ts,time_step);
185: TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
186: VecCopy(ts->vec_sol,th->X0);
187: VecCopy(ts->vec_dot,th->V0);
189: VecDestroy(&X1);
190: VecDestroy(&V1);
191: return(0);
192: }
194: static PetscErrorCode TSStep_Alpha(TS ts)
195: {
196: TS_Alpha *th = (TS_Alpha*)ts->data;
197: PetscInt rejections = 0;
198: PetscBool stageok,accept = PETSC_TRUE;
199: PetscReal next_time_step = ts->time_step;
203: PetscCitationsRegister(citation,&cited);
205: if (!ts->steprollback) {
206: if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
207: if (th->vec_dot_prev) { VecCopy(th->V0,th->vec_dot_prev); }
208: VecCopy(ts->vec_sol,th->X0);
209: VecCopy(ts->vec_dot,th->V0);
210: VecCopy(th->A1,th->A0);
211: }
213: th->status = TS_STEP_INCOMPLETE;
214: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
216: if (ts->steprestart) {
217: TSAlpha_Restart(ts,&stageok);
218: if (!stageok) goto reject_step;
219: }
221: TSAlpha_StageTime(ts);
222: VecCopy(th->X0,th->X1);
223: TSPreStage(ts,th->stage_time);
224: TSAlpha_SNESSolve(ts,NULL,th->X1);
225: TSPostStage(ts,th->stage_time,0,&th->Xa);
226: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
227: if (!stageok) goto reject_step;
229: th->status = TS_STEP_PENDING;
230: VecCopy(th->X1,ts->vec_sol);
231: VecCopy(th->V1,ts->vec_dot);
232: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
233: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
234: if (!accept) {
235: VecCopy(th->X0,ts->vec_sol);
236: VecCopy(th->V0,ts->vec_dot);
237: ts->time_step = next_time_step;
238: goto reject_step;
239: }
241: ts->ptime += ts->time_step;
242: ts->time_step = next_time_step;
243: break;
245: reject_step:
246: ts->reject++; accept = PETSC_FALSE;
247: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
248: ts->reason = TS_DIVERGED_STEP_REJECTED;
249: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
250: }
252: }
253: return(0);
254: }
256: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
257: {
258: TS_Alpha *th = (TS_Alpha*)ts->data;
259: Vec X = th->X1; /* X = solution */
260: Vec V = th->V1; /* V = solution */
261: Vec Y = th->vec_lte_work[0]; /* Y = X + LTE */
262: Vec Z = th->vec_lte_work[1]; /* Z = V + LTE */
263: PetscReal enormX,enormV,enormXa,enormVa,enormXr,enormVr;
267: if (!th->vec_sol_prev) {*wlte = -1; return(0);}
268: if (!th->vec_dot_prev) {*wlte = -1; return(0);}
269: if (!th->vec_lte_work[0]) {*wlte = -1; return(0);}
270: if (!th->vec_lte_work[1]) {*wlte = -1; return(0);}
271: if (ts->steprestart) {
272: /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
273: VecAXPY(Y,1,X);
274: VecAXPY(Z,1,V);
275: } else {
276: /* Compute LTE using backward differences with non-constant time step */
277: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
278: PetscReal a = 1 + h_prev/h;
279: PetscScalar scal[3]; Vec vecX[3],vecV[3];
280: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
281: vecX[0] = th->X1; vecX[1] = th->X0; vecX[2] = th->vec_sol_prev;
282: vecV[0] = th->V1; vecV[1] = th->V0; vecV[2] = th->vec_dot_prev;
283: VecCopy(X,Y);
284: VecMAXPY(Y,3,scal,vecX);
285: VecCopy(V,Z);
286: VecMAXPY(Z,3,scal,vecV);
287: }
288: /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
289: TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX,&enormXa,&enormXr);
290: TSErrorWeightedNorm(ts,V,Z,wnormtype,&enormV,&enormVa,&enormVr);
291: if (wnormtype == NORM_2)
292: *wlte = PetscSqrtReal(PetscSqr(enormX)/2 + PetscSqr(enormV)/2);
293: else
294: *wlte = PetscMax(enormX,enormV);
295: if (order) *order = 2;
296: return(0);
297: }
299: static PetscErrorCode TSRollBack_Alpha(TS ts)
300: {
301: TS_Alpha *th = (TS_Alpha*)ts->data;
305: VecCopy(th->X0,ts->vec_sol);
306: VecCopy(th->V0,ts->vec_dot);
307: return(0);
308: }
310: /*
311: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
312: {
313: TS_Alpha *th = (TS_Alpha*)ts->data;
314: PetscReal dt = t - ts->ptime;
318: VecCopy(ts->vec_dot,V);
319: VecAXPY(V,dt*(1-th->Gamma),th->A0);
320: VecAXPY(V,dt*th->Gamma,th->A1);
321: VecCopy(ts->vec_sol,X);
322: VecAXPY(X,dt,V);
323: VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
324: VecAXPY(X,dt*dt*th->Beta,th->A1);
325: return(0);
326: }
327: */
329: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
330: {
331: TS_Alpha *th = (TS_Alpha*)ts->data;
332: PetscReal ta = th->stage_time;
333: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
337: TSAlpha_StageVecs(ts,X);
338: /* F = Function(ta,Xa,Va,Aa) */
339: TSComputeI2Function(ts,ta,Xa,Va,Aa,F);
340: VecScale(F,th->scale_F);
341: return(0);
342: }
344: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,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;
349: PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
353: /* J,P = Jacobian(ta,Xa,Va,Aa) */
354: TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P);
355: return(0);
356: }
358: static PetscErrorCode TSReset_Alpha(TS ts)
359: {
360: TS_Alpha *th = (TS_Alpha*)ts->data;
364: VecDestroy(&th->X0);
365: VecDestroy(&th->Xa);
366: VecDestroy(&th->X1);
367: VecDestroy(&th->V0);
368: VecDestroy(&th->Va);
369: VecDestroy(&th->V1);
370: VecDestroy(&th->A0);
371: VecDestroy(&th->Aa);
372: VecDestroy(&th->A1);
373: VecDestroy(&th->vec_sol_prev);
374: VecDestroy(&th->vec_dot_prev);
375: VecDestroy(&th->vec_lte_work[0]);
376: VecDestroy(&th->vec_lte_work[1]);
377: return(0);
378: }
380: static PetscErrorCode TSDestroy_Alpha(TS ts)
381: {
385: TSReset_Alpha(ts);
386: PetscFree(ts->data);
388: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL);
389: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL);
390: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL);
391: return(0);
392: }
394: static PetscErrorCode TSSetUp_Alpha(TS ts)
395: {
396: TS_Alpha *th = (TS_Alpha*)ts->data;
397: PetscBool match;
401: VecDuplicate(ts->vec_sol,&th->X0);
402: VecDuplicate(ts->vec_sol,&th->Xa);
403: VecDuplicate(ts->vec_sol,&th->X1);
404: VecDuplicate(ts->vec_sol,&th->V0);
405: VecDuplicate(ts->vec_sol,&th->Va);
406: VecDuplicate(ts->vec_sol,&th->V1);
407: VecDuplicate(ts->vec_sol,&th->A0);
408: VecDuplicate(ts->vec_sol,&th->Aa);
409: VecDuplicate(ts->vec_sol,&th->A1);
411: TSGetAdapt(ts,&ts->adapt);
412: TSAdaptCandidatesClear(ts->adapt);
413: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
414: if (!match) {
415: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
416: VecDuplicate(ts->vec_sol,&th->vec_dot_prev);
417: VecDuplicate(ts->vec_sol,&th->vec_lte_work[0]);
418: VecDuplicate(ts->vec_sol,&th->vec_lte_work[1]);
419: }
421: TSGetSNES(ts,&ts->snes);
422: return(0);
423: }
425: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
426: {
427: TS_Alpha *th = (TS_Alpha*)ts->data;
431: PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
432: {
433: PetscBool flg;
434: PetscReal radius = 1;
435: PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg);
436: if (flg) {TSAlpha2SetRadius(ts,radius);}
437: PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL);
438: PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL);
439: PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL);
440: PetscOptionsReal("-ts_alpha_beta","Algorithmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL);
441: TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta);
442: }
443: PetscOptionsTail();
444: return(0);
445: }
447: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
448: {
449: TS_Alpha *th = (TS_Alpha*)ts->data;
450: PetscBool iascii;
454: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
455: if (iascii) {
456: 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);
457: }
458: return(0);
459: }
461: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)
462: {
463: PetscReal alpha_m,alpha_f,gamma,beta;
467: if (radius < 0 || radius > 1) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
468: alpha_m = (2-radius)/(1+radius);
469: alpha_f = 1/(1+radius);
470: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
471: beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta;
472: TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta);
473: return(0);
474: }
476: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
477: {
478: TS_Alpha *th = (TS_Alpha*)ts->data;
479: PetscReal tol = 100*PETSC_MACHINE_EPSILON;
480: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
483: th->Alpha_m = alpha_m;
484: th->Alpha_f = alpha_f;
485: th->Gamma = gamma;
486: th->Beta = beta;
487: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
488: return(0);
489: }
491: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
492: {
493: TS_Alpha *th = (TS_Alpha*)ts->data;
496: if (alpha_m) *alpha_m = th->Alpha_m;
497: if (alpha_f) *alpha_f = th->Alpha_f;
498: if (gamma) *gamma = th->Gamma;
499: if (beta) *beta = th->Beta;
500: return(0);
501: }
503: /*MC
504: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
505: for second-order systems
507: Level: beginner
509: References:
510: J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
511: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
512: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
514: .seealso: TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams()
515: M*/
516: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
517: {
518: TS_Alpha *th;
522: ts->ops->reset = TSReset_Alpha;
523: ts->ops->destroy = TSDestroy_Alpha;
524: ts->ops->view = TSView_Alpha;
525: ts->ops->setup = TSSetUp_Alpha;
526: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
527: ts->ops->step = TSStep_Alpha;
528: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
529: ts->ops->rollback = TSRollBack_Alpha;
530: /*ts->ops->interpolate = TSInterpolate_Alpha;*/
531: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
532: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
533: ts->default_adapt_type = TSADAPTNONE;
535: ts->usessnes = PETSC_TRUE;
537: PetscNewLog(ts,&th);
538: ts->data = (void*)th;
540: th->Alpha_m = 0.5;
541: th->Alpha_f = 0.5;
542: th->Gamma = 0.5;
543: th->Beta = 0.25;
544: th->order = 2;
546: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha);
547: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha);
548: PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha);
549: return(0);
550: }
552: /*@
553: TSAlpha2SetRadius - sets the desired spectral radius of the method
554: (i.e. high-frequency numerical damping)
556: Logically Collective on TS
558: The algorithmic parameters \alpha_m and \alpha_f of the
559: generalized-\alpha method can be computed in terms of a specified
560: spectral radius \rho in [0,1] for infinite time step in order to
561: control high-frequency numerical damping:
562: \alpha_m = (2-\rho)/(1+\rho)
563: \alpha_f = 1/(1+\rho)
565: Input Parameter:
566: + ts - timestepping context
567: - radius - the desired spectral radius
569: Options Database:
570: . -ts_alpha_radius <radius>
572: Level: intermediate
574: .seealso: TSAlpha2SetParams(), TSAlpha2GetParams()
575: @*/
576: PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius)
577: {
583: if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
584: PetscTryMethod(ts,"TSAlpha2SetRadius_C",(TS,PetscReal),(ts,radius));
585: return(0);
586: }
588: /*@
589: TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2
591: Logically Collective on TS
593: Second-order accuracy can be obtained so long as:
594: \gamma = 1/2 + alpha_m - alpha_f
595: \beta = 1/4 (1 + alpha_m - alpha_f)^2
597: Unconditional stability requires:
598: \alpha_m >= \alpha_f >= 1/2
601: Input Parameter:
602: + ts - timestepping context
603: . \alpha_m - algorithmic paramenter
604: . \alpha_f - algorithmic paramenter
605: . \gamma - algorithmic paramenter
606: - \beta - algorithmic paramenter
608: Options Database:
609: + -ts_alpha_alpha_m <alpha_m>
610: . -ts_alpha_alpha_f <alpha_f>
611: . -ts_alpha_gamma <gamma>
612: - -ts_alpha_beta <beta>
614: Note:
615: Use of this function is normally only required to hack TSALPHA2 to
616: use a modified integration scheme. Users should call
617: TSAlpha2SetRadius() to set the desired spectral radius of the methods
618: (i.e. high-frequency damping) in order so select optimal values for
619: these parameters.
621: Level: advanced
623: .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams()
624: @*/
625: PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)
626: {
635: PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta));
636: return(0);
637: }
639: /*@
640: TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2
642: Not Collective
644: Input Parameter:
645: . ts - timestepping context
647: Output Parameters:
648: + \alpha_m - algorithmic parameter
649: . \alpha_f - algorithmic parameter
650: . \gamma - algorithmic parameter
651: - \beta - algorithmic parameter
653: Note:
654: Use of this function is normally only required to hack TSALPHA2 to
655: use a modified integration scheme. Users should call
656: TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral
657: radius of the method) in order so select optimal values for these
658: parameters.
660: Level: advanced
662: .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams()
663: @*/
664: PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta)
665: {
674: PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta));
675: return(0);
676: }