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 Xold; 25 Vec X,Xdot; /* Storage for one stage */ 26 Vec res; /* DAE residuals */ 27 PetscTruth extrapolate; 28 PetscReal Theta; 29 PetscReal shift; 30 PetscReal stage_time; 31 } TS_Theta; 32 33 #undef __FUNCT__ 34 #define __FUNCT__ "TSStep_Theta" 35 static PetscErrorCode TSStep_Theta(TS ts,PetscInt *steps,PetscReal *ptime) 36 { 37 Vec sol = ts->vec_sol; 38 PetscErrorCode ierr; 39 PetscInt i,max_steps = ts->max_steps,its,lits; 40 TS_Theta *th = (TS_Theta*)ts->data; 41 42 PetscFunctionBegin; 43 *steps = -ts->steps; 44 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 45 46 for (i=0; i<max_steps; i++) { 47 if (ts->ptime + ts->time_step > ts->max_time) break; 48 th->stage_time = ts->ptime + th->Theta*ts->time_step; 49 th->shift = 1./(th->Theta*ts->time_step); 50 ts->ptime += ts->time_step; 51 52 ierr = VecCopy(sol,th->Xold);CHKERRQ(ierr); /* Used within function evalutaion */ 53 if (th->extrapolate) { 54 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,sol);CHKERRQ(ierr); 55 } else { 56 ierr = VecCopy(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 ierr = VecAXPY(sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 63 ts->steps++; 64 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 65 } 66 67 *steps += ts->steps; 68 *ptime = ts->ptime; 69 PetscFunctionReturn(0); 70 } 71 72 /*------------------------------------------------------------*/ 73 #undef __FUNCT__ 74 #define __FUNCT__ "TSDestroy_Theta" 75 static PetscErrorCode TSDestroy_Theta(TS ts) 76 { 77 TS_Theta *th = (TS_Theta*)ts->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = VecDestroy(th->Xold);CHKERRQ(ierr); 82 ierr = VecDestroy(th->X);CHKERRQ(ierr); 83 ierr = VecDestroy(th->Xdot);CHKERRQ(ierr); 84 ierr = VecDestroy(th->res);CHKERRQ(ierr); 85 ierr = PetscFree(th);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 /* 90 This defines the nonlinear equation that is to be solved with SNES 91 G(U) = F[t0+T*dt, U, (U-U0)*shift] = 0 92 */ 93 #undef __FUNCT__ 94 #define __FUNCT__ "TSThetaFunction" 95 static PetscErrorCode TSThetaFunction(SNES snes,Vec x,Vec y,void *ctx) 96 { 97 TS ts = (TS)ctx; 98 TS_Theta *th = (TS_Theta*)ts->data; 99 PetscErrorCode ierr; 100 101 PetscFunctionBegin; 102 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->Xold,x);CHKERRQ(ierr); 103 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "TSThetaJacobian" 109 static PetscErrorCode TSThetaJacobian(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,void *ctx) 110 { 111 TS ts = (TS)ctx; 112 TS_Theta *th = (TS_Theta*)ts->data; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 /* th->Xdot will have already been computed in TSThetaFunction */ 117 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "TSSetUp_Theta" 124 static PetscErrorCode TSSetUp_Theta(TS ts) 125 { 126 TS_Theta *th = (TS_Theta*)ts->data; 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 ierr = VecDuplicate(ts->vec_sol,&th->Xold);CHKERRQ(ierr); 131 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 132 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 133 ierr = VecDuplicate(ts->vec_sol,&th->res);CHKERRQ(ierr); 134 ierr = SNESSetFunction(ts->snes,th->res,TSThetaFunction,ts);CHKERRQ(ierr); 135 /* This is nasty. SNESSetFromOptions() is usually called in TSSetFromOptions(). With -snes_mf_operator, it will 136 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 137 context. Note that SNESSetFunction() normally has not been called before SNESSetFromOptions(), so when -snes_mf sets 138 the Jacobian user context to snes->funP, it will actually be NULL. This is not a problem because both snes->funP and 139 snes->jacP should be the TS. */ 140 { 141 Mat A,B; 142 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 143 void *ctx; 144 ierr = SNESGetJacobian(ts->snes,&A,&B,&func,&ctx);CHKERRQ(ierr); 145 ierr = SNESSetJacobian(ts->snes,A?A:ts->A,B?B:ts->B,func?func:&TSThetaJacobian,ctx?ctx:ts);CHKERRQ(ierr); 146 } 147 PetscFunctionReturn(0); 148 } 149 /*------------------------------------------------------------*/ 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "TSSetFromOptions_Theta" 153 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 154 { 155 TS_Theta *th = (TS_Theta*)ts->data; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscOptionsHead("Theta ODE solver options");CHKERRQ(ierr); 160 { 161 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 162 ierr = PetscOptionsTruth("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);CHKERRQ(ierr); 163 } 164 ierr = PetscOptionsTail();CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "TSView_Theta" 170 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 171 { 172 TS_Theta *th = (TS_Theta*)ts->data; 173 PetscTruth iascii; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 178 if (iascii) { 179 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 180 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");CHKERRQ(ierr); 181 } else { 182 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 183 } 184 PetscFunctionReturn(0); 185 } 186 187 /* ------------------------------------------------------------ */ 188 /*MC 189 TSTHETA - DAE solver using the implicit Theta method 190 191 Level: beginner 192 193 .seealso: TSCreate(), TS, TSSetType() 194 195 M*/ 196 EXTERN_C_BEGIN 197 #undef __FUNCT__ 198 #define __FUNCT__ "TSCreate_Theta" 199 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Theta(TS ts) 200 { 201 TS_Theta *th; 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 ts->ops->destroy = TSDestroy_Theta; 206 ts->ops->view = TSView_Theta; 207 ts->ops->setup = TSSetUp_Theta; 208 ts->ops->step = TSStep_Theta; 209 ts->ops->setfromoptions = TSSetFromOptions_Theta; 210 211 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 212 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 213 214 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 215 ts->data = (void*)th; 216 217 th->extrapolate = PETSC_TRUE; 218 th->Theta = 0.5; 219 220 PetscFunctionReturn(0); 221 } 222 EXTERN_C_END 223