1 /* 2 Code for Timestepping with explicit SSP. 3 */ 4 #include <private/tsimpl.h> /*I "petscts.h" I*/ 5 6 PetscFList TSSSPList = 0; 7 #define TSSSPType char* 8 9 #define TSSSPRKS2 "rks2" 10 #define TSSSPRKS3 "rks3" 11 #define TSSSPRK104 "rk104" 12 13 typedef struct { 14 PetscErrorCode (*onestep)(TS,PetscReal,PetscReal,Vec); 15 PetscInt nstages; 16 Vec *work; 17 PetscInt nwork; 18 PetscBool workout; 19 } TS_SSP; 20 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "SSPGetWorkVectors" 24 static PetscErrorCode SSPGetWorkVectors(TS ts,PetscInt n,Vec **work) 25 { 26 TS_SSP *ssp = (TS_SSP*)ts->data; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten"); 31 if (ssp->nwork < n) { 32 if (ssp->nwork > 0) { 33 ierr = VecDestroyVecs(ssp->nwork,&ssp->work);CHKERRQ(ierr); 34 } 35 ierr = VecDuplicateVecs(ts->vec_sol,n,&ssp->work);CHKERRQ(ierr); 36 ssp->nwork = n; 37 } 38 *work = ssp->work; 39 ssp->workout = PETSC_TRUE; 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "SSPRestoreWorkVectors" 45 static PetscErrorCode SSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work) 46 { 47 TS_SSP *ssp = (TS_SSP*)ts->data; 48 49 PetscFunctionBegin; 50 if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten"); 51 if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out"); 52 ssp->workout = PETSC_FALSE; 53 *work = PETSC_NULL; 54 PetscFunctionReturn(0); 55 } 56 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "SSPStep_RK_2" 60 /* Optimal second order SSP Runge-Kutta, low-storage, c_eff=(s-1)/s */ 61 /* Pseudocode 2 of Ketcheson 2008 */ 62 static PetscErrorCode SSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol) 63 { 64 TS_SSP *ssp = (TS_SSP*)ts->data; 65 Vec *work,F; 66 PetscInt i,s; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 s = ssp->nstages; 71 ierr = SSPGetWorkVectors(ts,2,&work);CHKERRQ(ierr); 72 F = work[1]; 73 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 74 for (i=0; i<s-1; i++) { 75 ierr = TSComputeRHSFunction(ts,t0+dt*(i/(s-1.)),work[0],F);CHKERRQ(ierr); 76 ierr = VecAXPY(work[0],dt/(s-1.),F);CHKERRQ(ierr); 77 } 78 ierr = TSComputeRHSFunction(ts,t0+dt,work[0],F);CHKERRQ(ierr); 79 ierr = VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);CHKERRQ(ierr); 80 ierr = SSPRestoreWorkVectors(ts,2,&work);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "SSPStep_RK_3" 86 /* Optimal third order SSP Runge-Kutta, low-storage, c_eff=(sqrt(s)-1)/sqrt(s), where sqrt(s) is an integer */ 87 /* Pseudocode 2 of Ketcheson 2008 */ 88 static PetscErrorCode SSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol) 89 { 90 TS_SSP *ssp = (TS_SSP*)ts->data; 91 Vec *work,F; 92 PetscInt i,s,n,r; 93 PetscReal c; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 s = ssp->nstages; 98 n = (PetscInt)(sqrt((PetscReal)s)+0.001); 99 r = s-n; 100 if (n*n != s) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for optimal third order schemes with %d stages, must be a square number at least 4",s); 101 ierr = SSPGetWorkVectors(ts,3,&work);CHKERRQ(ierr); 102 F = work[2]; 103 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 104 for (i=0; i<(n-1)*(n-2)/2; i++) { 105 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 106 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 107 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 108 } 109 ierr = VecCopy(work[0],work[1]);CHKERRQ(ierr); 110 for ( ; i<n*(n+1)/2-1; i++) { 111 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 112 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 113 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 114 } 115 { 116 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 117 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 118 ierr = VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);CHKERRQ(ierr); 119 i++; 120 } 121 for ( ; i<s; i++) { 122 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 123 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 124 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 125 } 126 ierr = VecCopy(work[0],sol);CHKERRQ(ierr); 127 ierr = SSPRestoreWorkVectors(ts,3,&work);CHKERRQ(ierr); 128 PetscFunctionReturn(0); 129 } 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "SSPStep_RK_10_4" 133 /* Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6 */ 134 /* SSPRK(10,4), Pseudocode 3 of Ketcheson 2008 */ 135 static PetscErrorCode SSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol) 136 { 137 const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1}; 138 Vec *work,F; 139 PetscInt i; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 ierr = SSPGetWorkVectors(ts,3,&work);CHKERRQ(ierr); 144 F = work[2]; 145 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 146 for (i=0; i<5; i++) { 147 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 148 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 149 } 150 ierr = VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);CHKERRQ(ierr); 151 ierr = VecAXPBY(work[0],15,-5,work[1]);CHKERRQ(ierr); 152 for ( ; i<9; i++) { 153 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 154 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 155 } 156 ierr = TSComputeRHSFunction(ts,t0+dt,work[0],F);CHKERRQ(ierr); 157 ierr = VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);CHKERRQ(ierr); 158 ierr = VecCopy(work[1],sol);CHKERRQ(ierr); 159 ierr = SSPRestoreWorkVectors(ts,3,&work);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "TSSetUp_SSP" 166 static PetscErrorCode TSSetUp_SSP(TS ts) 167 { 168 /* TS_SSP *ssp = (TS_SSP*)ts->data; */ 169 /* PetscErrorCode ierr; */ 170 171 PetscFunctionBegin; 172 PetscFunctionReturn(0); 173 } 174 175 #undef __FUNCT__ 176 #define __FUNCT__ "TSStep_SSP" 177 static PetscErrorCode TSStep_SSP(TS ts) 178 { 179 TS_SSP *ssp = (TS_SSP*)ts->data; 180 Vec sol = ts->vec_sol; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);CHKERRQ(ierr); 185 ts->ptime += ts->time_step; 186 ts->steps++; 187 PetscFunctionReturn(0); 188 } 189 /*------------------------------------------------------------*/ 190 #undef __FUNCT__ 191 #define __FUNCT__ "TSReset_SSP" 192 static PetscErrorCode TSReset_SSP(TS ts) 193 { 194 TS_SSP *ssp = (TS_SSP*)ts->data; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 if (ssp->work) {ierr = VecDestroyVecs(ssp->nwork,&ssp->work);CHKERRQ(ierr);} 199 ssp->nwork = 0; 200 ssp->workout = PETSC_FALSE; 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "TSDestroy_SSP" 206 static PetscErrorCode TSDestroy_SSP(TS ts) 207 { 208 PetscErrorCode ierr; 209 210 PetscFunctionBegin; 211 ierr = TSReset_SSP(ts);CHKERRQ(ierr); 212 ierr = PetscFree(ts->data);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 /*------------------------------------------------------------*/ 216 217 #undef __FUNCT__ 218 #define __FUNCT__ "TSSSPSetType" 219 static PetscErrorCode TSSSPSetType(TS ts,const TSSSPType type) 220 { 221 PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec); 222 TS_SSP *ssp = (TS_SSP*)ts->data; 223 224 PetscFunctionBegin; 225 ierr = PetscFListFind(TSSSPList,((PetscObject)ts)->comm,type,PETSC_TRUE,(void(**)(void))&r);CHKERRQ(ierr); 226 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type); 227 ssp->onestep = r; 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "TSSetFromOptions_SSP" 233 static PetscErrorCode TSSetFromOptions_SSP(TS ts) 234 { 235 char tname[256] = TSSSPRKS2; 236 TS_SSP *ssp = (TS_SSP*)ts->data; 237 PetscErrorCode ierr; 238 PetscBool flg; 239 240 PetscFunctionBegin; 241 ierr = PetscOptionsHead("SSP ODE solver options");CHKERRQ(ierr); 242 { 243 ierr = PetscOptionsList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);CHKERRQ(ierr); 244 if (flg) { 245 ierr = TSSSPSetType(ts,tname);CHKERRQ(ierr); 246 } 247 ierr = PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,PETSC_NULL);CHKERRQ(ierr); 248 } 249 ierr = PetscOptionsTail();CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "TSView_SSP" 255 static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer) 256 { 257 PetscFunctionBegin; 258 PetscFunctionReturn(0); 259 } 260 261 /* ------------------------------------------------------------ */ 262 263 /*MC 264 TSSSP - Explicit strong stability preserving ODE solver 265 266 Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation 267 bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's 268 scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties, 269 but they are usually formulated using a forward Euler time discretization or by coupling the space and time 270 discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very 271 difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the 272 semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a 273 time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP). 274 275 Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while 276 still being SSP. Some theoretical bounds 277 278 1. There are no explicit methods with c_eff > 1. 279 280 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0. 281 282 3. There are no implicit methods with order greater than 1 and c_eff > 2. 283 284 This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows 285 for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work 286 vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice. 287 288 Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104} 289 290 rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s 291 292 rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s 293 294 rk104: A 10-stage fourth order method. c_eff = 0.6 295 296 Level: beginner 297 298 References: 299 Ketcheson, Highly efficient strong stability preserving Runge-Kutta methods with low-storage implementations, SISC, 2008. 300 301 Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009. 302 303 .seealso: TSCreate(), TS, TSSetType() 304 305 M*/ 306 EXTERN_C_BEGIN 307 #undef __FUNCT__ 308 #define __FUNCT__ "TSCreate_SSP" 309 PetscErrorCode TSCreate_SSP(TS ts) 310 { 311 TS_SSP *ssp; 312 PetscErrorCode ierr; 313 314 PetscFunctionBegin; 315 if (!TSSSPList) { 316 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS2, "SSPStep_RK_2", (void(*)(void))SSPStep_RK_2);CHKERRQ(ierr); 317 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS3, "SSPStep_RK_3", (void(*)(void))SSPStep_RK_3);CHKERRQ(ierr); 318 ierr = PetscFListAdd(&TSSSPList,TSSSPRK104, "SSPStep_RK_10_4",(void(*)(void))SSPStep_RK_10_4);CHKERRQ(ierr); 319 } 320 321 ts->ops->setup = TSSetUp_SSP; 322 ts->ops->step = TSStep_SSP; 323 ts->ops->reset = TSReset_SSP; 324 ts->ops->destroy = TSDestroy_SSP; 325 ts->ops->setfromoptions = TSSetFromOptions_SSP; 326 ts->ops->view = TSView_SSP; 327 328 ierr = PetscNewLog(ts,TS_SSP,&ssp);CHKERRQ(ierr); 329 ts->data = (void*)ssp; 330 331 ierr = TSSSPSetType(ts,TSSSPRKS2);CHKERRQ(ierr); 332 ssp->nstages = 5; 333 PetscFunctionReturn(0); 334 } 335 EXTERN_C_END 336