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