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_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_ERR_ORDER,"Work vectors have not been gotten"); 54 if (*work != ssp->work) SETERRQ(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 = floor(sqrt(s)); 102 r = s-n; 103 if (n*n != s) SETERRQ1(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 TS_SSP *ssp = (TS_SSP*)ts->data; 141 const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1}; 142 Vec *work,F; 143 PetscInt i,s; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 s = ssp->nstages; 148 ierr = SSPGetWorkVectors(ts,3,&work);CHKERRQ(ierr); 149 F = work[2]; 150 ierr = VecCopy(sol,work[0]);CHKERRQ(ierr); 151 for (i=0; i<5; i++) { 152 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 153 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 154 } 155 ierr = VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);CHKERRQ(ierr); 156 ierr = VecAXPBY(work[0],15,-5,work[1]);CHKERRQ(ierr); 157 for ( ; i<9; i++) { 158 ierr = TSComputeRHSFunction(ts,t0+c[i],work[0],F);CHKERRQ(ierr); 159 ierr = VecAXPY(work[0],dt/6,F);CHKERRQ(ierr); 160 } 161 ierr = TSComputeRHSFunction(ts,t0+dt,work[0],F);CHKERRQ(ierr); 162 ierr = VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);CHKERRQ(ierr); 163 ierr = VecCopy(work[1],sol);CHKERRQ(ierr); 164 ierr = SSPRestoreWorkVectors(ts,3,&work);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "TSSetUp_SSP" 171 static PetscErrorCode TSSetUp_SSP(TS ts) 172 { 173 /* TS_SSP *ssp = (TS_SSP*)ts->data; */ 174 /* PetscErrorCode ierr; */ 175 176 PetscFunctionBegin; 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "TSStep_SSP" 182 static PetscErrorCode TSStep_SSP(TS ts,PetscInt *steps,PetscReal *ptime) 183 { 184 TS_SSP *ssp = (TS_SSP*)ts->data; 185 Vec sol = ts->vec_sol; 186 PetscErrorCode ierr; 187 PetscInt i,max_steps = ts->max_steps; 188 189 PetscFunctionBegin; 190 *steps = -ts->steps; 191 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 192 193 for (i=0; i<max_steps; i++) { 194 PetscReal dt = ts->time_step; 195 196 ts->ptime += dt; 197 ierr = (*ssp->onestep)(ts,ts->ptime-dt,dt,sol);CHKERRQ(ierr); 198 ts->steps++; 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_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 - ODE solver using the explicit SSP method 270 271 Level: beginner 272 273 .seealso: TSCreate(), TS, TSSetType() 274 275 M*/ 276 EXTERN_C_BEGIN 277 #undef __FUNCT__ 278 #define __FUNCT__ "TSCreate_SSP" 279 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_SSP(TS ts) 280 { 281 TS_SSP *ssp; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 if (!TSSSPList) { 286 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS2, "SSPStep_RK_2", (void(*)(void))SSPStep_RK_2);CHKERRQ(ierr); 287 ierr = PetscFListAdd(&TSSSPList,TSSSPRKS3, "SSPStep_RK_3", (void(*)(void))SSPStep_RK_3);CHKERRQ(ierr); 288 ierr = PetscFListAdd(&TSSSPList,TSSSPRK104, "SSPStep_RK_10_4",(void(*)(void))SSPStep_RK_10_4);CHKERRQ(ierr); 289 } 290 291 ts->ops->setup = TSSetUp_SSP; 292 ts->ops->step = TSStep_SSP; 293 ts->ops->destroy = TSDestroy_SSP; 294 ts->ops->setfromoptions = TSSetFromOptions_SSP; 295 ts->ops->view = TSView_SSP; 296 297 ierr = PetscNewLog(ts,TS_SSP,&ssp);CHKERRQ(ierr); 298 ts->data = (void*)ssp; 299 300 ierr = TSSSPSetType(ts,TSSSPRKS2);CHKERRQ(ierr); 301 ssp->nstages = 5; 302 PetscFunctionReturn(0); 303 } 304 EXTERN_C_END 305 306 307 308 309