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 xdot; 18 Vec *work; 19 PetscInt nwork; 20 PetscBool workout; 21 } TS_SSP; 22 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "SSPGetWorkVectors" 26 static PetscErrorCode SSPGetWorkVectors(TS ts,PetscInt n,Vec **work) 27 { 28 TS_SSP *ssp = (TS_SSP*)ts->data; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten"); 33 if (ssp->nwork < n) { 34 if (ssp->nwork > 0) { 35 ierr = VecDestroyVecs(ssp->nwork,&ssp->work);CHKERRQ(ierr); 36 } 37 ierr = VecDuplicateVecs(ts->vec_sol,n,&ssp->work);CHKERRQ(ierr); 38 ssp->nwork = n; 39 } 40 *work = ssp->work; 41 ssp->workout = PETSC_TRUE; 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "SSPRestoreWorkVectors" 47 static PetscErrorCode SSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work) 48 { 49 TS_SSP *ssp = (TS_SSP*)ts->data; 50 51 PetscFunctionBegin; 52 if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten"); 53 if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out"); 54 ssp->workout = PETSC_FALSE; 55 *work = PETSC_NULL; 56 PetscFunctionReturn(0); 57 } 58 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "SSPStep_RK_2" 62 /* Optimal second order SSP Runge-Kutta, low-storage, c_eff=(s-1)/s */ 63 /* Pseudocode 2 of Ketcheson 2008 */ 64 static PetscErrorCode SSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol) 65 { 66 TS_SSP *ssp = (TS_SSP*)ts->data; 67 Vec *work,F; 68 PetscInt i,s; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 s = ssp->nstages; 73 ierr = SSPGetWorkVectors(ts,2,&work);CHKERRQ(ierr); 74 F = work[1]; 75 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 76 for (i=0; i<s-1; i++) { 77 ierr = TSComputeRHSFunction(ts,t0+dt*(i/(s-1.)),work[0],F);CHKERRQ(ierr); 78 ierr = VecAXPY(work[0],dt/(s-1.),F);CHKERRQ(ierr); 79 } 80 ierr = TSComputeRHSFunction(ts,t0+dt,work[0],F);CHKERRQ(ierr); 81 ierr = VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);CHKERRQ(ierr); 82 ierr = SSPRestoreWorkVectors(ts,2,&work);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "SSPStep_RK_3" 88 /* Optimal third order SSP Runge-Kutta, low-storage, c_eff=(sqrt(s)-1)/sqrt(s), where sqrt(s) is an integer */ 89 /* Pseudocode 2 of Ketcheson 2008 */ 90 static PetscErrorCode SSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol) 91 { 92 TS_SSP *ssp = (TS_SSP*)ts->data; 93 Vec *work,F; 94 PetscInt i,s,n,r; 95 PetscReal c; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 s = ssp->nstages; 100 n = (PetscInt)(sqrt((PetscReal)s)+0.001); 101 r = s-n; 102 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); 103 ierr = SSPGetWorkVectors(ts,3,&work);CHKERRQ(ierr); 104 F = work[2]; 105 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 106 for (i=0; i<(n-1)*(n-2)/2; i++) { 107 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 108 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 109 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 110 } 111 ierr = VecCopy(work[0],work[1]);CHKERRQ(ierr); 112 for ( ; i<n*(n+1)/2-1; i++) { 113 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 114 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 115 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 116 } 117 { 118 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 119 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 120 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); 121 i++; 122 } 123 for ( ; i<s; i++) { 124 c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n); 125 ierr = TSComputeRHSFunction(ts,t0+c*dt,work[0],F);CHKERRQ(ierr); 126 ierr = VecAXPY(work[0],dt/r,F);CHKERRQ(ierr); 127 } 128 ierr = VecCopy(work[0],sol);CHKERRQ(ierr); 129 ierr = SSPRestoreWorkVectors(ts,3,&work);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "SSPStep_RK_10_4" 135 /* Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6 */ 136 /* SSPRK(10,4), Pseudocode 3 of Ketcheson 2008 */ 137 static PetscErrorCode SSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol) 138 { 139 const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1}; 140 Vec *work,F; 141 PetscInt i; 142 PetscErrorCode ierr; 143 144 PetscFunctionBegin; 145 ierr = SSPGetWorkVectors(ts,3,&work);CHKERRQ(ierr); 146 F = work[2]; 147 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 148 for (i=0; i<5; i++) { 149 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 150 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 151 } 152 ierr = VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);CHKERRQ(ierr); 153 ierr = VecAXPBY(work[0],15,-5,work[1]);CHKERRQ(ierr); 154 for ( ; i<9; i++) { 155 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 156 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 157 } 158 ierr = TSComputeRHSFunction(ts,t0+dt,work[0],F);CHKERRQ(ierr); 159 ierr = VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);CHKERRQ(ierr); 160 ierr = VecCopy(work[1],sol);CHKERRQ(ierr); 161 ierr = SSPRestoreWorkVectors(ts,3,&work);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "TSSetUp_SSP" 168 static PetscErrorCode TSSetUp_SSP(TS ts) 169 { 170 /* TS_SSP *ssp = (TS_SSP*)ts->data; */ 171 /* PetscErrorCode ierr; */ 172 173 PetscFunctionBegin; 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "TSStep_SSP" 179 static PetscErrorCode TSStep_SSP(TS ts,PetscInt *steps,PetscReal *ptime) 180 { 181 TS_SSP *ssp = (TS_SSP*)ts->data; 182 Vec sol = ts->vec_sol; 183 PetscErrorCode ierr; 184 PetscInt i,max_steps = ts->max_steps; 185 186 PetscFunctionBegin; 187 *steps = -ts->steps; 188 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 189 190 for (i=0; i<max_steps; i++) { 191 PetscReal dt = ts->time_step; 192 193 ierr = TSPreStep(ts);CHKERRQ(ierr); 194 ts->ptime += dt; 195 ierr = (*ssp->onestep)(ts,ts->ptime-dt,dt,sol);CHKERRQ(ierr); 196 ts->steps++; 197 ierr = TSPostStep(ts);CHKERRQ(ierr); 198 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 199 if (ts->ptime > ts->max_time) break; 200 } 201 202 *steps += ts->steps; 203 *ptime = ts->ptime; 204 PetscFunctionReturn(0); 205 } 206 /*------------------------------------------------------------*/ 207 #undef __FUNCT__ 208 #define __FUNCT__ "TSDestroy_SSP" 209 static PetscErrorCode TSDestroy_SSP(TS ts) 210 { 211 TS_SSP *ssp = (TS_SSP*)ts->data; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 if (ssp->work) {ierr = VecDestroyVecs(ssp->nwork,&ssp->work);CHKERRQ(ierr);} 216 ierr = PetscFree(ts->data);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 /*------------------------------------------------------------*/ 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "TSSSPSetType" 223 static PetscErrorCode TSSSPSetType(TS ts,const TSSSPType type) 224 { 225 PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec); 226 TS_SSP *ssp = (TS_SSP*)ts->data; 227 228 PetscFunctionBegin; 229 ierr = PetscFListFind(TSSSPList,((PetscObject)ts)->comm,type,(void(**)(void))&r);CHKERRQ(ierr); 230 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type); 231 ssp->onestep = r; 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "TSSetFromOptions_SSP" 237 static PetscErrorCode TSSetFromOptions_SSP(TS ts) 238 { 239 char tname[256] = TSSSPRKS2; 240 TS_SSP *ssp = (TS_SSP*)ts->data; 241 PetscErrorCode ierr; 242 PetscBool flg; 243 244 PetscFunctionBegin; 245 ierr = PetscOptionsHead("SSP ODE solver options");CHKERRQ(ierr); 246 { 247 ierr = PetscOptionsList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);CHKERRQ(ierr); 248 if (flg) { 249 ierr = TSSSPSetType(ts,tname);CHKERRQ(ierr); 250 } 251 ierr = PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,PETSC_NULL);CHKERRQ(ierr); 252 } 253 ierr = PetscOptionsTail();CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "TSView_SSP" 259 static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer) 260 { 261 PetscFunctionBegin; 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 .seealso: TSCreate(), TS, TSSetType() 303 304 M*/ 305 EXTERN_C_BEGIN 306 #undef __FUNCT__ 307 #define __FUNCT__ "TSCreate_SSP" 308 PetscErrorCode TSCreate_SSP(TS ts) 309 { 310 TS_SSP *ssp; 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 if (!TSSSPList) { 315 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS2, "SSPStep_RK_2", (void(*)(void))SSPStep_RK_2);CHKERRQ(ierr); 316 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS3, "SSPStep_RK_3", (void(*)(void))SSPStep_RK_3);CHKERRQ(ierr); 317 ierr = PetscFListAdd(&TSSSPList,TSSSPRK104, "SSPStep_RK_10_4",(void(*)(void))SSPStep_RK_10_4);CHKERRQ(ierr); 318 } 319 320 ts->ops->setup = TSSetUp_SSP; 321 ts->ops->step = TSStep_SSP; 322 ts->ops->destroy = TSDestroy_SSP; 323 ts->ops->setfromoptions = TSSetFromOptions_SSP; 324 ts->ops->view = TSView_SSP; 325 326 ierr = PetscNewLog(ts,TS_SSP,&ssp);CHKERRQ(ierr); 327 ts->data = (void*)ssp; 328 329 ierr = TSSSPSetType(ts,TSSSPRKS2);CHKERRQ(ierr); 330 ssp->nstages = 5; 331 PetscFunctionReturn(0); 332 } 333 EXTERN_C_END 334 335 336 337 338