1 /* 2 Code for timestepping with implicit Theta method 3 4 Notes: 5 This method can be applied to DAE. 6 7 This method is cast as a 1-stage implicit Runge-Kutta method. 8 9 Theta | Theta 10 ------------- 11 | 1 12 13 To apply a diagonally implicit RK method to DAE, the stage formula 14 15 X_i = x + h sum_j a_ij X'_j 16 17 is interpreted as a formula for X'_i in terms of X_i and known stuff (X'_j, j<i) 18 */ 19 #include <private/tsimpl.h> /*I "petscts.h" I*/ 20 21 typedef struct { 22 Vec X,Xdot; /* Storage for one stage */ 23 PetscBool extrapolate; 24 PetscReal Theta; 25 PetscReal shift; 26 PetscReal stage_time; 27 } TS_Theta; 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "TSStep_Theta" 31 static PetscErrorCode TSStep_Theta(TS ts) 32 { 33 TS_Theta *th = (TS_Theta*)ts->data; 34 PetscInt its,lits; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 ts->time_step = ts->next_time_step; 39 th->stage_time = ts->ptime + th->Theta*ts->time_step; 40 th->shift = 1./(th->Theta*ts->time_step); 41 42 if (th->extrapolate) { 43 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,ts->vec_sol);CHKERRQ(ierr); 44 } else { 45 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 46 } 47 ierr = SNESSolve(ts->snes,PETSC_NULL,th->X);CHKERRQ(ierr); 48 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 49 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 50 ts->nonlinear_its += its; ts->linear_its += lits; 51 52 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,th->X);CHKERRQ(ierr); 53 ierr = VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 54 ts->ptime += ts->time_step; 55 ts->next_time_step = ts->time_step; 56 ts->steps++; 57 PetscFunctionReturn(0); 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "TSInterpolate_Theta" 62 static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X) 63 { 64 TS_Theta *th = (TS_Theta*)ts->data; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 ierr = VecWAXPY(X,t-ts->ptime,ts->vec_sol,th->Xdot);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 /*------------------------------------------------------------*/ 73 #undef __FUNCT__ 74 #define __FUNCT__ "TSReset_Theta" 75 static PetscErrorCode TSReset_Theta(TS ts) 76 { 77 TS_Theta *th = (TS_Theta*)ts->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = VecDestroy(&th->X);CHKERRQ(ierr); 82 ierr = VecDestroy(&th->Xdot);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "TSDestroy_Theta" 88 static PetscErrorCode TSDestroy_Theta(TS ts) 89 { 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = TSReset_Theta(ts);CHKERRQ(ierr); 94 ierr = PetscFree(ts->data);CHKERRQ(ierr); 95 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","",PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_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 } 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 EXTERN_C_END 203 204 /* ------------------------------------------------------------ */ 205 /*MC 206 TSTHETA - DAE solver using the implicit Theta method 207 208 Level: beginner 209 210 .seealso: TSCreate(), TS, TSSetType() 211 212 M*/ 213 EXTERN_C_BEGIN 214 #undef __FUNCT__ 215 #define __FUNCT__ "TSCreate_Theta" 216 PetscErrorCode TSCreate_Theta(TS ts) 217 { 218 TS_Theta *th; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 ts->ops->reset = TSReset_Theta; 223 ts->ops->destroy = TSDestroy_Theta; 224 ts->ops->view = TSView_Theta; 225 ts->ops->setup = TSSetUp_Theta; 226 ts->ops->step = TSStep_Theta; 227 ts->ops->interpolate = TSInterpolate_Theta; 228 ts->ops->setfromoptions = TSSetFromOptions_Theta; 229 ts->ops->snesfunction = SNESTSFormFunction_Theta; 230 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 231 232 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 233 ts->data = (void*)th; 234 235 th->extrapolate = PETSC_FALSE; 236 th->Theta = 0.5; 237 238 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 239 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 EXTERN_C_END 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "TSThetaGetTheta" 246 /*@ 247 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 248 249 Not Collective 250 251 Input Parameter: 252 . ts - timestepping context 253 254 Output Parameter: 255 . theta - stage abscissa 256 257 Note: 258 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 259 260 Level: Advanced 261 262 .seealso: TSThetaSetTheta() 263 @*/ 264 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 265 { 266 PetscErrorCode ierr; 267 268 PetscFunctionBegin; 269 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 270 PetscValidPointer(theta,2); 271 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "TSThetaSetTheta" 277 /*@ 278 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 279 280 Not Collective 281 282 Input Parameter: 283 + ts - timestepping context 284 - theta - stage abscissa 285 286 Options Database: 287 . -ts_theta_theta <theta> 288 289 Level: Intermediate 290 291 .seealso: TSThetaGetTheta() 292 @*/ 293 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 294 { 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 299 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 /* 304 * TSBEULER and TSCN are straightforward specializations of TSTHETA. 305 * The creation functions for these specializations are below. 306 */ 307 308 #undef __FUNCT__ 309 #define __FUNCT__ "TSView_BEuler" 310 static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer) 311 { 312 PetscFunctionBegin; 313 PetscFunctionReturn(0); 314 } 315 316 /*MC 317 TSBEULER - ODE solver using the implicit backward Euler method 318 319 Level: beginner 320 321 .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA 322 323 M*/ 324 EXTERN_C_BEGIN 325 #undef __FUNCT__ 326 #define __FUNCT__ "TSCreate_BEuler" 327 PetscErrorCode TSCreate_BEuler(TS ts) 328 { 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 333 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); 334 ts->ops->view = TSView_BEuler; 335 PetscFunctionReturn(0); 336 } 337 EXTERN_C_END 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "TSView_CN" 341 static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer) 342 { 343 PetscFunctionBegin; 344 PetscFunctionReturn(0); 345 } 346 347 /*MC 348 TSCN - ODE solver using the implicit Crank-Nicolson method. 349 350 Level: beginner 351 352 Notes: 353 TSCN is equivalent to TSTHETA with Theta=0.5 354 355 .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA 356 357 M*/ 358 EXTERN_C_BEGIN 359 #undef __FUNCT__ 360 #define __FUNCT__ "TSCreate_CN" 361 PetscErrorCode TSCreate_CN(TS ts) 362 { 363 PetscErrorCode ierr; 364 365 PetscFunctionBegin; 366 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 367 ierr = TSThetaSetTheta(ts,0.5);CHKERRQ(ierr); 368 ts->ops->view = TSView_CN; 369 PetscFunctionReturn(0); 370 } 371 EXTERN_C_END 372