1 /* 2 Code for time stepping with the Partitioned Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) for nonsplit RHS multi-rate RK, 7 user should give the indexes for both slow and fast components; 8 2) The general system is written as 9 Usdot = Fs(t,Us,Uf) 10 Ufdot = Ff(t,Us,Uf) for split RHS multi-rate RK, 11 user should partioned RHS by themselves and also provide the indexes for both slow and fast components. 12 3) To correct The confusing terminology in the paper, we use 'slow method', 'slow buffer method' and 'fast method' to denote the methods applied to 'slow region', 'slow buffer region' and 'fast region' respectively. The 'slow method' in the original paper actually means the 'slow buffer method'. 13 4) Why does the buffer region have to be inside the slow region? The buffer region is treated with a slow method essentially. Applying the slow method to a region with a fast characteristic time scale is apparently not a good choice. 14 15 Reference: 16 Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007 17 18 */ 19 20 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 21 #include <petscdm.h> 22 23 static TSMPRKType TSMPRKDefault = TSMPRK2A22; 24 static PetscBool TSMPRKRegisterAllCalled; 25 static PetscBool TSMPRKPackageInitialized; 26 27 typedef struct _MPRKTableau *MPRKTableau; 28 struct _MPRKTableau { 29 char *name; 30 PetscInt order; /* Classical approximation order of the method i */ 31 PetscInt s; /* Number of stages */ 32 PetscInt np; /* Number of partitions */ 33 PetscReal *Af,*bf,*cf; /* Tableau for fast components */ 34 PetscReal *Amb,*bmb,*cmb; /* Tableau for medium components */ 35 PetscInt *rmb; /* Array of flags for repeated stages in medium method */ 36 PetscReal *Asb,*bsb,*csb; /* Tableau for slow components */ 37 PetscInt *rsb; /* Array of flags for repeated staged in slow method*/ 38 }; 39 typedef struct _MPRKTableauLink *MPRKTableauLink; 40 struct _MPRKTableauLink { 41 struct _MPRKTableau tab; 42 MPRKTableauLink next; 43 }; 44 static MPRKTableauLink MPRKTableauList; 45 46 typedef struct { 47 MPRKTableau tableau; 48 TSMPRKMultirateType mprkmtype; 49 Vec *Y; /* States computed during the step */ 50 Vec Ytmp; 51 Vec *YdotRHS_slow; /* Function evaluations by slow tableau for slow components */ 52 Vec *YdotRHS_slowbuffer; /* Function evaluations by slow tableau for slow components */ 53 Vec *YdotRHS_medium; /* Function evaluations by slow tableau for slow components */ 54 Vec *YdotRHS_mediumbuffer; /* Function evaluations by slow tableau for slow components */ 55 Vec *YdotRHS_fast; /* Function evaluations by fast tableau for fast components */ 56 PetscScalar *work_slow; /* Scalar work_slow by slow tableau */ 57 PetscScalar *work_slowbuffer; /* Scalar work_slow by slow tableau */ 58 PetscScalar *work_medium; /* Scalar work_slow by medium tableau */ 59 PetscScalar *work_mediumbuffer; /* Scalar work_slow by medium tableau */ 60 PetscScalar *work_fast; /* Scalar work_fast by fast tableau */ 61 PetscReal stage_time; 62 TSStepStatus status; 63 PetscReal ptime; 64 PetscReal time_step; 65 IS is_slow,is_slowbuffer,is_medium,is_mediumbuffer,is_fast; 66 TS subts_slow,subts_slowbuffer,subts_medium,subts_mediumbuffer,subts_fast; 67 } TS_MPRK; 68 69 static PetscErrorCode TSMPRKGenerateTableau2(PetscInt ratio,PetscInt s,const PetscReal Abase[],const PetscReal bbase[],PetscReal A1[],PetscReal b1[],PetscReal A2[],PetscReal b2[]) 70 { 71 PetscInt i,j,k,l; 72 73 PetscFunctionBegin; 74 for (k=0; k<ratio; k++) { 75 /* diagonal blocks */ 76 for (i=0; i<s; i++) 77 for (j=0; j<s; j++) { 78 A1[(k*s+i)*ratio*s+k*s+j] = Abase[i*s+j]; 79 A2[(k*s+i)*ratio*s+k*s+j] = Abase[i*s+j]/ratio; 80 } 81 /* off diagonal blocks */ 82 for (l=0; l<k; l++) 83 for (i=0; i<s; i++) 84 for (j=0; j<s; j++) 85 A2[(k*s+i)*ratio*s+l*s+j] = bbase[j]/ratio; 86 for (j=0; j<s; j++) { 87 b1[k*s+j] = bbase[j]/ratio; 88 b2[k*s+j] = bbase[j]/ratio; 89 } 90 } 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode TSMPRKGenerateTableau3(PetscInt ratio,PetscInt s,const PetscReal Abase[],const PetscReal bbase[],PetscReal A1[],PetscReal b1[],PetscReal A2[],PetscReal b2[],PetscReal A3[],PetscReal b3[]) 95 { 96 PetscInt i,j,k,l,m,n; 97 98 PetscFunctionBegin; 99 for (k=0; k<ratio; k++) { /* diagonal blocks of size ratio*s by ratio*s */ 100 for (l=0; l<ratio; l++) /* diagonal sub-blocks of size s by s */ 101 for (i=0; i<s; i++) 102 for (j=0; j<s; j++) { 103 A1[((k*ratio+l)*s+i)*ratio*ratio*s+(k*ratio+l)*s+j] = Abase[i*s+j]; 104 A2[((k*ratio+l)*s+i)*ratio*ratio*s+(k*ratio+l)*s+j] = Abase[i*s+j]/ratio; 105 A3[((k*ratio+l)*s+i)*ratio*ratio*s+(k*ratio+l)*s+j] = Abase[i*s+j]/ratio/ratio; 106 } 107 for (l=0; l<k; l++) /* off-diagonal blocks of size ratio*s by ratio*s */ 108 for (m=0; m<ratio; m++) 109 for (n=0; n<ratio; n++) 110 for (i=0; i<s; i++) 111 for (j=0; j<s; j++) { 112 A2[((k*ratio+m)*s+i)*ratio*ratio*s+(l*ratio+n)*s+j] = bbase[j]/ratio/ratio; 113 A3[((k*ratio+m)*s+i)*ratio*ratio*s+(l*ratio+n)*s+j] = bbase[j]/ratio/ratio; 114 } 115 for (m=0; m<ratio; m++) 116 for (n=0; n<m; n++) /* off-diagonal sub-blocks of size s by s in the diagonal blocks */ 117 for (i=0; i<s; i++) 118 for (j=0; j<s; j++) 119 A3[((k*ratio+m)*s+i)*ratio*ratio*s+(k*ratio+n)*s+j] = bbase[j]/ratio/ratio; 120 for (n=0; n<ratio; n++) 121 for (j=0; j<s; j++) { 122 b1[(k*ratio+n)*s+j] = bbase[j]/ratio/ratio; 123 b2[(k*ratio+n)*s+j] = bbase[j]/ratio/ratio; 124 b3[(k*ratio+n)*s+j] = bbase[j]/ratio/ratio; 125 } 126 } 127 PetscFunctionReturn(0); 128 } 129 130 /*MC 131 TSMPRK2A22 - Second Order Partitioned Runge Kutta scheme based on RK2A. 132 133 This method has four stages for slow and fast parts. The refinement factor of the stepsize is 2. 134 r = 2, np = 2 135 Options database: 136 . -ts_mprk_type 2a22 137 138 Level: advanced 139 140 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 141 M*/ 142 /*MC 143 TSMPRK2A23 - Second Order Partitioned Runge Kutta scheme based on RK2A. 144 145 This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 2. 146 r = 2, np = 3 147 Options database: 148 . -ts_mprk_type 2a23 149 150 Level: advanced 151 152 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 153 M*/ 154 /*MC 155 TSMPRK2A32 - Second Order Partitioned Runge Kutta scheme based on RK2A. 156 157 This method has four stages for slow and fast parts. The refinement factor of the stepsize is 3. 158 r = 3, np = 2 159 Options database: 160 . -ts_mprk_type 2a32 161 162 Level: advanced 163 164 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 165 M*/ 166 /*MC 167 TSMPRK2A33 - Second Order Partitioned Runge Kutta scheme based on RK2A. 168 169 This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 3. 170 r = 3, np = 3 171 Options database: 172 . -ts_mprk_type 2a33 173 174 Level: advanced 175 176 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 177 M*/ 178 /*MC 179 TSMPRK3P2M - Third Order Partitioned Runge Kutta scheme. 180 181 This method has eight stages for both slow and fast parts. 182 183 Options database: 184 . -ts_mprk_type pm3 (put here temporarily) 185 186 Level: advanced 187 188 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 189 M*/ 190 /*MC 191 TSMPRKP2 - Second Order Partitioned Runge Kutta scheme. 192 193 This method has five stages for both slow and fast parts. 194 195 Options database: 196 . -ts_mprk_type p2 197 198 Level: advanced 199 200 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 201 M*/ 202 /*MC 203 TSMPRKP3 - Third Order Partitioned Runge Kutta scheme. 204 205 This method has ten stages for both slow and fast parts. 206 207 Options database: 208 . -ts_mprk_type p3 209 210 Level: advanced 211 212 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 213 M*/ 214 215 /*@C 216 TSMPRKRegisterAll - Registers all of the Partirioned Runge-Kutta explicit methods in TSMPRK 217 218 Not Collective, but should be called by all processes which will need the schemes to be registered 219 220 Level: advanced 221 222 .keywords: TS, TSMPRK, register, all 223 224 .seealso: TSMPRKRegisterDestroy() 225 @*/ 226 PetscErrorCode TSMPRKRegisterAll(void) 227 { 228 PetscErrorCode ierr; 229 230 PetscFunctionBegin; 231 if (TSMPRKRegisterAllCalled) PetscFunctionReturn(0); 232 TSMPRKRegisterAllCalled = PETSC_TRUE; 233 234 #define RC PetscRealConstant 235 { 236 const PetscReal 237 Abase[2][2] = {{0,0}, 238 {RC(1.0),0}}, 239 bbase[2] = {RC(0.5),RC(0.5)}; 240 PetscReal 241 Asb[4][4] = {{0}},Af[4][4] = {{0}},bsb[4] = {0},bf[4] = {0}; 242 PetscInt 243 rsb[4] = {0,0,1,2}; 244 ierr = TSMPRKGenerateTableau2(2,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Af[0][0],bf);CHKERRQ(ierr); 245 ierr = TSMPRKRegister(TSMPRK2A22,2,4,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL);CHKERRQ(ierr); 246 } 247 { 248 const PetscReal 249 Abase[2][2] = {{0,0}, 250 {RC(1.0),0}}, 251 bbase[2] = {RC(0.5),RC(0.5)}; 252 PetscReal 253 Asb[8][8] = {{0}},Amb[8][8] = {{0}},Af[8][8] = {{0}},bsb[8] ={0},bmb[8] = {0},bf[8] = {0}; 254 PetscInt 255 rsb[8] = {0,0,1,2,1,2,1,2},rmb[8] = {0,0,1,2,0,0,5,6}; 256 ierr = TSMPRKGenerateTableau3(2,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Amb[0][0],bmb,&Af[0][0],bf);CHKERRQ(ierr); 257 ierr = TSMPRKRegister(TSMPRK2A23,2,8,&Asb[0][0],bsb,NULL,rsb,&Amb[0][0],bmb,NULL,rmb,&Af[0][0],bf,NULL);CHKERRQ(ierr); 258 } 259 { 260 const PetscReal 261 Abase[2][2] = {{0,0}, 262 {RC(1.0),0}}, 263 bbase[2] = {RC(0.5),RC(0.5)}; 264 PetscReal 265 Asb[6][6] = {{0}},Af[6][6] = {{0}},bsb[6] = {0},bf[6] = {0}; 266 PetscInt 267 rsb[6] = {0,0,1,2,1,2}; 268 ierr = TSMPRKGenerateTableau2(3,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Af[0][0],bf);CHKERRQ(ierr); 269 ierr = TSMPRKRegister(TSMPRK2A32,2,6,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL);CHKERRQ(ierr); 270 } 271 { 272 const PetscReal 273 Abase[2][2] = {{0,0}, 274 {RC(1.0),0}}, 275 bbase[2] = {RC(0.5),RC(0.5)}; 276 PetscReal 277 Asb[18][18] = {{0}},Amb[18][18] = {{0}},Af[18][18] = {{0}},bsb[18] ={0},bmb[18] = {0},bf[18] = {0}; 278 PetscInt 279 rsb[18] = {0,0,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2},rmb[18] = {0,0,1,2,1,2,0,0,7,8,7,8,0,0,13,14,13,14}; 280 ierr = TSMPRKGenerateTableau3(3,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Amb[0][0],bmb,&Af[0][0],bf);CHKERRQ(ierr); 281 ierr = TSMPRKRegister(TSMPRK2A33,2,18,&Asb[0][0],bsb,NULL,rsb,&Amb[0][0],bmb,NULL,rmb,&Af[0][0],bf,NULL);CHKERRQ(ierr); 282 } 283 /* 284 PetscReal 285 Asb[8][8] = {{Abase[0][0],Abase[0][1],0,0,0,0,0,0}, 286 {Abase[1][0],Abase[1][1],0,0,0,0,0,0}, 287 {0,0,Abase[0][0],Abase[0][1],0,0,0,0}, 288 {0,0,Abase[1][0],Abase[1][1],0,0,0,0}, 289 {0,0,0,0,Abase[0][0],Abase[0][1],0,0}, 290 {0,0,0,0,Abase[1][0],Abase[1][1],0,0}, 291 {0,0,0,0,0,0,Abase[0][0],Abase[0][1]}, 292 {0,0,0,0,0,0,Abase[1][0],Abase[1][1]}}, 293 Amb[8][8] = {{Abase[0][0]/m,Abase[0][1]/m,0,0,0,0,0,0}, 294 {Abase[1][0]/m,Abase[1][1]/m,0,0,0,0,0,0}, 295 {0,0,Abase[0][0]/m,Abase[0][1]/m,0,0,0,0}, 296 {0,0,Abase[1][0]/m,Abase[1][1]/m,0,0,0,0}, 297 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[0][0]/m,Abase[0][1]/m,0,0}, 298 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[1][0]/m,Abase[1][1]/m,0,0}, 299 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[0][0]/m,Abase[0][1]/m}, 300 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[1][0]/m,Abase[1][1]/m}}, 301 Af[8][8] = {{Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0,0,0}, 302 {Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0,0,0}, 303 {0,0,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0}, 304 {0,0,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0}, 305 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0}, 306 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0}, 307 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m,Abase[0][1]/m}, 308 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m,Abase[1][1]/m}}, 309 bsb[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m}, 310 bmb[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m,bbase[1]/m/m}, 311 bf[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m,bbase[0]/m/m,bbase[1]/m/m}, 312 */ 313 /*{ 314 const PetscReal 315 As[8][8] = {{0,0,0,0,0,0,0,0}, 316 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0}, 317 {RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0,0,0,0,0}, 318 {RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0,0}, 319 {0,0,0,0,0,0,0,0}, 320 {0,0,0,0,RC(1.0)/RC(2.0),0,0,0}, 321 {0,0,0,0,RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0}, 322 {0,0,0,0,RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0}}, 323 A[8][8] = {{0,0,0,0,0,0,0,0}, 324 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0}, 325 {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0}, 326 {RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0,0}, 327 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,0,0,0}, 328 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(4.0),0,0,0}, 329 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0}, 330 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0}}, 331 bs[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)}, 332 b[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)}; 333 ierr = TSMPRKRegister(TSMPRKPM3,3,8,&As[0][0],bs,NULL,&A[0][0],b,NULL);CHKERRQ(ierr); 334 }*/ 335 336 { 337 const PetscReal 338 Asb[5][5] = {{0,0,0,0,0}, 339 {RC(1.0)/RC(2.0),0,0,0,0}, 340 {RC(1.0)/RC(2.0),0,0,0,0}, 341 {RC(1.0),0,0,0,0}, 342 {RC(1.0),0,0,0,0}}, 343 Af[5][5] = {{0,0,0,0,0}, 344 {RC(1.0)/RC(2.0),0,0,0,0}, 345 {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),0,0,0}, 346 {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(2.0),0,0}, 347 {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),0}}, 348 bsb[5] = {RC(1.0)/RC(2.0),0,0,0,RC(1.0)/RC(2.0)}, 349 bf[5] = {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),0}; 350 const PetscInt 351 rsb[5] = {0,0,2,0,4}; 352 ierr = TSMPRKRegister(TSMPRKP2,2,5,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL);CHKERRQ(ierr); 353 } 354 355 { 356 const PetscReal 357 Asb[10][10] = {{0,0,0,0,0,0,0,0,0,0}, 358 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 359 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 360 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0,0,0}, 361 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0,0,0}, 362 {RC(-1.0)/RC(6.0),0,0,0,RC(2.0)/RC(3.0),0,0,0,0,0}, 363 {RC(1.0)/RC(12.0),0,0,0,RC(1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0}, 364 {RC(1.0)/RC(12.0),0,0,0,RC(1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0}, 365 {RC(1.0)/RC(3.0),0,0,0,RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0}, 366 {RC(1.0)/RC(3.0),0,0,0,RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0}}, 367 Af[10][10] = {{0,0,0,0,0,0,0,0,0,0}, 368 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 369 {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0,0,0}, 370 {RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0,0,0,0}, 371 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,0,0,0,0,0}, 372 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,0,0,0,0,0}, 373 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,RC(1.0)/RC(4.0),0,0,0,0}, 374 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0}, 375 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0}, 376 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0}}, 377 bsb[10] = {RC(1.0)/RC(6.0),0,0,0,RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),0,0,0,RC(1.0)/RC(6.0)}, 378 bf[10] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0}; 379 const PetscInt 380 rsb[10] = {0,0,2,0,4,0,0,7,0,9}; 381 ierr = TSMPRKRegister(TSMPRKP3,3,10,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL);CHKERRQ(ierr); 382 } 383 #undef RC 384 PetscFunctionReturn(0); 385 } 386 387 /*@C 388 TSMPRKRegisterDestroy - Frees the list of schemes that were registered by TSMPRKRegister(). 389 390 Not Collective 391 392 Level: advanced 393 394 .keywords: TSMPRK, register, destroy 395 .seealso: TSMPRKRegister(), TSMPRKRegisterAll() 396 @*/ 397 PetscErrorCode TSMPRKRegisterDestroy(void) 398 { 399 PetscErrorCode ierr; 400 MPRKTableauLink link; 401 402 PetscFunctionBegin; 403 while ((link = MPRKTableauList)) { 404 MPRKTableau t = &link->tab; 405 MPRKTableauList = link->next; 406 ierr = PetscFree3(t->Asb,t->bsb,t->csb);CHKERRQ(ierr); 407 ierr = PetscFree3(t->Amb,t->bmb,t->cmb);CHKERRQ(ierr); 408 ierr = PetscFree3(t->Af,t->bf,t->cf);CHKERRQ(ierr); 409 ierr = PetscFree2(t->rsb,t->rmb);CHKERRQ(ierr); 410 ierr = PetscFree (t->name);CHKERRQ(ierr); 411 ierr = PetscFree (link);CHKERRQ(ierr); 412 } 413 TSMPRKRegisterAllCalled = PETSC_FALSE; 414 PetscFunctionReturn(0); 415 } 416 417 /*@C 418 TSMPRKInitializePackage - This function initializes everything in the TSMPRK package. It is called 419 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_MPRK() 420 when using static libraries. 421 422 Level: developer 423 424 .keywords: TS, TSMPRK, initialize, package 425 .seealso: PetscInitialize() 426 @*/ 427 PetscErrorCode TSMPRKInitializePackage(void) 428 { 429 PetscErrorCode ierr; 430 431 PetscFunctionBegin; 432 if (TSMPRKPackageInitialized) PetscFunctionReturn(0); 433 TSMPRKPackageInitialized = PETSC_TRUE; 434 ierr = TSMPRKRegisterAll();CHKERRQ(ierr); 435 ierr = PetscRegisterFinalize(TSMPRKFinalizePackage);CHKERRQ(ierr); 436 PetscFunctionReturn(0); 437 } 438 439 /*@C 440 TSRKFinalizePackage - This function destroys everything in the TSMPRK package. It is 441 called from PetscFinalize(). 442 443 Level: developer 444 445 .keywords: Petsc, destroy, package 446 .seealso: PetscFinalize() 447 @*/ 448 PetscErrorCode TSMPRKFinalizePackage(void) 449 { 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 TSMPRKPackageInitialized = PETSC_FALSE; 454 ierr = TSMPRKRegisterDestroy();CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457 458 /*@C 459 TSMPRKRegister - register a MPRK scheme by providing the entries in the Butcher tableau 460 461 Not Collective, but the same schemes should be registered on all processes on which they will be used 462 463 Input Parameters: 464 + name - identifier for method 465 . order - approximation order of method 466 . s - number of stages, this is the dimension of the matrices below 467 . Af - stage coefficients for fast components(dimension s*s, row-major) 468 . bf - step completion table for fast components(dimension s) 469 . cf - abscissa for fast components(dimension s) 470 . As - stage coefficients for slow components(dimension s*s, row-major) 471 . bs - step completion table for slow components(dimension s) 472 - cs - abscissa for slow components(dimension s) 473 474 Notes: 475 Several MPRK methods are provided, this function is only needed to create new methods. 476 477 Level: advanced 478 479 .keywords: TS, register 480 481 .seealso: TSMPRK 482 @*/ 483 PetscErrorCode TSMPRKRegister(TSMPRKType name,PetscInt order,PetscInt s, 484 const PetscReal Asb[],const PetscReal bsb[],const PetscReal csb[],const PetscInt rsb[], 485 const PetscReal Amb[],const PetscReal bmb[],const PetscReal cmb[],const PetscInt rmb[], 486 const PetscReal Af[],const PetscReal bf[],const PetscReal cf[]) 487 { 488 MPRKTableauLink link; 489 MPRKTableau t; 490 PetscInt i,j; 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 PetscValidCharPointer(name,1); 495 PetscValidRealPointer(Asb,4); 496 if (bsb) PetscValidRealPointer(bsb,5); 497 if (csb) PetscValidRealPointer(csb,6); 498 if (rsb) PetscValidRealPointer(rsb,7); 499 if (Amb) PetscValidRealPointer(Amb,8); 500 if (bmb) PetscValidRealPointer(bmb,9); 501 if (cmb) PetscValidRealPointer(cmb,10); 502 if (rmb) PetscValidRealPointer(rmb,11); 503 PetscValidRealPointer(Af,12); 504 if (bf) PetscValidRealPointer(bf,8); 505 if (cf) PetscValidRealPointer(cf,9); 506 507 ierr = PetscNew(&link);CHKERRQ(ierr); 508 t = &link->tab; 509 510 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 511 t->order = order; 512 t->s = s; 513 t->np = 2; 514 ierr = PetscMalloc3(s*s,&t->Af,s,&t->bf,s,&t->cf);CHKERRQ(ierr); 515 ierr = PetscMemcpy(t->Af,Af,s*s*sizeof(Af[0]));CHKERRQ(ierr); 516 if (bf) { 517 ierr = PetscMemcpy(t->bf,bf,s*sizeof(bf[0]));CHKERRQ(ierr); 518 } else 519 for (i=0; i<s; i++) t->bf[i] = Af[(s-1)*s+i]; 520 if (cf) { 521 ierr = PetscMemcpy(t->cf,cf,s*sizeof(cf[0]));CHKERRQ(ierr); 522 } else { 523 for (i=0; i<s; i++) 524 for (j=0,t->cf[i]=0; j<s; j++) 525 t->cf[i] += Af[i*s+j]; 526 } 527 528 if (Amb) { 529 t->np = 3; 530 ierr = PetscMalloc3(s*s,&t->Amb,s,&t->bmb,s,&t->cmb);CHKERRQ(ierr); 531 ierr = PetscCalloc1(s,&t->rmb);CHKERRQ(ierr); 532 ierr = PetscMemcpy(t->Amb,Amb,s*s*sizeof(Amb[0]));CHKERRQ(ierr); 533 if (bmb) { 534 ierr = PetscMemcpy(t->bmb,bmb,s*sizeof(bmb[0]));CHKERRQ(ierr); 535 } else { 536 for (i=0; i<s; i++) t->bmb[i] = Amb[(s-1)*s+i]; 537 } 538 if (cmb) { 539 ierr = PetscMemcpy(t->cmb,cmb,s*sizeof(cmb[0]));CHKERRQ(ierr); 540 } else { 541 for (i=0; i<s; i++) 542 for (j=0,t->cmb[i]=0; j<s; j++) 543 t->cmb[i] += Amb[i*s+j]; 544 } 545 if (rmb) { 546 ierr = PetscMemcpy(t->rmb,rmb,s*sizeof(rmb[0]));CHKERRQ(ierr); 547 } 548 } 549 550 ierr = PetscMalloc3(s*s,&t->Asb,s,&t->bsb,s,&t->csb);CHKERRQ(ierr); 551 ierr = PetscMemcpy(t->Asb,Asb,s*s*sizeof(Asb[0]));CHKERRQ(ierr); 552 ierr = PetscCalloc1(s,&t->rsb);CHKERRQ(ierr); 553 if (bsb) { 554 ierr = PetscMemcpy(t->bsb,bsb,s*sizeof(bsb[0]));CHKERRQ(ierr); 555 } else 556 for (i=0; i<s; i++) t->bsb[i] = Asb[(s-1)*s+i]; 557 if (csb) { 558 ierr = PetscMemcpy(t->csb,csb,s*sizeof(csb[0]));CHKERRQ(ierr); 559 } else { 560 for (i=0; i<s; i++) 561 for (j=0,t->csb[i]=0; j<s; j++) 562 t->csb[i] += Asb[i*s+j]; 563 } 564 if (rsb) { 565 ierr = PetscMemcpy(t->rsb,rsb,s*sizeof(rsb[0]));CHKERRQ(ierr); 566 } 567 link->next = MPRKTableauList; 568 MPRKTableauList = link; 569 PetscFunctionReturn(0); 570 } 571 572 static PetscErrorCode TSMPRKSetSplits(TS ts) 573 { 574 TS_MPRK *mprk = (TS_MPRK*)ts->data; 575 MPRKTableau tab = mprk->tableau; 576 DM dm,subdm,newdm; 577 PetscErrorCode ierr; 578 579 PetscFunctionBegin; 580 ierr = TSRHSSplitGetSubTS(ts,"slow",&mprk->subts_slow);CHKERRQ(ierr); 581 ierr = TSRHSSplitGetSubTS(ts,"fast",&mprk->subts_fast);CHKERRQ(ierr); 582 if (!mprk->subts_slow || !mprk->subts_fast) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up the RHSFunctions for 'slow' and 'fast' components using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS"); 583 584 /* Only copy the DM */ 585 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 586 587 ierr = TSRHSSplitGetSubTS(ts,"slowbuffer",&mprk->subts_slowbuffer);CHKERRQ(ierr); 588 if (!mprk->subts_slowbuffer) { 589 mprk->subts_slowbuffer = mprk->subts_slow; 590 mprk->subts_slow = NULL; 591 } 592 if (mprk->subts_slow) { 593 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 594 ierr = TSGetDM(mprk->subts_slow,&subdm);CHKERRQ(ierr); 595 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 596 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 597 ierr = TSSetDM(mprk->subts_slow,newdm);CHKERRQ(ierr); 598 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 599 } 600 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 601 ierr = TSGetDM(mprk->subts_slowbuffer,&subdm);CHKERRQ(ierr); 602 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 603 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 604 ierr = TSSetDM(mprk->subts_slowbuffer,newdm);CHKERRQ(ierr); 605 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 606 607 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 608 ierr = TSGetDM(mprk->subts_fast,&subdm);CHKERRQ(ierr); 609 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 610 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 611 ierr = TSSetDM(mprk->subts_fast,newdm);CHKERRQ(ierr); 612 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 613 614 if (tab->np == 3) { 615 ierr = TSRHSSplitGetSubTS(ts,"medium",&mprk->subts_medium);CHKERRQ(ierr); 616 ierr = TSRHSSplitGetSubTS(ts,"mediumbuffer",&mprk->subts_mediumbuffer);CHKERRQ(ierr); 617 if (mprk->subts_medium && !mprk->subts_mediumbuffer) { 618 mprk->subts_mediumbuffer = mprk->subts_medium; 619 mprk->subts_medium = NULL; 620 } 621 if (mprk->subts_medium) { 622 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 623 ierr = TSGetDM(mprk->subts_medium,&subdm);CHKERRQ(ierr); 624 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 625 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 626 ierr = TSSetDM(mprk->subts_medium,newdm);CHKERRQ(ierr); 627 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 628 } 629 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 630 ierr = TSGetDM(mprk->subts_mediumbuffer,&subdm);CHKERRQ(ierr); 631 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 632 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 633 ierr = TSSetDM(mprk->subts_mediumbuffer,newdm);CHKERRQ(ierr); 634 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 635 } 636 PetscFunctionReturn(0); 637 } 638 639 /* 640 This if for nonsplit RHS MPRK 641 The step completion formula is 642 643 x1 = x0 + h b^T YdotRHS 644 645 */ 646 static PetscErrorCode TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool *done) 647 { 648 TS_MPRK *mprk = (TS_MPRK*)ts->data; 649 MPRKTableau tab = mprk->tableau; 650 Vec Xtmp = mprk->Ytmp,Xslow,Xfast; 651 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow; 652 PetscReal h = ts->time_step; 653 PetscInt s = tab->s,j; 654 PetscErrorCode ierr; 655 656 PetscFunctionBegin; 657 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 658 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 659 for (j=0; j<s; j++) ws[j] = h*tab->bsb[j]; 660 ierr = VecCopy(X,Xtmp);CHKERRQ(ierr); 661 ierr = VecMAXPY(Xtmp,s,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 662 ierr = VecGetSubVector(Xtmp,mprk->is_slow,&Xslow);CHKERRQ(ierr); 663 ierr = VecISCopy(X,mprk->is_slow,SCATTER_FORWARD,Xslow);CHKERRQ(ierr); 664 ierr = VecRestoreSubVector(Xtmp,mprk->is_slow,&Xslow);CHKERRQ(ierr); 665 666 /* Update fast part of X, note that the slow part has been changed but is simply discarded here */ 667 ierr = VecCopy(X,Xtmp);CHKERRQ(ierr); 668 ierr = VecMAXPY(Xtmp,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 669 ierr = VecGetSubVector(Xtmp,mprk->is_fast,&Xfast);CHKERRQ(ierr); 670 ierr = VecISCopy(X,mprk->is_fast,SCATTER_FORWARD,Xfast);CHKERRQ(ierr); 671 ierr = VecRestoreSubVector(Xtmp,mprk->is_fast,&Xfast);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 static PetscErrorCode TSStep_MPRK(TS ts) 676 { 677 TS_MPRK *mprk = (TS_MPRK*)ts->data; 678 Vec *Y = mprk->Y,Ytmp = mprk->Ytmp,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow; 679 Vec Yfast,Yslow; 680 MPRKTableau tab = mprk->tableau; 681 const PetscInt s = tab->s; 682 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 683 PetscScalar *wf = mprk->work_fast,*wsb = mprk->work_slowbuffer; 684 PetscInt i,j; 685 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 686 PetscErrorCode ierr; 687 688 PetscFunctionBegin; 689 for (i=0; i<s; i++) { 690 mprk->stage_time = t + h*cf[i]; 691 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 692 693 /* update the satge value for all components by slow and fast tableau respectively */ 694 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 695 ierr = VecCopy(ts->vec_sol,Ytmp);CHKERRQ(ierr); 696 ierr = VecMAXPY(Ytmp,i,wsb,YdotRHS_slow);CHKERRQ(ierr); 697 ierr = VecGetSubVector(Ytmp,mprk->is_slow,&Yslow);CHKERRQ(ierr); 698 ierr = VecISCopy(Y[i],mprk->is_slow,SCATTER_FORWARD,Yslow);CHKERRQ(ierr); 699 ierr = VecRestoreSubVector(Ytmp,mprk->is_slow,&Yslow);CHKERRQ(ierr); 700 701 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 702 ierr = VecCopy(ts->vec_sol,Ytmp);CHKERRQ(ierr); 703 ierr = VecMAXPY(Ytmp,i,wf,YdotRHS_fast);CHKERRQ(ierr); 704 ierr = VecGetSubVector(Ytmp,mprk->is_fast,&Yfast);CHKERRQ(ierr); 705 ierr = VecISCopy(Y[i],mprk->is_fast,SCATTER_FORWARD,Yfast);CHKERRQ(ierr); 706 ierr = VecRestoreSubVector(Ytmp,mprk->is_fast,&Yfast);CHKERRQ(ierr); 707 708 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 709 /* compute the stage RHS by fast and slow tableau respectively */ 710 ierr = TSComputeRHSFunction(ts,t+h*csb[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 711 ierr = TSComputeRHSFunction(ts,t+h*cf[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 712 } 713 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 714 ts->ptime += ts->time_step; 715 ts->time_step = next_time_step; 716 PetscFunctionReturn(0); 717 } 718 719 /* 720 This if for partitioned RHS MPRK 721 The step completion formula is 722 723 x1 = x0 + h b^T YdotRHS 724 725 */ 726 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 727 { 728 TS_MPRK *mprk = (TS_MPRK*)ts->data; 729 MPRKTableau tab = mprk->tableau; 730 Vec Xslow,Xfast,Xslowbuffer; /* subvectors for slow and fast componets in X respectively */ 731 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 732 PetscReal h = ts->time_step; 733 PetscInt s = tab->s,j,basestages; 734 PetscErrorCode ierr; 735 736 PetscFunctionBegin; 737 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 738 739 /* slow region */ 740 if (mprk->is_slow) { 741 basestages = 0; 742 for (j=0; j<s; j++) { 743 if (!tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->bsb[j]; 744 else ws[basestages++] = h*tab->bsb[j]; 745 } 746 ierr = VecGetSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 747 ierr = VecMAXPY(Xslow,basestages,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 748 ierr = VecRestoreSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 749 } 750 751 /* slow buffer region */ 752 for (j=0; j<s; j++) wsb[j] = h*tab->bsb[j]; 753 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 754 ierr = VecMAXPY(Xslowbuffer,s,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 755 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 756 757 if (tab->np == 3) { 758 Vec Xmedium,Xmediumbuffer; 759 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 760 /* medium region */ 761 if (mprk->is_medium) { 762 basestages = 0; 763 for (j=0; j<s; j++) { 764 if (!tab->rmb[j]) wm[tab->rmb[j]-1] += h*tab->bmb[j]; 765 else wm[basestages++] = h*tab->bmb[j]; 766 } 767 ierr = VecGetSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 768 ierr = VecMAXPY(Xmedium,basestages,wm,mprk->YdotRHS_medium);CHKERRQ(ierr); 769 ierr = VecRestoreSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 770 } 771 /* medium buffer region */ 772 for (j=0; j<s; j++) wmb[j] = h*tab->bmb[j]; 773 ierr = VecGetSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 774 ierr = VecMAXPY(Xmediumbuffer,s,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 775 ierr = VecRestoreSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 776 } 777 778 /* fast region */ 779 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 780 ierr = VecGetSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 781 ierr = VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 782 ierr = VecRestoreSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 783 PetscFunctionReturn(0); 784 } 785 786 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 787 { 788 TS_MPRK *mprk = (TS_MPRK*)ts->data; 789 MPRKTableau tab = mprk->tableau; 790 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 791 Vec Yslow,Yslowbuffer,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 792 PetscInt s = tab->s; 793 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 794 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 795 PetscInt i,j,basestages; 796 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 for (i=0; i<s; i++) { 801 mprk->stage_time = t + h*cf[i]; 802 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 803 /* calculate the stage value for fast and slow components respectively */ 804 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 805 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 806 807 /* slow buffer region */ 808 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 809 ierr = VecMAXPY(Yslowbuffer,i,wsb,YdotRHS_slowbuffer);CHKERRQ(ierr); 810 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 811 /* slow region */ 812 if (!tab->rsb[i] && mprk->is_slow) { /* not a repeated stage */ 813 basestages = 0; 814 for (j=0; j<i; j++) { 815 if (!tab->rsb[j]) ws[tab->rsb[j]-1] += h*Asb[i*s+j]; 816 else ws[basestages++] = wsb[j]; 817 } 818 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 819 ierr = VecMAXPY(Yslow,i,ws,YdotRHS_slow);CHKERRQ(ierr); 820 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 821 /* only depends on the slow buffer region */ 822 ierr = TSComputeRHSFunction(mprk->subts_slow,t+h*csb[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 823 } 824 825 /* fast region */ 826 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 827 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 828 ierr = VecMAXPY(Yfast,i,wf,YdotRHS_fast);CHKERRQ(ierr); 829 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 830 831 if (tab->np == 3) { 832 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 833 Vec Ymedium,Ymediumbuffer; 834 const PetscReal *Amb = tab->Amb,*cmb = tab->cmb; 835 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 836 837 for (j=0; j<i; j++) wmb[j] = h*Amb[i*s+j]; 838 /* medium buffer region */ 839 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 840 ierr = VecMAXPY(Ymediumbuffer,i,wmb,YdotRHS_mediumbuffer);CHKERRQ(ierr); 841 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 842 /* medium region */ 843 if (!tab->rmb[i] && mprk->is_medium) { /* not a repeated stage */ 844 basestages = 0; 845 for (j=0; j<i; j++) { 846 if (tab->rmb[j]) wm[tab->rmb[j]-1] += h*Amb[i*s+j]; 847 else wm[basestages++] = wmb[j]; 848 } 849 ierr = VecGetSubVector(Y[basestages],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 850 ierr = VecMAXPY(Ymedium,i,wm,YdotRHS_medium);CHKERRQ(ierr); 851 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 852 /* only depends on the medium buffer region and the slow buffer region */ 853 ierr = TSComputeRHSFunction(mprk->subts_medium,t+h*cmb[i],Y[i],YdotRHS_medium[i]);CHKERRQ(ierr); 854 } 855 /* must be computed after fast region and slow region are updated in Y */ 856 ierr = TSComputeRHSFunction(mprk->subts_mediumbuffer,t+h*cmb[i],Y[i],YdotRHS_mediumbuffer[i]);CHKERRQ(ierr); 857 } 858 /* must be computed after all regions are updated in Y */ 859 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[i]);CHKERRQ(ierr); 860 ierr = TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 861 } 862 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 863 ts->ptime += ts->time_step; 864 ts->time_step = next_time_step; 865 PetscFunctionReturn(0); 866 } 867 868 static PetscErrorCode TSMPRKTableauReset(TS ts) 869 { 870 TS_MPRK *mprk = (TS_MPRK*)ts->data; 871 MPRKTableau tab = mprk->tableau; 872 PetscErrorCode ierr; 873 874 PetscFunctionBegin; 875 if (!tab) PetscFunctionReturn(0); 876 ierr = PetscFree(mprk->work_fast);CHKERRQ(ierr); 877 ierr = PetscFree(mprk->work_slow);CHKERRQ(ierr); 878 ierr = VecDestroyVecs(tab->s,&mprk->Y);CHKERRQ(ierr); 879 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 880 ierr = VecDestroy(&mprk->Ytmp);CHKERRQ(ierr); 881 } 882 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 883 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 884 PetscFunctionReturn(0); 885 } 886 887 static PetscErrorCode TSReset_MPRK(TS ts) 888 { 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr); 893 PetscFunctionReturn(0); 894 } 895 896 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 897 { 898 PetscFunctionBegin; 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 903 { 904 PetscFunctionBegin; 905 PetscFunctionReturn(0); 906 } 907 908 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 909 { 910 PetscFunctionBegin; 911 PetscFunctionReturn(0); 912 } 913 914 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 915 { 916 PetscFunctionBegin; 917 PetscFunctionReturn(0); 918 } 919 920 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 921 { 922 TS_MPRK *mprk = (TS_MPRK*)ts->data; 923 MPRKTableau tab = mprk->tableau; 924 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y);CHKERRQ(ierr); 929 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 930 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 931 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 932 ierr = VecDuplicate(ts->vec_sol,&mprk->Ytmp);CHKERRQ(ierr); 933 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 934 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 935 } 936 if (mprk->mprkmtype == TSMPRKSPLIT) { 937 if (mprk->is_slow) { 938 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 939 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 940 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 941 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 942 } 943 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 944 ierr = VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 945 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 946 ierr = PetscMalloc1(tab->s,&mprk->work_slowbuffer);CHKERRQ(ierr); 947 948 if (tab->np == 3) { 949 if (mprk->is_medium) { 950 ierr = VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 951 ierr = VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 952 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 953 ierr = PetscMalloc1(tab->s,&mprk->work_medium);CHKERRQ(ierr); 954 } 955 ierr = VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 956 ierr = VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 957 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 958 ierr = PetscMalloc1(tab->s,&mprk->work_mediumbuffer);CHKERRQ(ierr); 959 } 960 ierr = VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 961 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 962 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 963 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 964 } 965 PetscFunctionReturn(0); 966 } 967 968 static PetscErrorCode TSSetUp_MPRK(TS ts) 969 { 970 TS_MPRK *mprk = (TS_MPRK*)ts->data; 971 MPRKTableau tab = mprk->tableau; 972 DM dm; 973 PetscErrorCode ierr; 974 975 PetscFunctionBegin; 976 ierr = TSRHSSplitGetIS(ts,"slow",&mprk->is_slow);CHKERRQ(ierr); 977 ierr = TSRHSSplitGetIS(ts,"fast",&mprk->is_fast);CHKERRQ(ierr); 978 if (!mprk->is_slow || !mprk->is_fast) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'fast' respectively in order to use the method '%s'",tab->name); 979 980 if (tab->np == 3) { 981 ierr = TSRHSSplitGetIS(ts,"medium",&mprk->is_medium);CHKERRQ(ierr); 982 if (!mprk->is_medium) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'medium' and 'fast' respectively in order to use the method '%s'",tab->name); 983 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 984 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 985 mprk->is_mediumbuffer = mprk->is_medium; 986 mprk->is_medium = NULL; 987 } 988 } 989 990 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 991 ierr = TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer);CHKERRQ(ierr); 992 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 993 mprk->is_slowbuffer = mprk->is_slow; 994 mprk->is_slow = NULL; 995 } 996 /* 997 if (!mprk->is_medium) { 998 mprk->is_medium = mprk->is_fast; 999 mprk->is_fast = NULL; 1000 } else { 1001 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 1002 } 1003 */ 1004 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1005 ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr); 1006 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1007 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1008 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1009 ierr = PetscTryMethod(ts,"TSMPRKSetSplits_C",(TS),(ts));CHKERRQ(ierr); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 /* construct a database to chose nonsplit RHS mutirate mprk method or split RHS MPRK method */ 1014 const char *const TSMPRKMultirateTypes[] = {"NONSPLIT","SPLIT","TSMPRKMultirateType","TSMPRK",0}; 1015 1016 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1017 { 1018 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1019 PetscErrorCode ierr; 1020 1021 PetscFunctionBegin; 1022 ierr = PetscOptionsHead(PetscOptionsObject,"PRK ODE solver options");CHKERRQ(ierr); 1023 { 1024 MPRKTableauLink link; 1025 PetscInt count,choice; 1026 PetscBool flg; 1027 const char **namelist; 1028 PetscInt mprkmtype = 0; 1029 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1030 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1031 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1032 ierr = PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1033 if (flg) {ierr = TSMPRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1034 ierr = PetscFree(namelist);CHKERRQ(ierr); 1035 ierr = PetscOptionsEList("-ts_mprk_multirate_type","Use Combined RHS Multirate or Partioned RHS Multirat MPRK method","TSMPRKSetMultirateType",TSMPRKMultirateTypes,2,TSMPRKMultirateTypes[0],&mprkmtype,&flg);CHKERRQ(ierr); 1036 if (flg) {ierr = TSMPRKSetMultirateType(ts,mprkmtype);CHKERRQ(ierr);} 1037 } 1038 ierr = PetscOptionsTail();CHKERRQ(ierr); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1043 { 1044 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1045 PetscBool iascii; 1046 PetscErrorCode ierr; 1047 1048 PetscFunctionBegin; 1049 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1050 if (iascii) { 1051 MPRKTableau tab = mprk->tableau; 1052 TSMPRKType mprktype; 1053 char fbuf[512]; 1054 char sbuf[512]; 1055 PetscInt i; 1056 ierr = TSMPRKGetType(ts,&mprktype);CHKERRQ(ierr); 1057 ierr = PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype);CHKERRQ(ierr); 1058 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1059 1060 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf);CHKERRQ(ierr); 1061 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf);CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," Af = \n");CHKERRQ(ierr); 1063 for (i=0; i<tab->s; i++) { 1064 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s]);CHKERRQ(ierr); 1065 ierr = PetscViewerASCIIPrintf(viewer," %s\n",fbuf);CHKERRQ(ierr); 1066 } 1067 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf);CHKERRQ(ierr); 1068 ierr = PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf);CHKERRQ(ierr); 1069 1070 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb);CHKERRQ(ierr); 1071 ierr = PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf);CHKERRQ(ierr); 1072 ierr = PetscViewerASCIIPrintf(viewer," Asb = \n");CHKERRQ(ierr); 1073 for (i=0; i<tab->s; i++) { 1074 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s]);CHKERRQ(ierr); 1075 ierr = PetscViewerASCIIPrintf(viewer," %s\n",sbuf);CHKERRQ(ierr); 1076 } 1077 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb);CHKERRQ(ierr); 1078 ierr = PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf);CHKERRQ(ierr); 1079 1080 if (tab->np == 3) { 1081 char mbuf[512]; 1082 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb);CHKERRQ(ierr); 1083 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf);CHKERRQ(ierr); 1084 ierr = PetscViewerASCIIPrintf(viewer," Amb = \n");CHKERRQ(ierr); 1085 for (i=0; i<tab->s; i++) { 1086 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s]);CHKERRQ(ierr); 1087 ierr = PetscViewerASCIIPrintf(viewer," %s\n",mbuf);CHKERRQ(ierr); 1088 } 1089 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb);CHKERRQ(ierr); 1090 ierr = PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf);CHKERRQ(ierr); 1091 } 1092 } 1093 PetscFunctionReturn(0); 1094 } 1095 1096 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1097 { 1098 PetscErrorCode ierr; 1099 TSAdapt adapt; 1100 1101 PetscFunctionBegin; 1102 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1103 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /*@C 1108 TSMPRKSetType - Set the type of MPRK scheme 1109 1110 Logically collective 1111 1112 Input Parameter: 1113 + ts - timestepping context 1114 - mprktype - type of MPRK-scheme 1115 1116 Options Database: 1117 . -ts_mprk_type - <pm2,p2,p3> 1118 1119 Level: intermediate 1120 1121 .seealso: TSMPRKGetType(), TSMPRK, TSMPRKType 1122 @*/ 1123 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1124 { 1125 PetscErrorCode ierr; 1126 1127 PetscFunctionBegin; 1128 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1129 PetscValidCharPointer(mprktype,2); 1130 ierr = PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype));CHKERRQ(ierr); 1131 PetscFunctionReturn(0); 1132 } 1133 1134 /*@C 1135 TSMPRKGetType - Get the type of MPRK scheme 1136 1137 Logically collective 1138 1139 Input Parameter: 1140 . ts - timestepping context 1141 1142 Output Parameter: 1143 . mprktype - type of MPRK-scheme 1144 1145 Level: intermediate 1146 1147 .seealso: TSMPRKGetType() 1148 @*/ 1149 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1150 { 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1155 ierr = PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype));CHKERRQ(ierr); 1156 PetscFunctionReturn(0); 1157 } 1158 1159 /*@C 1160 TSMPRKSetMultirateType - Set the type of MPRK multirate scheme 1161 1162 Logically collective 1163 1164 Input Parameter: 1165 + ts - timestepping context 1166 - mprkmtype - type of the multirate configuration 1167 1168 Options Database: 1169 . -ts_mprk_multirate_type - <nonsplit,split> 1170 1171 Level: intermediate 1172 @*/ 1173 PetscErrorCode TSMPRKSetMultirateType(TS ts, TSMPRKMultirateType mprkmtype) 1174 { 1175 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1176 PetscErrorCode ierr; 1177 1178 PetscFunctionBegin; 1179 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1180 switch(mprkmtype){ 1181 case TSMPRKNONSPLIT: 1182 ts->ops->step = TSStep_MPRK; 1183 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1184 break; 1185 case TSMPRKSPLIT: 1186 ts->ops->step = TSStep_MPRKSPLIT; 1187 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1188 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetSplits_C",TSMPRKSetSplits);CHKERRQ(ierr); 1189 break; 1190 default : 1191 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",mprkmtype); 1192 } 1193 mprk->mprkmtype = mprkmtype; 1194 PetscFunctionReturn(0); 1195 } 1196 1197 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1198 { 1199 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1200 1201 PetscFunctionBegin; 1202 *mprktype = mprk->tableau->name; 1203 PetscFunctionReturn(0); 1204 } 1205 1206 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1207 { 1208 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1209 PetscBool match; 1210 MPRKTableauLink link; 1211 PetscErrorCode ierr; 1212 1213 PetscFunctionBegin; 1214 if (mprk->tableau) { 1215 ierr = PetscStrcmp(mprk->tableau->name,mprktype,&match);CHKERRQ(ierr); 1216 if (match) PetscFunctionReturn(0); 1217 } 1218 for (link = MPRKTableauList; link; link=link->next) { 1219 ierr = PetscStrcmp(link->tab.name,mprktype,&match);CHKERRQ(ierr); 1220 if (match) { 1221 if (ts->setupcalled) {ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr);} 1222 mprk->tableau = &link->tab; 1223 if (ts->setupcalled) {ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr);} 1224 PetscFunctionReturn(0); 1225 } 1226 } 1227 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1232 { 1233 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1234 1235 PetscFunctionBegin; 1236 *ns = mprk->tableau->s; 1237 if (Y) *Y = mprk->Y; 1238 PetscFunctionReturn(0); 1239 } 1240 1241 static PetscErrorCode TSDestroy_MPRK(TS ts) 1242 { 1243 PetscErrorCode ierr; 1244 1245 PetscFunctionBegin; 1246 ierr = TSReset_MPRK(ts);CHKERRQ(ierr); 1247 if (ts->dm) { 1248 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1249 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1250 } 1251 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1252 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL);CHKERRQ(ierr); 1253 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL);CHKERRQ(ierr); 1254 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",NULL);CHKERRQ(ierr); 1255 PetscFunctionReturn(0); 1256 } 1257 1258 /*MC 1259 TSMPRK - ODE solver using Partitioned Runge-Kutta schemes 1260 1261 The user should provide the right hand side of the equation 1262 using TSSetRHSFunction(). 1263 1264 Notes: 1265 The default is TSMPRKPM2, it can be changed with TSRKSetType() or -ts_mprk_type 1266 1267 Level: beginner 1268 1269 .seealso: TSCreate(), TS, TSSetType(), TSMPRKSetType(), TSMPRKGetType(), TSMPRKType, TSMPRKRegister(), TSMPRKSetMultirateType() 1270 TSMPRKM2, TSMPRKM3, TSMPRKRFSMR3, TSMPRKRFSMR2 1271 1272 M*/ 1273 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1274 { 1275 TS_MPRK *mprk; 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 ierr = TSMPRKInitializePackage();CHKERRQ(ierr); 1280 1281 ts->ops->reset = TSReset_MPRK; 1282 ts->ops->destroy = TSDestroy_MPRK; 1283 ts->ops->view = TSView_MPRK; 1284 ts->ops->load = TSLoad_MPRK; 1285 ts->ops->setup = TSSetUp_MPRK; 1286 ts->ops->step = TSStep_MPRK; 1287 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1288 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1289 ts->ops->getstages = TSGetStages_MPRK; 1290 1291 ierr = PetscNewLog(ts,&mprk);CHKERRQ(ierr); 1292 ts->data = (void*)mprk; 1293 1294 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK);CHKERRQ(ierr); 1295 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK);CHKERRQ(ierr); 1296 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",TSMPRKSetMultirateType);CHKERRQ(ierr); 1297 1298 ierr = TSMPRKSetType(ts,TSMPRKDefault);CHKERRQ(ierr); 1299 PetscFunctionReturn(0); 1300 } 1301