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