1 /* 2 Code for timestepping with implicit generalized-\alpha method 3 for first 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{Jansen2000,\n" 10 " title = {A generalized-$\\alpha$ method for integrating the filtered {N}avier--{S}tokes equations with a stabilized finite element method},\n" 11 " author = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},\n" 12 " journal = {Computer Methods in Applied Mechanics and Engineering},\n" 13 " volume = {190},\n" 14 " number = {3--4},\n" 15 " pages = {305--319},\n" 16 " year = {2000},\n" 17 " issn = {0045-7825},\n" 18 " doi = {http://dx.doi.org/10.1016/S0045-7825(00)00203-6}\n}\n"; 19 20 typedef struct { 21 PetscReal stage_time; 22 PetscReal shift_V; 23 PetscReal scale_F; 24 Vec X0,Xa,X1; 25 Vec V0,Va,V1; 26 27 PetscReal Alpha_m; 28 PetscReal Alpha_f; 29 PetscReal Gamma; 30 PetscInt order; 31 32 PetscBool adapt; 33 Vec vec_sol_prev; 34 Vec vec_lte_work; 35 36 TSStepStatus status; 37 } TS_Alpha; 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "TSAlpha_StageTime" 41 static PetscErrorCode TSAlpha_StageTime(TS ts) 42 { 43 TS_Alpha *th = (TS_Alpha*)ts->data; 44 PetscReal t = ts->ptime; 45 PetscReal dt = ts->time_step; 46 PetscReal Alpha_m = th->Alpha_m; 47 PetscReal Alpha_f = th->Alpha_f; 48 PetscReal Gamma = th->Gamma; 49 50 PetscFunctionBegin; 51 th->stage_time = t + Alpha_f*dt; 52 th->shift_V = Alpha_m/(Alpha_f*Gamma*dt); 53 th->scale_F = 1/Alpha_f; 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "TSAlpha_StageVecs" 59 static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X) 60 { 61 TS_Alpha *th = (TS_Alpha*)ts->data; 62 Vec X1 = X, V1 = th->V1; 63 Vec Xa = th->Xa, Va = th->Va; 64 Vec X0 = th->X0, V0 = th->V0; 65 PetscReal dt = ts->time_step; 66 PetscReal Alpha_m = th->Alpha_m; 67 PetscReal Alpha_f = th->Alpha_f; 68 PetscReal Gamma = th->Gamma; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */ 73 ierr = VecWAXPY(V1,-1.0,X0,X1);CHKERRQ(ierr); 74 ierr = VecAXPBY(V1,1-1/Gamma,1/(Gamma*dt),V0);CHKERRQ(ierr); 75 /* Xa = X0 + Alpha_f*(X1-X0) */ 76 ierr = VecWAXPY(Xa,-1.0,X0,X1);CHKERRQ(ierr); 77 ierr = VecAYPX(Xa,Alpha_f,X0);CHKERRQ(ierr); 78 /* Va = V0 + Alpha_m*(V1-V0) */ 79 ierr = VecWAXPY(Va,-1.0,V0,V1);CHKERRQ(ierr); 80 ierr = VecAYPX(Va,Alpha_m,V0);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "TS_SNESSolve" 86 static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x) 87 { 88 PetscInt nits,lits; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = SNESSolve(ts->snes,b,x);CHKERRQ(ierr); 93 ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr); 94 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 95 ts->snes_its += nits; ts->ksp_its += lits; 96 PetscFunctionReturn(0); 97 } 98 99 /* 100 Compute a consistent initial state for the generalized-alpha method. 101 - Solve two successive backward Euler steps with halved time step. 102 - Compute the initial time derivative using backward differences. 103 - If using adaptivity, estimate the LTE of the initial step. 104 */ 105 #undef __FUNCT__ 106 #define __FUNCT__ "TSAlpha_Restart" 107 static PetscErrorCode TSAlpha_Restart(TS ts,PetscBool *initok) 108 { 109 TS_Alpha *th = (TS_Alpha*)ts->data; 110 PetscReal time_step; 111 PetscReal alpha_m,alpha_f,gamma; 112 Vec X0 = ts->vec_sol, X1, X2 = th->X1; 113 PetscBool stageok; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = VecDuplicate(X0,&X1);CHKERRQ(ierr); 118 119 /* Setup backward Euler with halved time step */ 120 ierr = TSAlphaGetParams(ts,&alpha_m,&alpha_f,&gamma);CHKERRQ(ierr); 121 ierr = TSAlphaSetParams(ts,1,1,1);CHKERRQ(ierr); 122 ierr = TSGetTimeStep(ts,&time_step);CHKERRQ(ierr); 123 ts->time_step = time_step/2; 124 ierr = TSAlpha_StageTime(ts);CHKERRQ(ierr); 125 th->stage_time = ts->ptime; 126 ierr = VecZeroEntries(th->V0);CHKERRQ(ierr); 127 128 /* First BE step, (t0,X0) -> (t1,X1) */ 129 th->stage_time += ts->time_step; 130 ierr = VecCopy(X0,th->X0);CHKERRQ(ierr); 131 ierr = TSPreStage(ts,th->stage_time);CHKERRQ(ierr); 132 ierr = VecCopy(th->X0,X1);CHKERRQ(ierr); 133 ierr = TS_SNESSolve(ts,NULL,X1);CHKERRQ(ierr); 134 ierr = TSPostStage(ts,th->stage_time,0,&X1);CHKERRQ(ierr); 135 ierr = TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);CHKERRQ(ierr); 136 if (!stageok) goto finally; 137 138 /* Second BE step, (t1,X1) -> (t2,X2) */ 139 th->stage_time += ts->time_step; 140 ierr = VecCopy(X1,th->X0);CHKERRQ(ierr); 141 ierr = TSPreStage(ts,th->stage_time);CHKERRQ(ierr); 142 ierr = VecCopy(th->X0,X2);CHKERRQ(ierr); 143 ierr = TS_SNESSolve(ts,NULL,X2);CHKERRQ(ierr); 144 ierr = TSPostStage(ts,th->stage_time,0,&X2);CHKERRQ(ierr); 145 ierr = TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X2,&stageok);CHKERRQ(ierr); 146 if (!stageok) goto finally; 147 148 /* Compute V0 ~ dX/dt at t0 with backward differences */ 149 ierr = VecZeroEntries(th->V0);CHKERRQ(ierr); 150 ierr = VecAXPY(th->V0,-3/ts->time_step,X0);CHKERRQ(ierr); 151 ierr = VecAXPY(th->V0,+4/ts->time_step,X1);CHKERRQ(ierr); 152 ierr = VecAXPY(th->V0,-1/ts->time_step,X2);CHKERRQ(ierr); 153 154 /* Rough, lower-order estimate LTE of the initial step */ 155 if (th->adapt) { 156 ierr = VecZeroEntries(th->vec_lte_work);CHKERRQ(ierr); 157 ierr = VecAXPY(th->vec_lte_work,+2,X2);CHKERRQ(ierr); 158 ierr = VecAXPY(th->vec_lte_work,-4,X1);CHKERRQ(ierr); 159 ierr = VecAXPY(th->vec_lte_work,+2,X0);CHKERRQ(ierr); 160 } 161 162 finally: 163 /* Revert TSAlpha to the initial state (t0,X0) */ 164 if (initok) *initok = stageok; 165 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 166 ierr = TSAlphaSetParams(ts,alpha_m,alpha_f,gamma);CHKERRQ(ierr); 167 ierr = VecCopy(ts->vec_sol,th->X0);CHKERRQ(ierr); 168 169 ierr = VecDestroy(&X1);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "TSStep_Alpha" 175 static PetscErrorCode TSStep_Alpha(TS ts) 176 { 177 TS_Alpha *th = (TS_Alpha*)ts->data; 178 PetscInt rejections = 0; 179 PetscBool stageok,accept = PETSC_TRUE; 180 PetscReal next_time_step = ts->time_step; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 185 186 if (!ts->steprollback) { 187 if (th->adapt) { ierr = VecCopy(th->X0,th->vec_sol_prev);CHKERRQ(ierr); } 188 ierr = VecCopy(ts->vec_sol,th->X0);CHKERRQ(ierr); 189 ierr = VecCopy(th->V1,th->V0);CHKERRQ(ierr); 190 } 191 192 th->status = TS_STEP_INCOMPLETE; 193 while (!ts->reason && th->status != TS_STEP_COMPLETE) { 194 195 if (ts->steprestart) { 196 ierr = TSAlpha_Restart(ts,&stageok);CHKERRQ(ierr); 197 if (!stageok) goto reject_step; 198 } 199 200 ierr = TSAlpha_StageTime(ts);CHKERRQ(ierr); 201 ierr = VecCopy(th->X0,th->X1);CHKERRQ(ierr); 202 ierr = TSPreStage(ts,th->stage_time);CHKERRQ(ierr); 203 ierr = TS_SNESSolve(ts,NULL,th->X1);CHKERRQ(ierr); 204 ierr = TSPostStage(ts,th->stage_time,0,&th->Xa);CHKERRQ(ierr); 205 ierr = TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);CHKERRQ(ierr); 206 if (!stageok) goto reject_step; 207 208 th->status = TS_STEP_PENDING; 209 ierr = VecCopy(th->X1,ts->vec_sol);CHKERRQ(ierr); 210 ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 211 th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 212 if (!accept) { 213 ierr = VecCopy(th->X0,ts->vec_sol);CHKERRQ(ierr); 214 ts->time_step = next_time_step; 215 goto reject_step; 216 } 217 218 ts->ptime += ts->time_step; 219 ts->time_step = next_time_step; 220 break; 221 222 reject_step: 223 ts->reject++; accept = PETSC_FALSE; 224 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 225 ts->reason = TS_DIVERGED_STEP_REJECTED; 226 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 227 } 228 229 } 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "TSEvaluateWLTE_Alpha" 235 static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte) 236 { 237 TS_Alpha *th = (TS_Alpha*)ts->data; 238 Vec X = th->X1; /* X = solution */ 239 Vec Y = th->vec_lte_work; /* Y = X + LTE */ 240 PetscErrorCode ierr; 241 242 PetscFunctionBegin; 243 if (ts->steprestart) { 244 /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */ 245 ierr = VecAXPY(Y,1,X);CHKERRQ(ierr); 246 } else { 247 /* Compute LTE using backward differences with non-constant time step */ 248 PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev; 249 PetscReal a = 1 + h_prev/h; 250 PetscScalar scal[3]; Vec vecs[3]; 251 scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1)); 252 vecs[0] = th->X1; vecs[1] = th->X0; vecs[2] = th->vec_sol_prev; 253 ierr = VecCopy(X,Y);CHKERRQ(ierr); 254 ierr = VecMAXPY(Y,3,scal,vecs);CHKERRQ(ierr); 255 } 256 ierr = TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte);CHKERRQ(ierr); 257 if (order) *order = 2; 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "TSRollBack_Alpha" 263 static PetscErrorCode TSRollBack_Alpha(TS ts) 264 { 265 TS_Alpha *th = (TS_Alpha*)ts->data; 266 PetscErrorCode ierr; 267 268 PetscFunctionBegin; 269 ierr = VecCopy(th->X0,ts->vec_sol);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "TSInterpolate_Alpha" 275 static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X) 276 { 277 TS_Alpha *th = (TS_Alpha*)ts->data; 278 PetscReal dt = t - ts->ptime; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 283 ierr = VecAXPY(X,th->Gamma*dt,th->V1);CHKERRQ(ierr); 284 ierr = VecAXPY(X,(1-th->Gamma)*dt,th->V0);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 #undef __FUNCT__ 289 #define __FUNCT__ "SNESTSFormFunction_Alpha" 290 static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts) 291 { 292 TS_Alpha *th = (TS_Alpha*)ts->data; 293 PetscReal ta = th->stage_time; 294 Vec Xa = th->Xa, Va = th->Va; 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 ierr = TSAlpha_StageVecs(ts,X);CHKERRQ(ierr); 299 /* F = Function(ta,Xa,Va) */ 300 ierr = TSComputeIFunction(ts,ta,Xa,Va,F,PETSC_FALSE);CHKERRQ(ierr); 301 ierr = VecScale(F,th->scale_F);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "SNESTSFormJacobian_Alpha" 307 static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts) 308 { 309 TS_Alpha *th = (TS_Alpha*)ts->data; 310 PetscReal ta = th->stage_time; 311 Vec Xa = th->Xa, Va = th->Va; 312 PetscReal dVdX = th->shift_V; 313 PetscErrorCode ierr; 314 315 PetscFunctionBegin; 316 /* J,P = Jacobian(ta,Xa,Va) */ 317 ierr = TSComputeIJacobian(ts,ta,Xa,Va,dVdX,J,P,PETSC_FALSE);CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 } 320 321 #undef __FUNCT__ 322 #define __FUNCT__ "TSReset_Alpha" 323 static PetscErrorCode TSReset_Alpha(TS ts) 324 { 325 TS_Alpha *th = (TS_Alpha*)ts->data; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 ierr = VecDestroy(&th->X0);CHKERRQ(ierr); 330 ierr = VecDestroy(&th->Xa);CHKERRQ(ierr); 331 ierr = VecDestroy(&th->X1);CHKERRQ(ierr); 332 ierr = VecDestroy(&th->V0);CHKERRQ(ierr); 333 ierr = VecDestroy(&th->Va);CHKERRQ(ierr); 334 ierr = VecDestroy(&th->V1);CHKERRQ(ierr); 335 ierr = VecDestroy(&th->vec_sol_prev);CHKERRQ(ierr); 336 ierr = VecDestroy(&th->vec_lte_work);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "TSDestroy_Alpha" 342 static PetscErrorCode TSDestroy_Alpha(TS ts) 343 { 344 PetscErrorCode ierr; 345 346 PetscFunctionBegin; 347 ierr = TSReset_Alpha(ts);CHKERRQ(ierr); 348 ierr = PetscFree(ts->data);CHKERRQ(ierr); 349 350 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaUseAdapt_C",NULL);CHKERRQ(ierr); 351 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",NULL);CHKERRQ(ierr); 352 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",NULL);CHKERRQ(ierr); 353 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",NULL);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356 357 #undef __FUNCT__ 358 #define __FUNCT__ "TSSetUp_Alpha" 359 static PetscErrorCode TSSetUp_Alpha(TS ts) 360 { 361 TS_Alpha *th = (TS_Alpha*)ts->data; 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 ierr = VecDuplicate(ts->vec_sol,&th->X0);CHKERRQ(ierr); 366 ierr = VecDuplicate(ts->vec_sol,&th->Xa);CHKERRQ(ierr); 367 ierr = VecDuplicate(ts->vec_sol,&th->X1);CHKERRQ(ierr); 368 ierr = VecDuplicate(ts->vec_sol,&th->V0);CHKERRQ(ierr); 369 ierr = VecDuplicate(ts->vec_sol,&th->Va);CHKERRQ(ierr); 370 ierr = VecDuplicate(ts->vec_sol,&th->V1);CHKERRQ(ierr); 371 372 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 373 ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 374 if (!th->adapt) { 375 ierr = TSAdaptSetType(ts->adapt,TSADAPTNONE);CHKERRQ(ierr); 376 } else { 377 ierr = VecDuplicate(ts->vec_sol,&th->vec_sol_prev);CHKERRQ(ierr); 378 ierr = VecDuplicate(ts->vec_sol,&th->vec_lte_work);CHKERRQ(ierr); 379 if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) 380 ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP; 381 } 382 383 ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr); 384 PetscFunctionReturn(0); 385 } 386 387 #undef __FUNCT__ 388 #define __FUNCT__ "TSSetFromOptions_Alpha" 389 static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts) 390 { 391 TS_Alpha *th = (TS_Alpha*)ts->data; 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 ierr = PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");CHKERRQ(ierr); 396 { 397 PetscBool flg; 398 PetscReal radius = 1; 399 PetscBool adapt = th->adapt; 400 ierr = PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlphaSetRadius",radius,&radius,&flg);CHKERRQ(ierr); 401 if (flg) {ierr = TSAlphaSetRadius(ts,radius);CHKERRQ(ierr);} 402 ierr = PetscOptionsReal("-ts_alpha_alpha_m","Algoritmic parameter alpha_m","TSAlphaSetParams",th->Alpha_m,&th->Alpha_m,NULL);CHKERRQ(ierr); 403 ierr = PetscOptionsReal("-ts_alpha_alpha_f","Algoritmic parameter alpha_f","TSAlphaSetParams",th->Alpha_f,&th->Alpha_f,NULL);CHKERRQ(ierr); 404 ierr = PetscOptionsReal("-ts_alpha_gamma","Algoritmic parameter gamma","TSAlphaSetParams",th->Gamma,&th->Gamma,NULL);CHKERRQ(ierr); 405 ierr = TSAlphaSetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma);CHKERRQ(ierr); 406 ierr = PetscOptionsBool("-ts_alpha_adapt","Use time-step adaptivity with the Alpha method","TSAlpha2UseAdapt",adapt,&adapt,&flg);CHKERRQ(ierr); 407 if (flg) {ierr = TSAlphaUseAdapt(ts,adapt);CHKERRQ(ierr);} 408 } 409 ierr = PetscOptionsTail();CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "TSView_Alpha" 415 static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer) 416 { 417 TS_Alpha *th = (TS_Alpha*)ts->data; 418 PetscBool iascii; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 423 if (iascii) { 424 ierr = PetscViewerASCIIPrintf(viewer," Alpha_m=%g, Alpha_f=%g, Gamma=%g\n",(double)th->Alpha_m,(double)th->Alpha_f,(double)th->Gamma);CHKERRQ(ierr); 425 } 426 if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);} 427 if (ts->snes) {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);} 428 PetscFunctionReturn(0); 429 } 430 431 #undef __FUNCT__ 432 #define __FUNCT__ "TSAlphaUseAdapt_Alpha" 433 static PetscErrorCode TSAlphaUseAdapt_Alpha(TS ts,PetscBool use) 434 { 435 TS_Alpha *th = (TS_Alpha*)ts->data; 436 437 PetscFunctionBegin; 438 if (use == th->adapt) PetscFunctionReturn(0); 439 if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()"); 440 th->adapt = use; 441 PetscFunctionReturn(0); 442 } 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "TSAlphaSetRadius_Alpha" 446 static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts,PetscReal radius) 447 { 448 PetscReal alpha_m,alpha_f,gamma; 449 PetscErrorCode ierr; 450 451 PetscFunctionBegin; 452 if (radius < 0 || radius > 1) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius); 453 alpha_m = (PetscReal)0.5*(3-radius)/(1+radius); 454 alpha_f = 1/(1+radius); 455 gamma = (PetscReal)0.5 + alpha_m - alpha_f; 456 ierr = TSAlphaSetParams(ts,alpha_m,alpha_f,gamma);CHKERRQ(ierr); 457 PetscFunctionReturn(0); 458 } 459 460 #undef __FUNCT__ 461 #define __FUNCT__ "TSAlphaSetParams_Alpha" 462 static PetscErrorCode TSAlphaSetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma) 463 { 464 TS_Alpha *th = (TS_Alpha*)ts->data; 465 PetscReal tol = 100*PETSC_MACHINE_EPSILON; 466 PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma; 467 468 PetscFunctionBegin; 469 th->Alpha_m = alpha_m; 470 th->Alpha_f = alpha_f; 471 th->Gamma = gamma; 472 th->order = (PetscAbsReal(res) < tol) ? 2 : 1; 473 PetscFunctionReturn(0); 474 } 475 476 #undef __FUNCT__ 477 #define __FUNCT__ "TSAlphaGetParams_Alpha" 478 static PetscErrorCode TSAlphaGetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma) 479 { 480 TS_Alpha *th = (TS_Alpha*)ts->data; 481 482 PetscFunctionBegin; 483 if (alpha_m) *alpha_m = th->Alpha_m; 484 if (alpha_f) *alpha_f = th->Alpha_f; 485 if (gamma) *gamma = th->Gamma; 486 PetscFunctionReturn(0); 487 } 488 489 /*MC 490 TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method 491 for first-order systems 492 493 Level: beginner 494 495 References: 496 K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha 497 method for integrating the filtered Navier-Stokes equations with a 498 stabilized finite element method", Computer Methods in Applied 499 Mechanics and Engineering, 190, 305-319, 2000. 500 DOI: 10.1016/S0045-7825(00)00203-6. 501 502 J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural 503 Dynamics with Improved Numerical Dissipation: The Generalized-alpha 504 Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993. 505 506 .seealso: TS, TSCreate(), TSSetType(), TSAlphaSetRadius(), TSAlphaSetParams() 507 M*/ 508 #undef __FUNCT__ 509 #define __FUNCT__ "TSCreate_Alpha" 510 PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts) 511 { 512 TS_Alpha *th; 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 ts->ops->reset = TSReset_Alpha; 517 ts->ops->destroy = TSDestroy_Alpha; 518 ts->ops->view = TSView_Alpha; 519 ts->ops->setup = TSSetUp_Alpha; 520 ts->ops->setfromoptions = TSSetFromOptions_Alpha; 521 ts->ops->step = TSStep_Alpha; 522 ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha; 523 ts->ops->rollback = TSRollBack_Alpha; 524 ts->ops->interpolate = TSInterpolate_Alpha; 525 ts->ops->snesfunction = SNESTSFormFunction_Alpha; 526 ts->ops->snesjacobian = SNESTSFormJacobian_Alpha; 527 528 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 529 ts->data = (void*)th; 530 531 th->Alpha_m = 0.5; 532 th->Alpha_f = 0.5; 533 th->Gamma = 0.5; 534 th->order = 2; 535 536 th->adapt = PETSC_FALSE; 537 538 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaUseAdapt_C",TSAlphaUseAdapt_Alpha);CHKERRQ(ierr); 539 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",TSAlphaSetRadius_Alpha);CHKERRQ(ierr); 540 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",TSAlphaSetParams_Alpha);CHKERRQ(ierr); 541 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",TSAlphaGetParams_Alpha);CHKERRQ(ierr); 542 PetscFunctionReturn(0); 543 } 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "TSAlphaUseAdapt" 547 /*@ 548 TSAlphaUseAdapt - Use time-step adaptivity with the Alpha method 549 550 Logically Collective on TS 551 552 Input Parameter: 553 + ts - timestepping context 554 - use - flag to use adaptivity 555 556 Options Database: 557 . -ts_alpha_adapt 558 559 Level: intermediate 560 561 .seealso: TSAdapt, TSADAPTBASIC 562 @*/ 563 PetscErrorCode TSAlphaUseAdapt(TS ts,PetscBool use) 564 { 565 PetscErrorCode ierr; 566 567 PetscFunctionBegin; 568 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 569 PetscValidLogicalCollectiveBool(ts,use,2); 570 ierr = PetscTryMethod(ts,"TSAlphaUseAdapt_C",(TS,PetscBool),(ts,use));CHKERRQ(ierr); 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "TSAlphaSetRadius" 576 /*@ 577 TSAlphaSetRadius - sets the desired spectral radius of the method 578 (i.e. high-frequency numerical damping) 579 580 Logically Collective on TS 581 582 The algorithmic parameters \alpha_m and \alpha_f of the 583 generalized-\alpha method can be computed in terms of a specified 584 spectral radius \rho in [0,1] for infinite time step in order to 585 control high-frequency numerical damping: 586 \alpha_m = 0.5*(3-\rho)/(1+\rho) 587 \alpha_f = 1/(1+\rho) 588 589 Input Parameter: 590 + ts - timestepping context 591 - radius - the desired spectral radius 592 593 Options Database: 594 . -ts_alpha_radius <radius> 595 596 Level: intermediate 597 598 .seealso: TSAlphaSetParams(), TSAlphaGetParams() 599 @*/ 600 PetscErrorCode TSAlphaSetRadius(TS ts,PetscReal radius) 601 { 602 PetscErrorCode ierr; 603 604 PetscFunctionBegin; 605 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 606 PetscValidLogicalCollectiveReal(ts,radius,2); 607 if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius); 608 ierr = PetscTryMethod(ts,"TSAlphaSetRadius_C",(TS,PetscReal),(ts,radius));CHKERRQ(ierr); 609 PetscFunctionReturn(0); 610 } 611 612 #undef __FUNCT__ 613 #define __FUNCT__ "TSAlphaSetParams" 614 /*@ 615 TSAlphaSetParams - sets the algorithmic parameters for TSALPHA 616 617 Logically Collective on TS 618 619 Second-order accuracy can be obtained so long as: 620 \gamma = 0.5 + alpha_m - alpha_f 621 622 Unconditional stability requires: 623 \alpha_m >= \alpha_f >= 0.5 624 625 Backward Euler method is recovered with: 626 \alpha_m = \alpha_f = gamma = 1 627 628 Input Parameter: 629 + ts - timestepping context 630 . \alpha_m - algorithmic paramenter 631 . \alpha_f - algorithmic paramenter 632 - \gamma - algorithmic paramenter 633 634 Options Database: 635 + -ts_alpha_alpha_m <alpha_m> 636 . -ts_alpha_alpha_f <alpha_f> 637 - -ts_alpha_gamma <gamma> 638 639 Note: 640 Use of this function is normally only required to hack TSALPHA to 641 use a modified integration scheme. Users should call 642 TSAlphaSetRadius() to set the desired spectral radius of the methods 643 (i.e. high-frequency damping) in order so select optimal values for 644 these parameters. 645 646 Level: advanced 647 648 .seealso: TSAlphaSetRadius(), TSAlphaGetParams() 649 @*/ 650 PetscErrorCode TSAlphaSetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma) 651 { 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 656 PetscValidLogicalCollectiveReal(ts,alpha_m,2); 657 PetscValidLogicalCollectiveReal(ts,alpha_f,3); 658 PetscValidLogicalCollectiveReal(ts,gamma,4); 659 ierr = PetscTryMethod(ts,"TSAlphaSetParams_C",(TS,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma));CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "TSAlphaGetParams" 665 /*@ 666 TSAlphaGetParams - gets the algorithmic parameters for TSALPHA 667 668 Not Collective 669 670 Input Parameter: 671 . ts - timestepping context 672 673 Output Parameters: 674 + \alpha_m - algorithmic parameter 675 . \alpha_f - algorithmic parameter 676 - \gamma - algorithmic parameter 677 678 Note: 679 Use of this function is normally only required to hack TSALPHA to 680 use a modified integration scheme. Users should call 681 TSAlphaSetRadius() to set the high-frequency damping (i.e. spectral 682 radius of the method) in order so select optimal values for these 683 parameters. 684 685 Level: advanced 686 687 .seealso: TSAlphaSetRadius(), TSAlphaSetParams() 688 @*/ 689 PetscErrorCode TSAlphaGetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma) 690 { 691 PetscErrorCode ierr; 692 693 PetscFunctionBegin; 694 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 695 if (alpha_m) PetscValidRealPointer(alpha_m,2); 696 if (alpha_f) PetscValidRealPointer(alpha_f,3); 697 if (gamma) PetscValidRealPointer(gamma,4); 698 ierr = PetscUseMethod(ts,"TSAlphaGetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma));CHKERRQ(ierr); 699 PetscFunctionReturn(0); 700 } 701