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