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,PetscInt *steps,PetscReal *ptime) 178 { 179 TS_SSP *ssp = (TS_SSP*)ts->data; 180 Vec sol = ts->vec_sol; 181 PetscInt i; 182 PetscErrorCode ierr; 183 184 PetscFunctionBegin; 185 *steps = -ts->steps; 186 *ptime = ts->ptime; 187 188 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 189 190 for (i=0; i<ts->max_steps; i++) { 191 if (ts->ptime + ts->time_step > ts->max_time) break; 192 ierr = TSPreStep(ts);CHKERRQ(ierr); 193 194 ierr = (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);CHKERRQ(ierr); 195 196 ts->ptime += ts->time_step; 197 ts->steps++; 198 199 ierr = TSPostStep(ts);CHKERRQ(ierr); 200 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 201 } 202 203 *steps += ts->steps; 204 *ptime = ts->ptime; 205 PetscFunctionReturn(0); 206 } 207 /*------------------------------------------------------------*/ 208 #undef __FUNCT__ 209 #define __FUNCT__ "TSReset_SSP" 210 static PetscErrorCode TSReset_SSP(TS ts) 211 { 212 TS_SSP *ssp = (TS_SSP*)ts->data; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 if (ssp->work) {ierr = VecDestroyVecs(ssp->nwork,&ssp->work);CHKERRQ(ierr);} 217 ssp->nwork = 0; 218 ssp->workout = PETSC_FALSE; 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "TSDestroy_SSP" 224 static PetscErrorCode TSDestroy_SSP(TS ts) 225 { 226 TS_SSP *ssp = (TS_SSP*)ts->data; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 ierr = TSReset_SSP(ts);CHKERRQ(ierr); 231 ierr = PetscFree(ts->data);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 /*------------------------------------------------------------*/ 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "TSSSPSetType" 238 static PetscErrorCode TSSSPSetType(TS ts,const TSSSPType type) 239 { 240 PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec); 241 TS_SSP *ssp = (TS_SSP*)ts->data; 242 243 PetscFunctionBegin; 244 ierr = PetscFListFind(TSSSPList,((PetscObject)ts)->comm,type,PETSC_TRUE,(void(**)(void))&r);CHKERRQ(ierr); 245 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type); 246 ssp->onestep = r; 247 PetscFunctionReturn(0); 248 } 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "TSSetFromOptions_SSP" 252 static PetscErrorCode TSSetFromOptions_SSP(TS ts) 253 { 254 char tname[256] = TSSSPRKS2; 255 TS_SSP *ssp = (TS_SSP*)ts->data; 256 PetscErrorCode ierr; 257 PetscBool flg; 258 259 PetscFunctionBegin; 260 ierr = PetscOptionsHead("SSP ODE solver options");CHKERRQ(ierr); 261 { 262 ierr = PetscOptionsList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);CHKERRQ(ierr); 263 if (flg) { 264 ierr = TSSSPSetType(ts,tname);CHKERRQ(ierr); 265 } 266 ierr = PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,PETSC_NULL);CHKERRQ(ierr); 267 } 268 ierr = PetscOptionsTail();CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "TSView_SSP" 274 static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer) 275 { 276 PetscFunctionBegin; 277 PetscFunctionReturn(0); 278 } 279 280 /* ------------------------------------------------------------ */ 281 282 /*MC 283 TSSSP - Explicit strong stability preserving ODE solver 284 285 Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation 286 bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's 287 scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties, 288 but they are usually formulated using a forward Euler time discretization or by coupling the space and time 289 discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very 290 difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the 291 semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a 292 time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP). 293 294 Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while 295 still being SSP. Some theoretical bounds 296 297 1. There are no explicit methods with c_eff > 1. 298 299 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0. 300 301 3. There are no implicit methods with order greater than 1 and c_eff > 2. 302 303 This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows 304 for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work 305 vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice. 306 307 Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104} 308 309 rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s 310 311 rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s 312 313 rk104: A 10-stage fourth order method. c_eff = 0.6 314 315 Level: beginner 316 317 .seealso: TSCreate(), TS, TSSetType() 318 319 M*/ 320 EXTERN_C_BEGIN 321 #undef __FUNCT__ 322 #define __FUNCT__ "TSCreate_SSP" 323 PetscErrorCode TSCreate_SSP(TS ts) 324 { 325 TS_SSP *ssp; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 if (!TSSSPList) { 330 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS2, "SSPStep_RK_2", (void(*)(void))SSPStep_RK_2);CHKERRQ(ierr); 331 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS3, "SSPStep_RK_3", (void(*)(void))SSPStep_RK_3);CHKERRQ(ierr); 332 ierr = PetscFListAdd(&TSSSPList,TSSSPRK104, "SSPStep_RK_10_4",(void(*)(void))SSPStep_RK_10_4);CHKERRQ(ierr); 333 } 334 335 ts->ops->setup = TSSetUp_SSP; 336 ts->ops->step = TSStep_SSP; 337 ts->ops->reset = TSReset_SSP; 338 ts->ops->destroy = TSDestroy_SSP; 339 ts->ops->setfromoptions = TSSetFromOptions_SSP; 340 ts->ops->view = TSView_SSP; 341 342 ierr = PetscNewLog(ts,TS_SSP,&ssp);CHKERRQ(ierr); 343 ts->data = (void*)ssp; 344 345 ierr = TSSSPSetType(ts,TSSSPRKS2);CHKERRQ(ierr); 346 ssp->nstages = 5; 347 PetscFunctionReturn(0); 348 } 349 EXTERN_C_END 350