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);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);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 if (ts->problem_type == TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Only for nonlinear problems"); 140 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 141 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 /*------------------------------------------------------------*/ 145 146 #undef __FUNCT__ 147 #define __FUNCT__ "TSSetFromOptions_Theta" 148 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 149 { 150 TS_Theta *th = (TS_Theta*)ts->data; 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 ierr = PetscOptionsHead("Theta ODE solver options");CHKERRQ(ierr); 155 { 156 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 157 ierr = PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);CHKERRQ(ierr); 158 } 159 ierr = PetscOptionsTail();CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "TSView_Theta" 165 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 166 { 167 TS_Theta *th = (TS_Theta*)ts->data; 168 PetscBool iascii; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 173 if (iascii) { 174 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 175 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");CHKERRQ(ierr); 176 } else { 177 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 178 } 179 PetscFunctionReturn(0); 180 } 181 182 EXTERN_C_BEGIN 183 #undef __FUNCT__ 184 #define __FUNCT__ "TSThetaGetTheta_Theta" 185 PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta) 186 { 187 TS_Theta *th = (TS_Theta*)ts->data; 188 189 PetscFunctionBegin; 190 *theta = th->Theta; 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNCT__ 195 #define __FUNCT__ "TSThetaSetTheta_Theta" 196 PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta) 197 { 198 TS_Theta *th = (TS_Theta*)ts->data; 199 200 PetscFunctionBegin; 201 if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta); 202 th->Theta = theta; 203 PetscFunctionReturn(0); 204 } 205 EXTERN_C_END 206 207 /* ------------------------------------------------------------ */ 208 /*MC 209 TSTHETA - DAE solver using the implicit Theta method 210 211 Level: beginner 212 213 .seealso: TSCreate(), TS, TSSetType() 214 215 M*/ 216 EXTERN_C_BEGIN 217 #undef __FUNCT__ 218 #define __FUNCT__ "TSCreate_Theta" 219 PetscErrorCode TSCreate_Theta(TS ts) 220 { 221 TS_Theta *th; 222 PetscErrorCode ierr; 223 SNES snes; 224 225 PetscFunctionBegin; 226 ts->ops->reset = TSReset_Theta; 227 ts->ops->destroy = TSDestroy_Theta; 228 ts->ops->view = TSView_Theta; 229 ts->ops->setup = TSSetUp_Theta; 230 ts->ops->step = TSStep_Theta; 231 ts->ops->setfromoptions = TSSetFromOptions_Theta; 232 ts->ops->snesfunction = SNESTSFormFunction_Theta; 233 ts->ops->snesjacobian = SNESTSFormJacobian_Theta; 234 235 ts->problem_type = TS_NONLINEAR; 236 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 237 238 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 239 ts->data = (void*)th; 240 241 th->extrapolate = PETSC_FALSE; 242 th->Theta = 0.5; 243 244 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);CHKERRQ(ierr); 245 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 EXTERN_C_END 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "TSThetaGetTheta" 252 /*@ 253 TSThetaGetTheta - Get the abscissa of the stage in (0,1]. 254 255 Not Collective 256 257 Input Parameter: 258 . ts - timestepping context 259 260 Output Parameter: 261 . theta - stage abscissa 262 263 Note: 264 Use of this function is normally only required to hack TSTHETA to use a modified integration scheme. 265 266 Level: Advanced 267 268 .seealso: TSThetaSetTheta() 269 @*/ 270 PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta) 271 { 272 PetscErrorCode ierr; 273 274 PetscFunctionBegin; 275 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 276 PetscValidPointer(theta,2); 277 ierr = PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "TSThetaSetTheta" 283 /*@ 284 TSThetaSetTheta - Set the abscissa of the stage in (0,1]. 285 286 Not Collective 287 288 Input Parameter: 289 + ts - timestepping context 290 - theta - stage abscissa 291 292 Options Database: 293 . -ts_theta_theta <theta> 294 295 Level: Intermediate 296 297 .seealso: TSThetaGetTheta() 298 @*/ 299 PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta) 300 { 301 PetscErrorCode ierr; 302 303 PetscFunctionBegin; 304 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 305 ierr = PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308