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