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