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 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 249 } 250 ierr = PetscOptionsTail();CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNCT__ 255 #define __FUNCT__ "TSView_SSP" 256 static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer) 257 { 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 /* ------------------------------------------------------------ */ 266 267 /*MC 268 TSSSP - Explicit strong stability preserving ODE solver 269 270 Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation 271 bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's 272 scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties, 273 but they are usually formulated using a forward Euler time discretization or by coupling the space and time 274 discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very 275 difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the 276 semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a 277 time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP). 278 279 Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while 280 still being SSP. Some theoretical bounds 281 282 1. There are no explicit methods with c_eff > 1. 283 284 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0. 285 286 3. There are no implicit methods with order greater than 1 and c_eff > 2. 287 288 This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows 289 for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work 290 vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice. 291 292 Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104} 293 294 rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s 295 296 rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s 297 298 rk104: A 10-stage fourth order method. c_eff = 0.6 299 300 Level: beginner 301 302 References: 303 Ketcheson, Highly efficient strong stability preserving Runge-Kutta methods with low-storage implementations, SISC, 2008. 304 305 Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009. 306 307 .seealso: TSCreate(), TS, TSSetType() 308 309 M*/ 310 EXTERN_C_BEGIN 311 #undef __FUNCT__ 312 #define __FUNCT__ "TSCreate_SSP" 313 PetscErrorCode TSCreate_SSP(TS ts) 314 { 315 TS_SSP *ssp; 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 if (!TSSSPList) { 320 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS2, "SSPStep_RK_2", (void(*)(void))SSPStep_RK_2);CHKERRQ(ierr); 321 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS3, "SSPStep_RK_3", (void(*)(void))SSPStep_RK_3);CHKERRQ(ierr); 322 ierr = PetscFListAdd(&TSSSPList,TSSSPRK104, "SSPStep_RK_10_4",(void(*)(void))SSPStep_RK_10_4);CHKERRQ(ierr); 323 } 324 325 ts->ops->setup = TSSetUp_SSP; 326 ts->ops->step = TSStep_SSP; 327 ts->ops->reset = TSReset_SSP; 328 ts->ops->destroy = TSDestroy_SSP; 329 ts->ops->setfromoptions = TSSetFromOptions_SSP; 330 ts->ops->view = TSView_SSP; 331 332 ierr = PetscNewLog(ts,TS_SSP,&ssp);CHKERRQ(ierr); 333 ts->data = (void*)ssp; 334 335 ierr = TSSSPSetType(ts,TSSSPRKS2);CHKERRQ(ierr); 336 ssp->nstages = 5; 337 PetscFunctionReturn(0); 338 } 339 EXTERN_C_END 340