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 } else { 164 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 165 } 166 PetscFunctionReturn(0); 167 } 168 169 EXTERN_C_BEGIN 170 #undef __FUNCT__ 171 #define __FUNCT__ "TSThetaGetTheta_Theta" 172 PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta) 173 { 174 TS_Theta *th = (TS_Theta*)ts->data; 175 176 PetscFunctionBegin; 177 *theta = th->Theta; 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "TSThetaSetTheta_Theta" 183 PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta) 184 { 185 TS_Theta *th = (TS_Theta*)ts->data; 186 187 PetscFunctionBegin; 188 if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta); 189 th->Theta = theta; 190 PetscFunctionReturn(0); 191 } 192 EXTERN_C_END 193 194 /* ------------------------------------------------------------ */ 195 /*MC 196 TSTHETA - DAE solver using the implicit Theta method 197 198 Level: beginner 199 200 .seealso: TSCreate(), TS, TSSetType() 201 202 M*/ 203 EXTERN_C_BEGIN 204 #undef __FUNCT__ 205 #define __FUNCT__ "TSCreate_Theta" 206 PetscErrorCode TSCreate_Theta(TS ts) 207 { 208 TS_Theta *th; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 ts->ops->reset = TSReset_Theta; 213 ts->ops->destroy = TSDestroy_Theta; 214 ts->ops->view = TSView_Theta; 215 ts->ops->setup = TSSetUp_Theta; 216 ts->ops->step = TSStep_Theta; 217 ts->ops->setfromoptions = TSSetFromOptions_Theta; 218 ts->ops->snesfunction = SNESTSFormFunction_Theta; 219 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 220 221 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 222 ts->data = (void*)th; 223 224 th->extrapolate = PETSC_FALSE; 225 th->Theta = 0.5; 226 227 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 228 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 EXTERN_C_END 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "TSThetaGetTheta" 235 /*@ 236 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 237 238 Not Collective 239 240 Input Parameter: 241 . ts - timestepping context 242 243 Output Parameter: 244 . theta - stage abscissa 245 246 Note: 247 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 248 249 Level: Advanced 250 251 .seealso: TSThetaSetTheta() 252 @*/ 253 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 254 { 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 259 PetscValidPointer(theta,2); 260 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "TSThetaSetTheta" 266 /*@ 267 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 268 269 Not Collective 270 271 Input Parameter: 272 + ts - timestepping context 273 - theta - stage abscissa 274 275 Options Database: 276 . -ts_theta_theta <theta> 277 278 Level: Intermediate 279 280 .seealso: TSThetaGetTheta() 281 @*/ 282 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 283 { 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 288 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 289 PetscFunctionReturn(0); 290 } 291 292 /* 293 * TSBEULER and TSCN are straightforward specializations of TSTHETA. 294 * The creation functions for these specializations are below. 295 */ 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "TSView_BEuler" 299 static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer) 300 { 301 PetscFunctionBegin; 302 PetscFunctionReturn(0); 303 } 304 305 /*MC 306 TSBEULER - ODE solver using the implicit backward Euler method 307 308 Level: beginner 309 310 .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA 311 312 M*/ 313 EXTERN_C_BEGIN 314 #undef __FUNCT__ 315 #define __FUNCT__ "TSCreate_BEuler" 316 PetscErrorCode TSCreate_BEuler(TS ts) 317 { 318 PetscErrorCode ierr; 319 320 PetscFunctionBegin; 321 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 322 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); 323 ts->ops->view = TSView_BEuler; 324 PetscFunctionReturn(0); 325 } 326 EXTERN_C_END 327 328 #undef __FUNCT__ 329 #define __FUNCT__ "TSView_CN" 330 static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer) 331 { 332 PetscFunctionBegin; 333 PetscFunctionReturn(0); 334 } 335 336 /*MC 337 TSCN - ODE solver using the implicit Crank-Nicolson method. 338 339 Level: beginner 340 341 Notes: 342 TSCN is equivalent to TSTHETA with Theta=0.5 343 344 .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA 345 346 M*/ 347 EXTERN_C_BEGIN 348 #undef __FUNCT__ 349 #define __FUNCT__ "TSCreate_CN" 350 PetscErrorCode TSCreate_CN(TS ts) 351 { 352 PetscErrorCode ierr; 353 354 PetscFunctionBegin; 355 ierr = TSCreate_Theta(ts);CHKERRQ(ierr); 356 ierr = TSThetaSetTheta(ts,0.5);CHKERRQ(ierr); 357 ts->ops->view = TSView_CN; 358 PetscFunctionReturn(0); 359 } 360 EXTERN_C_END 361