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 CHKERRQ(VecWAXPY(A1,-1.0,X0,X1)); 77 CHKERRQ(VecAXPY (A1,-dt,V0)); 78 CHKERRQ(VecAXPBY(A1,-(1-2*Beta)/(2*Beta),1/(dt*dt*Beta),A0)); 79 /* V1 = ... */ 80 CHKERRQ(VecWAXPY(V1,(1.0-Gamma)/Gamma,A0,A1)); 81 CHKERRQ(VecAYPX (V1,dt*Gamma,V0)); 82 /* Xa = X0 + Alpha_f*(X1-X0) */ 83 CHKERRQ(VecWAXPY(Xa,-1.0,X0,X1)); 84 CHKERRQ(VecAYPX (Xa,Alpha_f,X0)); 85 /* Va = V0 + Alpha_f*(V1-V0) */ 86 CHKERRQ(VecWAXPY(Va,-1.0,V0,V1)); 87 CHKERRQ(VecAYPX (Va,Alpha_f,V0)); 88 /* Aa = A0 + Alpha_m*(A1-A0) */ 89 CHKERRQ(VecWAXPY(Aa,-1.0,A0,A1)); 90 CHKERRQ(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 CHKERRQ(SNESSolve(ts->snes,b,x)); 100 CHKERRQ(SNESGetIterationNumber(ts->snes,&nits)); 101 CHKERRQ(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 CHKERRQ(VecDuplicate(X0,&X1)); 123 CHKERRQ(VecDuplicate(V0,&V1)); 124 125 /* Setup backward Euler with halved time step */ 126 CHKERRQ(TSAlpha2GetParams(ts,&alpha_m,&alpha_f,&gamma,&beta)); 127 CHKERRQ(TSAlpha2SetParams(ts,1,1,1,0.5)); 128 CHKERRQ(TSGetTimeStep(ts,&time_step)); 129 ts->time_step = time_step/2; 130 CHKERRQ(TSAlpha_StageTime(ts)); 131 th->stage_time = ts->ptime; 132 CHKERRQ(VecZeroEntries(th->A0)); 133 134 /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */ 135 th->stage_time += ts->time_step; 136 CHKERRQ(VecCopy(X0,th->X0)); 137 CHKERRQ(VecCopy(V0,th->V0)); 138 CHKERRQ(TSPreStage(ts,th->stage_time)); 139 CHKERRQ(VecCopy(th->X0,X1)); 140 CHKERRQ(TSAlpha_SNESSolve(ts,NULL,X1)); 141 CHKERRQ(VecCopy(th->V1,V1)); 142 CHKERRQ(TSPostStage(ts,th->stage_time,0,&X1)); 143 CHKERRQ(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 CHKERRQ(VecCopy(X1,th->X0)); 149 CHKERRQ(VecCopy(V1,th->V0)); 150 CHKERRQ(TSPreStage(ts,th->stage_time)); 151 CHKERRQ(VecCopy(th->X0,X2)); 152 CHKERRQ(TSAlpha_SNESSolve(ts,NULL,X2)); 153 CHKERRQ(VecCopy(th->V1,V2)); 154 CHKERRQ(TSPostStage(ts,th->stage_time,0,&X2)); 155 CHKERRQ(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 CHKERRQ(VecZeroEntries(th->A0)); 160 CHKERRQ(VecAXPY(th->A0,-3/ts->time_step,V0)); 161 CHKERRQ(VecAXPY(th->A0,+4/ts->time_step,V1)); 162 CHKERRQ(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 CHKERRQ(VecZeroEntries(th->vec_lte_work[0])); 167 CHKERRQ(VecAXPY(th->vec_lte_work[0],+2,X2)); 168 CHKERRQ(VecAXPY(th->vec_lte_work[0],-4,X1)); 169 CHKERRQ(VecAXPY(th->vec_lte_work[0],+2,X0)); 170 } 171 if (th->vec_lte_work[1]) { 172 CHKERRQ(VecZeroEntries(th->vec_lte_work[1])); 173 CHKERRQ(VecAXPY(th->vec_lte_work[1],+2,V2)); 174 CHKERRQ(VecAXPY(th->vec_lte_work[1],-4,V1)); 175 CHKERRQ(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 CHKERRQ(TSSetTimeStep(ts,time_step)); 182 CHKERRQ(TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta)); 183 CHKERRQ(VecCopy(ts->vec_sol,th->X0)); 184 CHKERRQ(VecCopy(ts->vec_dot,th->V0)); 185 186 CHKERRQ(VecDestroy(&X1)); 187 CHKERRQ(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 CHKERRQ(PetscCitationsRegister(citation,&cited)); 200 201 if (!ts->steprollback) { 202 if (th->vec_sol_prev) CHKERRQ(VecCopy(th->X0,th->vec_sol_prev)); 203 if (th->vec_dot_prev) CHKERRQ(VecCopy(th->V0,th->vec_dot_prev)); 204 CHKERRQ(VecCopy(ts->vec_sol,th->X0)); 205 CHKERRQ(VecCopy(ts->vec_dot,th->V0)); 206 CHKERRQ(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 CHKERRQ(TSAlpha_Restart(ts,&stageok)); 214 if (!stageok) goto reject_step; 215 } 216 217 CHKERRQ(TSAlpha_StageTime(ts)); 218 CHKERRQ(VecCopy(th->X0,th->X1)); 219 CHKERRQ(TSPreStage(ts,th->stage_time)); 220 CHKERRQ(TSAlpha_SNESSolve(ts,NULL,th->X1)); 221 CHKERRQ(TSPostStage(ts,th->stage_time,0,&th->Xa)); 222 CHKERRQ(TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok)); 223 if (!stageok) goto reject_step; 224 225 th->status = TS_STEP_PENDING; 226 CHKERRQ(VecCopy(th->X1,ts->vec_sol)); 227 CHKERRQ(VecCopy(th->V1,ts->vec_dot)); 228 CHKERRQ(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 CHKERRQ(VecCopy(th->X0,ts->vec_sol)); 232 CHKERRQ(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 CHKERRQ(PetscInfo(ts,"Step=%D, step rejections %D 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 CHKERRQ(VecAXPY(Y,1,X)); 269 CHKERRQ(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 CHKERRQ(VecCopy(X,Y)); 279 CHKERRQ(VecMAXPY(Y,3,scal,vecX)); 280 CHKERRQ(VecCopy(V,Z)); 281 CHKERRQ(VecMAXPY(Z,3,scal,vecV)); 282 } 283 /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */ 284 CHKERRQ(TSErrorWeightedNorm(ts,X,Y,wnormtype,&enormX,&enormXa,&enormXr)); 285 CHKERRQ(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 CHKERRQ(VecCopy(th->X0,ts->vec_sol)); 300 CHKERRQ(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 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 CHKERRQ(VecCopy(ts->vec_dot,V)); 313 CHKERRQ(VecAXPY(V,dt*(1-th->Gamma),th->A0)); 314 CHKERRQ(VecAXPY(V,dt*th->Gamma,th->A1)); 315 CHKERRQ(VecCopy(ts->vec_sol,X)); 316 CHKERRQ(VecAXPY(X,dt,V)); 317 CHKERRQ(VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0)); 318 CHKERRQ(VecAXPY(X,dt*dt*th->Beta,th->A1)); 319 PetscFunctionReturn(0); 320 } 321 */ 322 323 static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts) 324 { 325 TS_Alpha *th = (TS_Alpha*)ts->data; 326 PetscReal ta = th->stage_time; 327 Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa; 328 329 PetscFunctionBegin; 330 CHKERRQ(TSAlpha_StageVecs(ts,X)); 331 /* F = Function(ta,Xa,Va,Aa) */ 332 CHKERRQ(TSComputeI2Function(ts,ta,Xa,Va,Aa,F)); 333 CHKERRQ(VecScale(F,th->scale_F)); 334 PetscFunctionReturn(0); 335 } 336 337 static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts) 338 { 339 TS_Alpha *th = (TS_Alpha*)ts->data; 340 PetscReal ta = th->stage_time; 341 Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa; 342 PetscReal dVdX = th->shift_V, dAdX = th->shift_A; 343 344 PetscFunctionBegin; 345 /* J,P = Jacobian(ta,Xa,Va,Aa) */ 346 CHKERRQ(TSComputeI2Jacobian(ts,ta,Xa,Va,Aa,dVdX,dAdX,J,P)); 347 PetscFunctionReturn(0); 348 } 349 350 static PetscErrorCode TSReset_Alpha(TS ts) 351 { 352 TS_Alpha *th = (TS_Alpha*)ts->data; 353 354 PetscFunctionBegin; 355 CHKERRQ(VecDestroy(&th->X0)); 356 CHKERRQ(VecDestroy(&th->Xa)); 357 CHKERRQ(VecDestroy(&th->X1)); 358 CHKERRQ(VecDestroy(&th->V0)); 359 CHKERRQ(VecDestroy(&th->Va)); 360 CHKERRQ(VecDestroy(&th->V1)); 361 CHKERRQ(VecDestroy(&th->A0)); 362 CHKERRQ(VecDestroy(&th->Aa)); 363 CHKERRQ(VecDestroy(&th->A1)); 364 CHKERRQ(VecDestroy(&th->vec_sol_prev)); 365 CHKERRQ(VecDestroy(&th->vec_dot_prev)); 366 CHKERRQ(VecDestroy(&th->vec_lte_work[0])); 367 CHKERRQ(VecDestroy(&th->vec_lte_work[1])); 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode TSDestroy_Alpha(TS ts) 372 { 373 PetscFunctionBegin; 374 CHKERRQ(TSReset_Alpha(ts)); 375 CHKERRQ(PetscFree(ts->data)); 376 377 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",NULL)); 378 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",NULL)); 379 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",NULL)); 380 PetscFunctionReturn(0); 381 } 382 383 static PetscErrorCode TSSetUp_Alpha(TS ts) 384 { 385 TS_Alpha *th = (TS_Alpha*)ts->data; 386 PetscBool match; 387 388 PetscFunctionBegin; 389 CHKERRQ(VecDuplicate(ts->vec_sol,&th->X0)); 390 CHKERRQ(VecDuplicate(ts->vec_sol,&th->Xa)); 391 CHKERRQ(VecDuplicate(ts->vec_sol,&th->X1)); 392 CHKERRQ(VecDuplicate(ts->vec_sol,&th->V0)); 393 CHKERRQ(VecDuplicate(ts->vec_sol,&th->Va)); 394 CHKERRQ(VecDuplicate(ts->vec_sol,&th->V1)); 395 CHKERRQ(VecDuplicate(ts->vec_sol,&th->A0)); 396 CHKERRQ(VecDuplicate(ts->vec_sol,&th->Aa)); 397 CHKERRQ(VecDuplicate(ts->vec_sol,&th->A1)); 398 399 CHKERRQ(TSGetAdapt(ts,&ts->adapt)); 400 CHKERRQ(TSAdaptCandidatesClear(ts->adapt)); 401 CHKERRQ(PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match)); 402 if (!match) { 403 CHKERRQ(VecDuplicate(ts->vec_sol,&th->vec_sol_prev)); 404 CHKERRQ(VecDuplicate(ts->vec_sol,&th->vec_dot_prev)); 405 CHKERRQ(VecDuplicate(ts->vec_sol,&th->vec_lte_work[0])); 406 CHKERRQ(VecDuplicate(ts->vec_sol,&th->vec_lte_work[1])); 407 } 408 409 CHKERRQ(TSGetSNES(ts,&ts->snes)); 410 PetscFunctionReturn(0); 411 } 412 413 static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts) 414 { 415 TS_Alpha *th = (TS_Alpha*)ts->data; 416 417 PetscFunctionBegin; 418 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options")); 419 { 420 PetscBool flg; 421 PetscReal radius = 1; 422 CHKERRQ(PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlpha2SetRadius",radius,&radius,&flg)); 423 if (flg) CHKERRQ(TSAlpha2SetRadius(ts,radius)); 424 CHKERRQ(PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlpha2SetParams",th->Alpha_m,&th->Alpha_m,NULL)); 425 CHKERRQ(PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlpha2SetParams",th->Alpha_f,&th->Alpha_f,NULL)); 426 CHKERRQ(PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlpha2SetParams",th->Gamma,&th->Gamma,NULL)); 427 CHKERRQ(PetscOptionsReal("-ts_alpha_beta","Algorithmic parameter beta","TSAlpha2SetParams",th->Beta,&th->Beta,NULL)); 428 CHKERRQ(TSAlpha2SetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma,th->Beta)); 429 } 430 CHKERRQ(PetscOptionsTail()); 431 PetscFunctionReturn(0); 432 } 433 434 static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer) 435 { 436 TS_Alpha *th = (TS_Alpha*)ts->data; 437 PetscBool iascii; 438 439 PetscFunctionBegin; 440 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 441 if (iascii) { 442 CHKERRQ(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)); 443 } 444 PetscFunctionReturn(0); 445 } 446 447 static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius) 448 { 449 PetscReal alpha_m,alpha_f,gamma,beta; 450 451 PetscFunctionBegin; 452 PetscCheckFalse(radius < 0 || radius > 1,PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius); 453 alpha_m = (2-radius)/(1+radius); 454 alpha_f = 1/(1+radius); 455 gamma = (PetscReal)0.5 + alpha_m - alpha_f; 456 beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f); beta *= beta; 457 CHKERRQ(TSAlpha2SetParams(ts,alpha_m,alpha_f,gamma,beta)); 458 PetscFunctionReturn(0); 459 } 460 461 static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta) 462 { 463 TS_Alpha *th = (TS_Alpha*)ts->data; 464 PetscReal tol = 100*PETSC_MACHINE_EPSILON; 465 PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma; 466 467 PetscFunctionBegin; 468 th->Alpha_m = alpha_m; 469 th->Alpha_f = alpha_f; 470 th->Gamma = gamma; 471 th->Beta = beta; 472 th->order = (PetscAbsReal(res) < tol) ? 2 : 1; 473 PetscFunctionReturn(0); 474 } 475 476 static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta) 477 { 478 TS_Alpha *th = (TS_Alpha*)ts->data; 479 480 PetscFunctionBegin; 481 if (alpha_m) *alpha_m = th->Alpha_m; 482 if (alpha_f) *alpha_f = th->Alpha_f; 483 if (gamma) *gamma = th->Gamma; 484 if (beta) *beta = th->Beta; 485 PetscFunctionReturn(0); 486 } 487 488 /*MC 489 TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method 490 for second-order systems 491 492 Level: beginner 493 494 References: 495 . * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural 496 Dynamics with Improved Numerical Dissipation: The Generalized-alpha 497 Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993. 498 499 .seealso: TS, TSCreate(), TSSetType(), TSAlpha2SetRadius(), TSAlpha2SetParams() 500 M*/ 501 PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts) 502 { 503 TS_Alpha *th; 504 505 PetscFunctionBegin; 506 ts->ops->reset = TSReset_Alpha; 507 ts->ops->destroy = TSDestroy_Alpha; 508 ts->ops->view = TSView_Alpha; 509 ts->ops->setup = TSSetUp_Alpha; 510 ts->ops->setfromoptions = TSSetFromOptions_Alpha; 511 ts->ops->step = TSStep_Alpha; 512 ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha; 513 ts->ops->rollback = TSRollBack_Alpha; 514 /*ts->ops->interpolate = TSInterpolate_Alpha;*/ 515 ts->ops->snesfunction = SNESTSFormFunction_Alpha; 516 ts->ops->snesjacobian = SNESTSFormJacobian_Alpha; 517 ts->default_adapt_type = TSADAPTNONE; 518 519 ts->usessnes = PETSC_TRUE; 520 521 CHKERRQ(PetscNewLog(ts,&th)); 522 ts->data = (void*)th; 523 524 th->Alpha_m = 0.5; 525 th->Alpha_f = 0.5; 526 th->Gamma = 0.5; 527 th->Beta = 0.25; 528 th->order = 2; 529 530 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetRadius_C",TSAlpha2SetRadius_Alpha)); 531 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2SetParams_C",TSAlpha2SetParams_Alpha)); 532 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSAlpha2GetParams_C",TSAlpha2GetParams_Alpha)); 533 PetscFunctionReturn(0); 534 } 535 536 /*@ 537 TSAlpha2SetRadius - sets the desired spectral radius of the method 538 (i.e. high-frequency numerical damping) 539 540 Logically Collective on TS 541 542 The algorithmic parameters \alpha_m and \alpha_f of the 543 generalized-\alpha method can be computed in terms of a specified 544 spectral radius \rho in [0,1] for infinite time step in order to 545 control high-frequency numerical damping: 546 \alpha_m = (2-\rho)/(1+\rho) 547 \alpha_f = 1/(1+\rho) 548 549 Input Parameters: 550 + ts - timestepping context 551 - radius - the desired spectral radius 552 553 Options Database: 554 . -ts_alpha_radius <radius> - set the desired spectral radius 555 556 Level: intermediate 557 558 .seealso: TSAlpha2SetParams(), TSAlpha2GetParams() 559 @*/ 560 PetscErrorCode TSAlpha2SetRadius(TS ts,PetscReal radius) 561 { 562 PetscFunctionBegin; 563 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 564 PetscValidLogicalCollectiveReal(ts,radius,2); 565 PetscCheckFalse(radius < 0 || radius > 1,((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius); 566 CHKERRQ(PetscTryMethod(ts,"TSAlpha2SetRadius_C",(TS,PetscReal),(ts,radius))); 567 PetscFunctionReturn(0); 568 } 569 570 /*@ 571 TSAlpha2SetParams - sets the algorithmic parameters for TSALPHA2 572 573 Logically Collective on TS 574 575 Second-order accuracy can be obtained so long as: 576 \gamma = 1/2 + alpha_m - alpha_f 577 \beta = 1/4 (1 + alpha_m - alpha_f)^2 578 579 Unconditional stability requires: 580 \alpha_m >= \alpha_f >= 1/2 581 582 Input Parameters: 583 + ts - timestepping context 584 . \alpha_m - algorithmic parameter 585 . \alpha_f - algorithmic parameter 586 . \gamma - algorithmic parameter 587 - \beta - algorithmic parameter 588 589 Options Database: 590 + -ts_alpha_alpha_m <alpha_m> - set alpha_m 591 . -ts_alpha_alpha_f <alpha_f> - set alpha_f 592 . -ts_alpha_gamma <gamma> - set gamma 593 - -ts_alpha_beta <beta> - set beta 594 595 Note: 596 Use of this function is normally only required to hack TSALPHA2 to 597 use a modified integration scheme. Users should call 598 TSAlpha2SetRadius() to set the desired spectral radius of the methods 599 (i.e. high-frequency damping) in order so select optimal values for 600 these parameters. 601 602 Level: advanced 603 604 .seealso: TSAlpha2SetRadius(), TSAlpha2GetParams() 605 @*/ 606 PetscErrorCode TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta) 607 { 608 PetscFunctionBegin; 609 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 610 PetscValidLogicalCollectiveReal(ts,alpha_m,2); 611 PetscValidLogicalCollectiveReal(ts,alpha_f,3); 612 PetscValidLogicalCollectiveReal(ts,gamma,4); 613 PetscValidLogicalCollectiveReal(ts,beta,5); 614 CHKERRQ(PetscTryMethod(ts,"TSAlpha2SetParams_C",(TS,PetscReal,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma,beta))); 615 PetscFunctionReturn(0); 616 } 617 618 /*@ 619 TSAlpha2GetParams - gets the algorithmic parameters for TSALPHA2 620 621 Not Collective 622 623 Input Parameter: 624 . ts - timestepping context 625 626 Output Parameters: 627 + \alpha_m - algorithmic parameter 628 . \alpha_f - algorithmic parameter 629 . \gamma - algorithmic parameter 630 - \beta - algorithmic parameter 631 632 Note: 633 Use of this function is normally only required to hack TSALPHA2 to 634 use a modified integration scheme. Users should call 635 TSAlpha2SetRadius() to set the high-frequency damping (i.e. spectral 636 radius of the method) in order so select optimal values for these 637 parameters. 638 639 Level: advanced 640 641 .seealso: TSAlpha2SetRadius(), TSAlpha2SetParams() 642 @*/ 643 PetscErrorCode TSAlpha2GetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma,PetscReal *beta) 644 { 645 PetscFunctionBegin; 646 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 647 if (alpha_m) PetscValidRealPointer(alpha_m,2); 648 if (alpha_f) PetscValidRealPointer(alpha_f,3); 649 if (gamma) PetscValidRealPointer(gamma,4); 650 if (beta) PetscValidRealPointer(beta,5); 651 CHKERRQ(PetscUseMethod(ts,"TSAlpha2GetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma,beta))); 652 PetscFunctionReturn(0); 653 } 654