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 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 Vec sol = ts->vec_sol; 37 PetscErrorCode ierr; 38 PetscInt i,max_steps = ts->max_steps,its,lits; 39 TS_Theta *th = (TS_Theta*)ts->data; 40 41 PetscFunctionBegin; 42 *steps = -ts->steps; 43 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 44 45 for (i=0; i<max_steps; i++) { 46 if (ts->ptime + ts->time_step > ts->max_time) break; 47 th->stage_time = ts->ptime + th->Theta*ts->time_step; 48 th->shift = 1./(th->Theta*ts->time_step); 49 ts->ptime += ts->time_step; 50 51 ierr = VecCopy(sol,th->Xold);CHKERRQ(ierr); /* Used within function evalutaion */ 52 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,sol);CHKERRQ(ierr); 53 ierr = SNESSolve(ts->snes,PETSC_NULL,th->X);CHKERRQ(ierr); 54 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 55 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 56 ts->nonlinear_its += its; ts->linear_its += lits; 57 ierr = VecAXPY(sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 58 ts->steps++; 59 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 60 } 61 62 *steps += ts->steps; 63 *ptime = ts->ptime; 64 PetscFunctionReturn(0); 65 } 66 67 /*------------------------------------------------------------*/ 68 #undef __FUNCT__ 69 #define __FUNCT__ "TSDestroy_Theta" 70 static PetscErrorCode TSDestroy_Theta(TS ts) 71 { 72 TS_Theta *th = (TS_Theta*)ts->data; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 ierr = VecDestroy(th->Xold);CHKERRQ(ierr); 77 ierr = VecDestroy(th->X);CHKERRQ(ierr); 78 ierr = VecDestroy(th->Xdot);CHKERRQ(ierr); 79 ierr = VecDestroy(th->res);CHKERRQ(ierr); 80 ierr = PetscFree(th);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 /* 85 This defines the nonlinear equation that is to be solved with SNES 86 G(U) = F[t0+T*dt, U, (U-U0)*shift] = 0 87 */ 88 #undef __FUNCT__ 89 #define __FUNCT__ "TSThetaFunction" 90 static PetscErrorCode TSThetaFunction(SNES snes,Vec x,Vec y,void *ctx) 91 { 92 TS ts = (TS)ctx; 93 TS_Theta *th = (TS_Theta*)ts->data; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->Xold,x);CHKERRQ(ierr); 98 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "TSThetaJacobian" 104 static PetscErrorCode TSThetaJacobian(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,void *ctx) 105 { 106 TS ts = (TS)ctx; 107 TS_Theta *th = (TS_Theta*)ts->data; 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 /* th->Xdot will have already been computed in TSThetaFunction */ 112 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "TSSetUp_Theta" 119 static PetscErrorCode TSSetUp_Theta(TS ts) 120 { 121 TS_Theta *th = (TS_Theta*)ts->data; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 ierr = VecDuplicate(ts->vec_sol,&th->Xold);CHKERRQ(ierr); 126 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 127 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 128 ierr = VecDuplicate(ts->vec_sol,&th->res);CHKERRQ(ierr); 129 ierr = SNESSetFunction(ts->snes,th->res,TSThetaFunction,ts);CHKERRQ(ierr); 130 /* This is nasty. SNESSetFromOptions() is usually called in TSSetFromOptions(). With -snes_mf_operator, it will 131 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 132 context. Note that SNESSetFunction() normally has not been called before SNESSetFromOptions(), so when -snes_mf sets 133 the Jacobian user context to snes->funP, it will actually be NULL. This is not a problem because both snes->funP and 134 snes->jacP should be the TS. */ 135 { 136 Mat A,B; 137 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 138 void *ctx; 139 ierr = SNESGetJacobian(ts->snes,&A,&B,&func,&ctx);CHKERRQ(ierr); 140 ierr = SNESSetJacobian(ts->snes,A?A:ts->A,B?B:ts->B,func?func:&TSThetaJacobian,ctx?ctx:ts);CHKERRQ(ierr); 141 } 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("SUNDIALS 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 } 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 PetscTruth iascii; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 172 if (iascii) { 173 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 174 } else { 175 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 /* ------------------------------------------------------------ */ 181 /*MC 182 TS_THETA - DAE solver using the implicit Theta method 183 184 Level: beginner 185 186 .seealso: TSCreate(), TS, TSSetType() 187 188 M*/ 189 EXTERN_C_BEGIN 190 #undef __FUNCT__ 191 #define __FUNCT__ "TSCreate_Theta" 192 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Theta(TS ts) 193 { 194 TS_Theta *th; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 ts->ops->destroy = TSDestroy_Theta; 199 ts->ops->view = TSView_Theta; 200 ts->ops->setup = TSSetUp_Theta; 201 ts->ops->step = TSStep_Theta; 202 ts->ops->setfromoptions = TSSetFromOptions_Theta; 203 204 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 205 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 206 207 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 208 ts->data = (void*)th; 209 210 th->Theta = 0.5; 211 212 PetscFunctionReturn(0); 213 } 214 EXTERN_C_END 215