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,PetscInt *steps,PetscReal *ptime) 32 { 33 TS_Theta *th = (TS_Theta*)ts->data; 34 PetscInt i,its,lits; 35 PetscErrorCode ierr; 36 37 PetscFunctionBegin; 38 *steps = -ts->steps; 39 *ptime = ts->ptime; 40 41 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 42 43 for (i=0; i<ts->max_steps; i++) { 44 if (ts->ptime + ts->time_step > ts->max_time) break; 45 ierr = TSPreStep(ts);CHKERRQ(ierr); 46 47 th->stage_time = ts->ptime + th->Theta*ts->time_step; 48 th->shift = 1./(th->Theta*ts->time_step); 49 50 if (th->extrapolate) { 51 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,ts->vec_sol);CHKERRQ(ierr); 52 } else { 53 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 54 } 55 ierr = SNESSolve(ts->snes,PETSC_NULL,th->X);CHKERRQ(ierr); 56 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 57 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 58 ts->nonlinear_its += its; ts->linear_its += lits; 59 60 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,th->X);CHKERRQ(ierr); 61 ierr = VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 62 ts->ptime += ts->time_step; 63 ts->steps++; 64 65 ierr = TSPostStep(ts);CHKERRQ(ierr); 66 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 67 } 68 69 *steps += ts->steps; 70 *ptime = ts->ptime; 71 PetscFunctionReturn(0); 72 } 73 74 /*------------------------------------------------------------*/ 75 #undef __FUNCT__ 76 #define __FUNCT__ "TSReset_Theta" 77 static PetscErrorCode TSReset_Theta(TS ts) 78 { 79 TS_Theta *th = (TS_Theta*)ts->data; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 ierr = VecDestroy(&th->X);CHKERRQ(ierr); 84 ierr = VecDestroy(&th->Xdot);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "TSDestroy_Theta" 90 static PetscErrorCode TSDestroy_Theta(TS ts) 91 { 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = TSReset_Theta(ts);CHKERRQ(ierr); 96 ierr = PetscFree(ts->data);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 } else { 176 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 177 } 178 PetscFunctionReturn(0); 179 } 180 181 EXTERN_C_BEGIN 182 #undef __FUNCT__ 183 #define __FUNCT__ "TSThetaGetTheta_Theta" 184 PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta) 185 { 186 TS_Theta *th = (TS_Theta*)ts->data; 187 188 PetscFunctionBegin; 189 *theta = th->Theta; 190 PetscFunctionReturn(0); 191 } 192 193 #undef __FUNCT__ 194 #define __FUNCT__ "TSThetaSetTheta_Theta" 195 PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta) 196 { 197 TS_Theta *th = (TS_Theta*)ts->data; 198 199 PetscFunctionBegin; 200 if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta); 201 th->Theta = theta; 202 PetscFunctionReturn(0); 203 } 204 EXTERN_C_END 205 206 /* ------------------------------------------------------------ */ 207 /*MC 208 TSTHETA - DAE solver using the implicit Theta method 209 210 Level: beginner 211 212 .seealso: TSCreate(), TS, TSSetType() 213 214 M*/ 215 EXTERN_C_BEGIN 216 #undef __FUNCT__ 217 #define __FUNCT__ "TSCreate_Theta" 218 PetscErrorCode TSCreate_Theta(TS ts) 219 { 220 TS_Theta *th; 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 ts->ops->reset = TSReset_Theta; 225 ts->ops->destroy = TSDestroy_Theta; 226 ts->ops->view = TSView_Theta; 227 ts->ops->setup = TSSetUp_Theta; 228 ts->ops->step = TSStep_Theta; 229 ts->ops->setfromoptions = TSSetFromOptions_Theta; 230 ts->ops->snesfunction = SNESTSFormFunction_Theta; 231 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 232 233 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 234 ts->data = (void*)th; 235 236 th->extrapolate = PETSC_FALSE; 237 th->Theta = 0.5; 238 239 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 240 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 EXTERN_C_END 244 245 #undef __FUNCT__ 246 #define __FUNCT__ "TSThetaGetTheta" 247 /*@ 248 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 249 250 Not Collective 251 252 Input Parameter: 253 . ts - timestepping context 254 255 Output Parameter: 256 . theta - stage abscissa 257 258 Note: 259 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 260 261 Level: Advanced 262 263 .seealso: TSThetaSetTheta() 264 @*/ 265 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 266 { 267 PetscErrorCode ierr; 268 269 PetscFunctionBegin; 270 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 271 PetscValidPointer(theta,2); 272 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "TSThetaSetTheta" 278 /*@ 279 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 280 281 Not Collective 282 283 Input Parameter: 284 + ts - timestepping context 285 - theta - stage abscissa 286 287 Options Database: 288 . -ts_theta_theta <theta> 289 290 Level: Intermediate 291 292 .seealso: TSThetaGetTheta() 293 @*/ 294 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 295 { 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 300 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303