1 /* 2 Code for timestepping with implicit Theta method 3 */ 4 #include <private/tsimpl.h> /*I "petscts.h" I*/ 5 6 typedef struct { 7 Vec X,Xdot; /* Storage for one stage */ 8 Vec affine; /* Affine vector needed for residual at beginning of step */ 9 PetscBool extrapolate; 10 PetscBool endpoint; 11 PetscReal Theta; 12 PetscReal shift; 13 PetscReal stage_time; 14 } TS_Theta; 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "TSStep_Theta" 18 static PetscErrorCode TSStep_Theta(TS ts) 19 { 20 TS_Theta *th = (TS_Theta*)ts->data; 21 PetscInt its,lits; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ts->time_step = ts->next_time_step; 26 th->stage_time = ts->ptime + (th->endpoint ? 1. : th->Theta)*ts->time_step; 27 th->shift = 1./(th->Theta*ts->time_step); 28 29 if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */ 30 ierr = VecZeroEntries(th->Xdot);CHKERRQ(ierr); 31 if (!th->affine) {ierr = VecDuplicate(ts->vec_sol,&th->affine);CHKERRQ(ierr);} 32 ierr = TSComputeIFunction(ts,ts->ptime,ts->vec_sol,th->Xdot,th->affine,PETSC_FALSE);CHKERRQ(ierr); 33 ierr = VecScale(th->affine,(th->Theta-1.)/th->Theta);CHKERRQ(ierr); 34 } 35 if (th->extrapolate) { 36 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,ts->vec_sol);CHKERRQ(ierr); 37 } else { 38 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 39 } 40 ierr = SNESSolve(ts->snes,th->affine,th->X);CHKERRQ(ierr); 41 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 42 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 43 ts->nonlinear_its += its; ts->linear_its += lits; 44 45 if (th->endpoint) { 46 ierr = VecCopy(th->X,ts->vec_sol);CHKERRQ(ierr); 47 } else { 48 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,th->X);CHKERRQ(ierr); 49 ierr = VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 50 } 51 ts->ptime += ts->time_step; 52 ts->next_time_step = ts->time_step; 53 ts->steps++; 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "TSInterpolate_Theta" 59 static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X) 60 { 61 TS_Theta *th = (TS_Theta*)ts->data; 62 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 66 ierr = VecWAXPY(X,t-ts->ptime,th->X,th->Xdot);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 /*------------------------------------------------------------*/ 71 #undef __FUNCT__ 72 #define __FUNCT__ "TSReset_Theta" 73 static PetscErrorCode TSReset_Theta(TS ts) 74 { 75 TS_Theta *th = (TS_Theta*)ts->data; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 ierr = VecDestroy(&th->X);CHKERRQ(ierr); 80 ierr = VecDestroy(&th->Xdot);CHKERRQ(ierr); 81 ierr = VecDestroy(&th->affine);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "TSDestroy_Theta" 87 static PetscErrorCode TSDestroy_Theta(TS ts) 88 { 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = TSReset_Theta(ts);CHKERRQ(ierr); 93 ierr = PetscFree(ts->data);CHKERRQ(ierr); 94 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","",PETSC_NULL);CHKERRQ(ierr); 95 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","",PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetEndpoint_C","",PETSC_NULL);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 100 /* 101 This defines the nonlinear equation that is to be solved with SNES 102 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 103 */ 104 #undef __FUNCT__ 105 #define __FUNCT__ "SNESTSFormFunction_Theta" 106 static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts) 107 { 108 TS_Theta *th = (TS_Theta*)ts->data; 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,x);CHKERRQ(ierr); 113 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y,PETSC_FALSE);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "SNESTSFormJacobian_Theta" 119 static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts) 120 { 121 TS_Theta *th = (TS_Theta*)ts->data; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */ 126 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str,PETSC_FALSE);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "TSSetUp_Theta" 133 static PetscErrorCode TSSetUp_Theta(TS ts) 134 { 135 TS_Theta *th = (TS_Theta*)ts->data; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 140 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 /*------------------------------------------------------------*/ 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "TSSetFromOptions_Theta" 147 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 148 { 149 TS_Theta *th = (TS_Theta*)ts->data; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 ierr = PetscOptionsHead("Theta ODE solver options");CHKERRQ(ierr); 154 { 155 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 156 ierr = PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);CHKERRQ(ierr); 157 ierr = PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,PETSC_NULL);CHKERRQ(ierr); 158 } 159 ierr = PetscOptionsTail();CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "TSView_Theta" 165 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 166 { 167 TS_Theta *th = (TS_Theta*)ts->data; 168 PetscBool iascii; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 173 if (iascii) { 174 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 175 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");CHKERRQ(ierr); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 EXTERN_C_BEGIN 181 #undef __FUNCT__ 182 #define __FUNCT__ "TSThetaGetTheta_Theta" 183 PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta) 184 { 185 TS_Theta *th = (TS_Theta*)ts->data; 186 187 PetscFunctionBegin; 188 *theta = th->Theta; 189 PetscFunctionReturn(0); 190 } 191 192 #undef __FUNCT__ 193 #define __FUNCT__ "TSThetaSetTheta_Theta" 194 PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta) 195 { 196 TS_Theta *th = (TS_Theta*)ts->data; 197 198 PetscFunctionBegin; 199 if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta); 200 th->Theta = theta; 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "TSThetaSetEndpoint_Theta" 206 PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg) 207 { 208 TS_Theta *th = (TS_Theta*)ts->data; 209 210 PetscFunctionBegin; 211 th->endpoint = flg; 212 PetscFunctionReturn(0); 213 } 214 EXTERN_C_END 215 216 /* ------------------------------------------------------------ */ 217 /*MC 218 TSTHETA - DAE solver using the implicit Theta method 219 220 Level: beginner 221 222 Notes: 223 This method can be applied to DAE. 224 225 This method is cast as a 1-stage implicit Runge-Kutta method. 226 227 .vb 228 Theta | Theta 229 ------------- 230 | 1 231 .ve 232 233 For the default Theta=0.5, this is also known as the implicit midpoint rule. 234 235 When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit: 236 237 .vb 238 0 | 0 0 239 1 | 1-Theta Theta 240 ------------------- 241 | 1-Theta Theta 242 .ve 243 244 For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN). 245 246 To apply a diagonally implicit RK method to DAE, the stage formula 247 248 $ Y_i = X + h sum_j a_ij Y'_j 249 250 is interpreted as a formula for Y'_i in terms of Y_i and known stuff (Y'_j, j<i) 251 252 .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint() 253 254 M*/ 255 EXTERN_C_BEGIN 256 #undef __FUNCT__ 257 #define __FUNCT__ "TSCreate_Theta" 258 PetscErrorCode TSCreate_Theta(TS ts) 259 { 260 TS_Theta *th; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 ts->ops->reset = TSReset_Theta; 265 ts->ops->destroy = TSDestroy_Theta; 266 ts->ops->view = TSView_Theta; 267 ts->ops->setup = TSSetUp_Theta; 268 ts->ops->step = TSStep_Theta; 269 ts->ops->interpolate = TSInterpolate_Theta; 270 ts->ops->setfromoptions = TSSetFromOptions_Theta; 271 ts->ops->snesfunction = SNESTSFormFunction_Theta; 272 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 273 274 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 275 ts->data = (void*)th; 276 277 th->extrapolate = PETSC_FALSE; 278 th->Theta = 0.5; 279 280 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 281 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 282 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetEndpoint_C","TSThetaSetEndpoint_Theta",TSThetaSetEndpoint_Theta);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 EXTERN_C_END 286 287 #undef __FUNCT__ 288 #define __FUNCT__ "TSThetaGetTheta" 289 /*@ 290 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 291 292 Not Collective 293 294 Input Parameter: 295 . ts - timestepping context 296 297 Output Parameter: 298 . theta - stage abscissa 299 300 Note: 301 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 302 303 Level: Advanced 304 305 .seealso: TSThetaSetTheta() 306 @*/ 307 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 308 { 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 313 PetscValidPointer(theta,2); 314 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "TSThetaSetTheta" 320 /*@ 321 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 322 323 Not Collective 324 325 Input Parameter: 326 + ts - timestepping context 327 - theta - stage abscissa 328 329 Options Database: 330 . -ts_theta_theta <theta> 331 332 Level: Intermediate 333 334 .seealso: TSThetaGetTheta() 335 @*/ 336 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 337 { 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 342 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "TSThetaSetEndpoint" 348 /*@ 349 TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule). 350 351 Not Collective 352 353 Input Parameter: 354 + ts - timestepping context 355 - flg - PETSC_TRUE to use the endpoint variant 356 357 Options Database: 358 . -ts_theta_endpoint <flg> 359 360 Level: Intermediate 361 362 .seealso: TSTHETA, TSCN 363 @*/ 364 PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg) 365 { 366 PetscErrorCode ierr; 367 368 PetscFunctionBegin; 369 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 370 ierr = PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 371 PetscFunctionReturn(0); 372 } 373 374 /* 375 * TSBEULER and TSCN are straightforward specializations of TSTHETA. 376 * The creation functions for these specializations are below. 377 */ 378 379 #undef __FUNCT__ 380 #define __FUNCT__ "TSView_BEuler" 381 static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer) 382 { 383 PetscFunctionBegin; 384 PetscFunctionReturn(0); 385 } 386 387 /*MC 388 TSBEULER - ODE solver using the implicit backward Euler method 389 390 Level: beginner 391 392 .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA 393 394 M*/ 395 EXTERN_C_BEGIN 396 #undef __FUNCT__ 397 #define __FUNCT__ "TSCreate_BEuler" 398 PetscErrorCode TSCreate_BEuler(TS ts) 399 { 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 404 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); 405 ts->ops->view = TSView_BEuler; 406 PetscFunctionReturn(0); 407 } 408 EXTERN_C_END 409 410 #undef __FUNCT__ 411 #define __FUNCT__ "TSView_CN" 412 static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer) 413 { 414 PetscFunctionBegin; 415 PetscFunctionReturn(0); 416 } 417 418 /*MC 419 TSCN - ODE solver using the implicit Crank-Nicolson method. 420 421 Level: beginner 422 423 Notes: 424 TSCN is equivalent to TSTHETA with Theta=0.5 425 426 .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA 427 428 M*/ 429 EXTERN_C_BEGIN 430 #undef __FUNCT__ 431 #define __FUNCT__ "TSCreate_CN" 432 PetscErrorCode TSCreate_CN(TS ts) 433 { 434 PetscErrorCode ierr; 435 436 PetscFunctionBegin; 437 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 438 ierr = TSThetaSetTheta(ts,0.5);CHKERRQ(ierr); 439 ierr = TSThetaSetEndpoint(ts,PETSC_TRUE);CHKERRQ(ierr); 440 ts->ops->view = TSView_CN; 441 PetscFunctionReturn(0); 442 } 443 EXTERN_C_END 444