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 /*------------------------------------------------------------*/ 61 #undef __FUNCT__ 62 #define __FUNCT__ "TSReset_Theta" 63 static PetscErrorCode TSReset_Theta(TS ts) 64 { 65 TS_Theta *th = (TS_Theta*)ts->data; 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 ierr = VecDestroy(&th->X);CHKERRQ(ierr); 70 ierr = VecDestroy(&th->Xdot);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "TSDestroy_Theta" 76 static PetscErrorCode TSDestroy_Theta(TS ts) 77 { 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = TSReset_Theta(ts);CHKERRQ(ierr); 82 ierr = PetscFree(ts->data);CHKERRQ(ierr); 83 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","",PETSC_NULL);CHKERRQ(ierr); 84 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","",PETSC_NULL);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 /* 89 This defines the nonlinear equation that is to be solved with SNES 90 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 91 */ 92 #undef __FUNCT__ 93 #define __FUNCT__ "SNESTSFormFunction_Theta" 94 static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts) 95 { 96 TS_Theta *th = (TS_Theta*)ts->data; 97 PetscErrorCode ierr; 98 99 PetscFunctionBegin; 100 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,x);CHKERRQ(ierr); 101 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y,PETSC_FALSE);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "SNESTSFormJacobian_Theta" 107 static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts) 108 { 109 TS_Theta *th = (TS_Theta*)ts->data; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */ 114 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str,PETSC_FALSE);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "TSSetUp_Theta" 121 static PetscErrorCode TSSetUp_Theta(TS ts) 122 { 123 TS_Theta *th = (TS_Theta*)ts->data; 124 PetscErrorCode ierr; 125 126 PetscFunctionBegin; 127 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 128 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 /*------------------------------------------------------------*/ 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "TSSetFromOptions_Theta" 135 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 136 { 137 TS_Theta *th = (TS_Theta*)ts->data; 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 ierr = PetscOptionsHead("Theta ODE solver options");CHKERRQ(ierr); 142 { 143 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 144 ierr = PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);CHKERRQ(ierr); 145 } 146 ierr = PetscOptionsTail();CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "TSView_Theta" 152 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 153 { 154 TS_Theta *th = (TS_Theta*)ts->data; 155 PetscBool iascii; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 160 if (iascii) { 161 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 162 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");CHKERRQ(ierr); 163 } 164 PetscFunctionReturn(0); 165 } 166 167 EXTERN_C_BEGIN 168 #undef __FUNCT__ 169 #define __FUNCT__ "TSThetaGetTheta_Theta" 170 PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta) 171 { 172 TS_Theta *th = (TS_Theta*)ts->data; 173 174 PetscFunctionBegin; 175 *theta = th->Theta; 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "TSThetaSetTheta_Theta" 181 PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta) 182 { 183 TS_Theta *th = (TS_Theta*)ts->data; 184 185 PetscFunctionBegin; 186 if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta); 187 th->Theta = theta; 188 PetscFunctionReturn(0); 189 } 190 EXTERN_C_END 191 192 /* ------------------------------------------------------------ */ 193 /*MC 194 TSTHETA - DAE solver using the implicit Theta method 195 196 Level: beginner 197 198 .seealso: TSCreate(), TS, TSSetType() 199 200 M*/ 201 EXTERN_C_BEGIN 202 #undef __FUNCT__ 203 #define __FUNCT__ "TSCreate_Theta" 204 PetscErrorCode TSCreate_Theta(TS ts) 205 { 206 TS_Theta *th; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ts->ops->reset = TSReset_Theta; 211 ts->ops->destroy = TSDestroy_Theta; 212 ts->ops->view = TSView_Theta; 213 ts->ops->setup = TSSetUp_Theta; 214 ts->ops->step = TSStep_Theta; 215 ts->ops->setfromoptions = TSSetFromOptions_Theta; 216 ts->ops->snesfunction = SNESTSFormFunction_Theta; 217 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 218 219 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 220 ts->data = (void*)th; 221 222 th->extrapolate = PETSC_FALSE; 223 th->Theta = 0.5; 224 225 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 226 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 EXTERN_C_END 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "TSThetaGetTheta" 233 /*@ 234 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 235 236 Not Collective 237 238 Input Parameter: 239 . ts - timestepping context 240 241 Output Parameter: 242 . theta - stage abscissa 243 244 Note: 245 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 246 247 Level: Advanced 248 249 .seealso: TSThetaSetTheta() 250 @*/ 251 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 252 { 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 257 PetscValidPointer(theta,2); 258 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 #undef __FUNCT__ 263 #define __FUNCT__ "TSThetaSetTheta" 264 /*@ 265 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 266 267 Not Collective 268 269 Input Parameter: 270 + ts - timestepping context 271 - theta - stage abscissa 272 273 Options Database: 274 . -ts_theta_theta <theta> 275 276 Level: Intermediate 277 278 .seealso: TSThetaGetTheta() 279 @*/ 280 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 281 { 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 286 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 /* 291 * TSBEULER and TSCN are straightforward specializations of TSTHETA. 292 * The creation functions for these specializations are below. 293 */ 294 295 #undef __FUNCT__ 296 #define __FUNCT__ "TSView_BEuler" 297 static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer) 298 { 299 PetscFunctionBegin; 300 PetscFunctionReturn(0); 301 } 302 303 /*MC 304 TSBEULER - ODE solver using the implicit backward Euler method 305 306 Level: beginner 307 308 .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA 309 310 M*/ 311 EXTERN_C_BEGIN 312 #undef __FUNCT__ 313 #define __FUNCT__ "TSCreate_BEuler" 314 PetscErrorCode TSCreate_BEuler(TS ts) 315 { 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 320 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); 321 ts->ops->view = TSView_BEuler; 322 PetscFunctionReturn(0); 323 } 324 EXTERN_C_END 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "TSView_CN" 328 static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer) 329 { 330 PetscFunctionBegin; 331 PetscFunctionReturn(0); 332 } 333 334 /*MC 335 TSCN - ODE solver using the implicit Crank-Nicolson method. 336 337 Level: beginner 338 339 Notes: 340 TSCN is equivalent to TSTHETA with Theta=0.5 341 342 .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA 343 344 M*/ 345 EXTERN_C_BEGIN 346 #undef __FUNCT__ 347 #define __FUNCT__ "TSCreate_CN" 348 PetscErrorCode TSCreate_CN(TS ts) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 354 ierr = TSThetaSetTheta(ts,0.5);CHKERRQ(ierr); 355 ts->ops->view = TSView_CN; 356 PetscFunctionReturn(0); 357 } 358 EXTERN_C_END 359