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