1 #define PETSCTS_DLL 2 3 /* 4 Code for timestepping with implicit Theta method 5 6 Notes: 7 This method can be applied to DAE. 8 9 This method is cast as a 1-stage implicit Runge-Kutta method. 10 11 Theta | Theta 12 ------------- 13 | 1 14 15 To apply a diagonally implicit RK method to DAE, the stage formula 16 17 X_i = x + h sum_j a_ij X'_j 18 19 is interpreted as a formula for X'_i in terms of X_i and known stuff (X'_j, j<i) 20 */ 21 #include "private/tsimpl.h" /*I "petscts.h" I*/ 22 23 typedef struct { 24 Vec X,Xdot; /* Storage for one stage */ 25 Vec res; /* DAE residuals */ 26 PetscTruth extrapolate; 27 PetscReal Theta; 28 PetscReal shift; 29 PetscReal stage_time; 30 } TS_Theta; 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "TSStep_Theta" 34 static PetscErrorCode TSStep_Theta(TS ts,PetscInt *steps,PetscReal *ptime) 35 { 36 TS_Theta *th = (TS_Theta*)ts->data; 37 PetscInt i,max_steps = ts->max_steps,its,lits; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 *steps = -ts->steps; 42 *ptime = ts->ptime; 43 44 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 45 46 for (i=0; i<max_steps; i++) { 47 if (ts->ptime + ts->time_step > ts->max_time) break; 48 ierr = TSPreStep(ts);CHKERRQ(ierr); 49 50 th->stage_time = ts->ptime + th->Theta*ts->time_step; 51 th->shift = 1./(th->Theta*ts->time_step); 52 53 if (th->extrapolate) { 54 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,ts->vec_sol);CHKERRQ(ierr); 55 } else { 56 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 57 } 58 ierr = SNESSolve(ts->snes,PETSC_NULL,th->X);CHKERRQ(ierr); 59 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 60 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 61 ts->nonlinear_its += its; ts->linear_its += lits; 62 63 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,th->X);CHKERRQ(ierr); 64 ierr = VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 65 ts->ptime += ts->time_step; 66 ts->steps++; 67 68 ierr = TSPostStep(ts);CHKERRQ(ierr); 69 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 70 } 71 72 *steps += ts->steps; 73 *ptime = ts->ptime; 74 PetscFunctionReturn(0); 75 } 76 77 /*------------------------------------------------------------*/ 78 #undef __FUNCT__ 79 #define __FUNCT__ "TSDestroy_Theta" 80 static PetscErrorCode TSDestroy_Theta(TS ts) 81 { 82 TS_Theta *th = (TS_Theta*)ts->data; 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 if (th->X) {ierr = VecDestroy(th->X);CHKERRQ(ierr);} 87 if (th->Xdot) {ierr = VecDestroy(th->Xdot);CHKERRQ(ierr);} 88 if (th->res) {ierr = VecDestroy(th->res);CHKERRQ(ierr);} 89 ierr = PetscFree(th);CHKERRQ(ierr); 90 PetscFunctionReturn(0); 91 } 92 93 /* 94 This defines the nonlinear equation that is to be solved with SNES 95 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 96 */ 97 #undef __FUNCT__ 98 #define __FUNCT__ "TSThetaFunction" 99 static PetscErrorCode TSThetaFunction(SNES snes,Vec x,Vec y,void *ctx) 100 { 101 TS ts = (TS)ctx; 102 TS_Theta *th = (TS_Theta*)ts->data; 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,x);CHKERRQ(ierr); 107 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "TSThetaJacobian" 113 static PetscErrorCode TSThetaJacobian(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,void *ctx) 114 { 115 TS ts = (TS)ctx; 116 TS_Theta *th = (TS_Theta*)ts->data; 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 /* th->Xdot has already been computed in TSThetaFunction (SNES guarantees this) */ 121 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "TSSetUp_Theta" 128 static PetscErrorCode TSSetUp_Theta(TS ts) 129 { 130 TS_Theta *th = (TS_Theta*)ts->data; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 if (ts->problem_type == TS_LINEAR) { 135 SETERRQ(PETSC_ERR_ARG_WRONG,"Only for nonlinear problems"); 136 } 137 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 138 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 139 ierr = VecDuplicate(ts->vec_sol,&th->res);CHKERRQ(ierr); 140 ierr = SNESSetFunction(ts->snes,th->res,TSThetaFunction,ts);CHKERRQ(ierr); 141 /* This is nasty. SNESSetFromOptions() is usually called in TSSetFromOptions(). With -snes_mf_operator, it will 142 replace A and we don't want to mess with that. With -snes_mf, A and B will be replaced as well as the function and 143 context. Note that SNESSetFunction() normally has not been called before SNESSetFromOptions(), so when -snes_mf sets 144 the Jacobian user context to snes->funP, it will actually be NULL. This is not a problem because both snes->funP and 145 snes->jacP should be the TS. */ 146 { 147 Mat A,B; 148 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 149 void *ctx; 150 ierr = SNESGetJacobian(ts->snes,&A,&B,&func,&ctx);CHKERRQ(ierr); 151 ierr = SNESSetJacobian(ts->snes,A?A:ts->A,B?B:ts->B,func?func:&TSThetaJacobian,ctx?ctx:ts);CHKERRQ(ierr); 152 } 153 PetscFunctionReturn(0); 154 } 155 /*------------------------------------------------------------*/ 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "TSSetFromOptions_Theta" 159 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 160 { 161 TS_Theta *th = (TS_Theta*)ts->data; 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 ierr = PetscOptionsHead("Theta ODE solver options");CHKERRQ(ierr); 166 { 167 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 168 ierr = PetscOptionsTruth("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);CHKERRQ(ierr); 169 } 170 ierr = PetscOptionsTail();CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "TSView_Theta" 176 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 177 { 178 TS_Theta *th = (TS_Theta*)ts->data; 179 PetscTruth iascii; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 184 if (iascii) { 185 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 186 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");CHKERRQ(ierr); 187 } else { 188 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 189 } 190 PetscFunctionReturn(0); 191 } 192 193 /* ------------------------------------------------------------ */ 194 /*MC 195 TSTHETA - DAE solver using the implicit Theta method 196 197 Level: beginner 198 199 .seealso: TSCreate(), TS, TSSetType() 200 201 M*/ 202 EXTERN_C_BEGIN 203 #undef __FUNCT__ 204 #define __FUNCT__ "TSCreate_Theta" 205 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Theta(TS ts) 206 { 207 TS_Theta *th; 208 PetscErrorCode ierr; 209 210 PetscFunctionBegin; 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 217 ts->problem_type = TS_NONLINEAR; 218 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 219 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 220 221 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 222 ts->data = (void*)th; 223 224 th->extrapolate = PETSC_TRUE; 225 th->Theta = 0.5; 226 227 PetscFunctionReturn(0); 228 } 229 EXTERN_C_END 230