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 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 Multirate 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 Multirate 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 Multirate 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 Multirate 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 Multirate 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 Multirate 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,2,2,1,&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,2,2,2,&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,2,3,1,&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,2,3,3,&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,1,1,&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,5,2,1,&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 TSMPRKFinalizePackage - 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 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 .keywords: TS, register 482 483 .seealso: TSMPRK 484 @*/ 485 PetscErrorCode TSMPRKRegister(TSMPRKType name,PetscInt order, 486 PetscInt sbase,PetscInt ratio1,PetscInt ratio2, 487 const PetscReal Asb[],const PetscReal bsb[],const PetscReal csb[],const PetscInt rsb[], 488 const PetscReal Amb[],const PetscReal bmb[],const PetscReal cmb[],const PetscInt rmb[], 489 const PetscReal Af[],const PetscReal bf[],const PetscReal cf[]) 490 { 491 MPRKTableauLink link; 492 MPRKTableau t; 493 PetscInt s,i,j; 494 PetscErrorCode ierr; 495 496 PetscFunctionBegin; 497 PetscValidCharPointer(name,1); 498 PetscValidRealPointer(Asb,4); 499 if (bsb) PetscValidRealPointer(bsb,7); 500 if (csb) PetscValidRealPointer(csb,8); 501 if (rsb) PetscValidRealPointer(rsb,9); 502 if (Amb) PetscValidRealPointer(Amb,10); 503 if (bmb) PetscValidRealPointer(bmb,11); 504 if (cmb) PetscValidRealPointer(cmb,12); 505 if (rmb) PetscValidRealPointer(rmb,13); 506 PetscValidRealPointer(Af,14); 507 if (bf) PetscValidRealPointer(bf,15); 508 if (cf) PetscValidRealPointer(cf,16); 509 510 ierr = PetscNew(&link);CHKERRQ(ierr); 511 t = &link->tab; 512 513 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 514 s = sbase*ratio1*ratio2; /* this is the dimension of the matrices below */ 515 t->order = order; 516 t->sbase = sbase; 517 t->s = s; 518 t->np = 2; 519 520 ierr = PetscMalloc3(s*s,&t->Af,s,&t->bf,s,&t->cf);CHKERRQ(ierr); 521 ierr = PetscMemcpy(t->Af,Af,s*s*sizeof(Af[0]));CHKERRQ(ierr); 522 if (bf) { 523 ierr = PetscMemcpy(t->bf,bf,s*sizeof(bf[0]));CHKERRQ(ierr); 524 } else 525 for (i=0; i<s; i++) t->bf[i] = Af[(s-1)*s+i]; 526 if (cf) { 527 ierr = PetscMemcpy(t->cf,cf,s*sizeof(cf[0]));CHKERRQ(ierr); 528 } else { 529 for (i=0; i<s; i++) 530 for (j=0,t->cf[i]=0; j<s; j++) 531 t->cf[i] += Af[i*s+j]; 532 } 533 534 if (Amb) { 535 t->np = 3; 536 ierr = PetscMalloc3(s*s,&t->Amb,s,&t->bmb,s,&t->cmb);CHKERRQ(ierr); 537 ierr = PetscCalloc1(s,&t->rmb);CHKERRQ(ierr); 538 ierr = PetscMemcpy(t->Amb,Amb,s*s*sizeof(Amb[0]));CHKERRQ(ierr); 539 if (bmb) { 540 ierr = PetscMemcpy(t->bmb,bmb,s*sizeof(bmb[0]));CHKERRQ(ierr); 541 } else { 542 for (i=0; i<s; i++) t->bmb[i] = Amb[(s-1)*s+i]; 543 } 544 if (cmb) { 545 ierr = PetscMemcpy(t->cmb,cmb,s*sizeof(cmb[0]));CHKERRQ(ierr); 546 } else { 547 for (i=0; i<s; i++) 548 for (j=0,t->cmb[i]=0; j<s; j++) 549 t->cmb[i] += Amb[i*s+j]; 550 } 551 if (rmb) { 552 ierr = PetscMemcpy(t->rmb,rmb,s*sizeof(rmb[0]));CHKERRQ(ierr); 553 } 554 } 555 556 ierr = PetscMalloc3(s*s,&t->Asb,s,&t->bsb,s,&t->csb);CHKERRQ(ierr); 557 ierr = PetscMemcpy(t->Asb,Asb,s*s*sizeof(Asb[0]));CHKERRQ(ierr); 558 ierr = PetscCalloc1(s,&t->rsb);CHKERRQ(ierr); 559 if (bsb) { 560 ierr = PetscMemcpy(t->bsb,bsb,s*sizeof(bsb[0]));CHKERRQ(ierr); 561 } else 562 for (i=0; i<s; i++) t->bsb[i] = Asb[(s-1)*s+i]; 563 if (csb) { 564 ierr = PetscMemcpy(t->csb,csb,s*sizeof(csb[0]));CHKERRQ(ierr); 565 } else { 566 for (i=0; i<s; i++) 567 for (j=0,t->csb[i]=0; j<s; j++) 568 t->csb[i] += Asb[i*s+j]; 569 } 570 if (rsb) { 571 ierr = PetscMemcpy(t->rsb,rsb,s*sizeof(rsb[0]));CHKERRQ(ierr); 572 } 573 link->next = MPRKTableauList; 574 MPRKTableauList = link; 575 PetscFunctionReturn(0); 576 } 577 578 static PetscErrorCode TSMPRKSetSplits(TS ts) 579 { 580 TS_MPRK *mprk = (TS_MPRK*)ts->data; 581 MPRKTableau tab = mprk->tableau; 582 DM dm,subdm,newdm; 583 PetscErrorCode ierr; 584 585 PetscFunctionBegin; 586 ierr = TSRHSSplitGetSubTS(ts,"slow",&mprk->subts_slow);CHKERRQ(ierr); 587 ierr = TSRHSSplitGetSubTS(ts,"fast",&mprk->subts_fast);CHKERRQ(ierr); 588 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"); 589 590 /* Only copy the DM */ 591 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 592 593 ierr = TSRHSSplitGetSubTS(ts,"slowbuffer",&mprk->subts_slowbuffer);CHKERRQ(ierr); 594 if (!mprk->subts_slowbuffer) { 595 mprk->subts_slowbuffer = mprk->subts_slow; 596 mprk->subts_slow = NULL; 597 } 598 if (mprk->subts_slow) { 599 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 600 ierr = TSGetDM(mprk->subts_slow,&subdm);CHKERRQ(ierr); 601 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 602 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 603 ierr = TSSetDM(mprk->subts_slow,newdm);CHKERRQ(ierr); 604 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 605 } 606 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 607 ierr = TSGetDM(mprk->subts_slowbuffer,&subdm);CHKERRQ(ierr); 608 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 609 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 610 ierr = TSSetDM(mprk->subts_slowbuffer,newdm);CHKERRQ(ierr); 611 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 612 613 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 614 ierr = TSGetDM(mprk->subts_fast,&subdm);CHKERRQ(ierr); 615 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 616 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 617 ierr = TSSetDM(mprk->subts_fast,newdm);CHKERRQ(ierr); 618 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 619 620 if (tab->np == 3) { 621 ierr = TSRHSSplitGetSubTS(ts,"medium",&mprk->subts_medium);CHKERRQ(ierr); 622 ierr = TSRHSSplitGetSubTS(ts,"mediumbuffer",&mprk->subts_mediumbuffer);CHKERRQ(ierr); 623 if (mprk->subts_medium && !mprk->subts_mediumbuffer) { 624 mprk->subts_mediumbuffer = mprk->subts_medium; 625 mprk->subts_medium = NULL; 626 } 627 if (mprk->subts_medium) { 628 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 629 ierr = TSGetDM(mprk->subts_medium,&subdm);CHKERRQ(ierr); 630 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 631 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 632 ierr = TSSetDM(mprk->subts_medium,newdm);CHKERRQ(ierr); 633 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 634 } 635 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 636 ierr = TSGetDM(mprk->subts_mediumbuffer,&subdm);CHKERRQ(ierr); 637 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 638 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 639 ierr = TSSetDM(mprk->subts_mediumbuffer,newdm);CHKERRQ(ierr); 640 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 641 } 642 PetscFunctionReturn(0); 643 } 644 645 /* 646 This if for nonsplit RHS MPRK 647 The step completion formula is 648 649 x1 = x0 + h b^T YdotRHS 650 651 */ 652 static PetscErrorCode TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool *done) 653 { 654 TS_MPRK *mprk = (TS_MPRK*)ts->data; 655 MPRKTableau tab = mprk->tableau; 656 PetscScalar *wf = mprk->work_fast; 657 PetscReal h = ts->time_step; 658 PetscInt s = tab->s,j; 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 663 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 664 ierr = VecMAXPY(X,s,wf,mprk->YdotRHS);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 static PetscErrorCode TSStep_MPRK(TS ts) 669 { 670 TS_MPRK *mprk = (TS_MPRK*)ts->data; 671 Vec *Y = mprk->Y,*YdotRHS = mprk->YdotRHS,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 672 Vec Yslow,Yslowbuffer,Yfast; 673 MPRKTableau tab = mprk->tableau; 674 const PetscInt s = tab->s; 675 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 676 PetscScalar *wf = mprk->work_fast,*wsb = mprk->work_slowbuffer; 677 PetscInt i,j; 678 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 679 PetscErrorCode ierr; 680 681 PetscFunctionBegin; 682 for (i=0; i<s; i++) { 683 mprk->stage_time = t + h*cf[i]; 684 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 685 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 686 687 /* slow buffer region */ 688 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 689 for (j=0; j<i; j++) { 690 ierr = VecGetSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j]);CHKERRQ(ierr); 691 } 692 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 693 ierr = VecMAXPY(Yslowbuffer,i,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 694 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 695 for (j=0; j<i; j++) { 696 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j]);CHKERRQ(ierr); 697 } 698 /* slow region */ 699 if (mprk->is_slow) { 700 for (j=0; j<i; j++) { 701 ierr = VecGetSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j]);CHKERRQ(ierr); 702 } 703 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 704 ierr = VecMAXPY(Yslow,i,wsb,mprk->YdotRHS_slow);CHKERRQ(ierr); 705 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 706 for (j=0; j<i; j++) { 707 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j]);CHKERRQ(ierr); 708 } 709 } 710 711 /* fast region */ 712 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 713 for (j=0; j<i; j++) { 714 ierr = VecGetSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j]);CHKERRQ(ierr); 715 } 716 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 717 ierr = VecMAXPY(Yfast,i,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 718 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 719 for (j=0; j<i; j++) { 720 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j]);CHKERRQ(ierr); 721 } 722 if (tab->np == 3) { 723 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 724 Vec Ymedium,Ymediumbuffer; 725 PetscScalar *wmb = mprk->work_mediumbuffer; 726 727 for (j=0; j<i; j++) wmb[j] = h*tab->Amb[i*s+j]; 728 /* medium region */ 729 if (mprk->is_medium) { 730 for (j=0; j<i; j++) { 731 ierr = VecGetSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 732 } 733 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 734 ierr = VecMAXPY(Ymedium,i,wmb,mprk->YdotRHS_medium);CHKERRQ(ierr); 735 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 736 for (j=0; j<i; j++) { 737 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 738 } 739 } 740 /* medium buffer region */ 741 for (j=0; j<i; j++) { 742 ierr = VecGetSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 743 } 744 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 745 ierr = VecMAXPY(Ymediumbuffer,i,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 746 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 747 for (j=0; j<i; j++) { 748 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 749 } 750 } 751 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 752 /* compute the stage RHS by fast and slow tableau respectively */ 753 ierr = TSComputeRHSFunction(ts,t+h*csb[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 754 } 755 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 756 ts->ptime += ts->time_step; 757 ts->time_step = next_time_step; 758 PetscFunctionReturn(0); 759 } 760 761 /* 762 This if for the case when split RHS is used 763 The step completion formula is 764 x1 = x0 + h b^T YdotRHS 765 */ 766 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 767 { 768 TS_MPRK *mprk = (TS_MPRK*)ts->data; 769 MPRKTableau tab = mprk->tableau; 770 Vec Xslow,Xfast,Xslowbuffer; /* subvectors for slow and fast componets in X respectively */ 771 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 772 PetscReal h = ts->time_step; 773 PetscInt s = tab->s,j,computedstages; 774 PetscErrorCode ierr; 775 776 PetscFunctionBegin; 777 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 778 779 /* slow region */ 780 if (mprk->is_slow) { 781 computedstages = 0; 782 for (j=0; j<s; j++) { 783 if (tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->bsb[j]; 784 else ws[computedstages++] = h*tab->bsb[j]; 785 } 786 ierr = VecGetSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 787 ierr = VecMAXPY(Xslow,computedstages,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 788 ierr = VecRestoreSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 789 } 790 791 if (tab->np == 3 && mprk->is_medium) { 792 computedstages = 0; 793 for (j=0; j<s; j++) { 794 if (tab->rmb[j]) wsb[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bsb[j]; 795 else wsb[computedstages++] = h*tab->bsb[j]; 796 } 797 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 798 ierr = VecMAXPY(Xslowbuffer,computedstages,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 799 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 800 } else { 801 /* slow buffer region */ 802 for (j=0; j<s; j++) wsb[j] = h*tab->bsb[j]; 803 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 804 ierr = VecMAXPY(Xslowbuffer,s,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 805 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 806 } 807 if (tab->np == 3) { 808 Vec Xmedium,Xmediumbuffer; 809 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 810 /* medium region and slow buffer region */ 811 if (mprk->is_medium) { 812 computedstages = 0; 813 for (j=0; j<s; j++) { 814 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bmb[j]; 815 else wm[computedstages++] = h*tab->bmb[j]; 816 } 817 ierr = VecGetSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 818 ierr = VecMAXPY(Xmedium,computedstages,wm,mprk->YdotRHS_medium);CHKERRQ(ierr); 819 ierr = VecRestoreSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 820 } 821 /* medium buffer region */ 822 for (j=0; j<s; j++) wmb[j] = h*tab->bmb[j]; 823 ierr = VecGetSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 824 ierr = VecMAXPY(Xmediumbuffer,s,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 825 ierr = VecRestoreSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 826 } 827 /* fast region */ 828 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 829 ierr = VecGetSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 830 ierr = VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 831 ierr = VecRestoreSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 832 PetscFunctionReturn(0); 833 } 834 835 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 836 { 837 TS_MPRK *mprk = (TS_MPRK*)ts->data; 838 MPRKTableau tab = mprk->tableau; 839 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 840 Vec Yslow,Yslowbuffer,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 841 PetscInt s = tab->s; 842 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 843 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 844 PetscInt i,j,computedstages; 845 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 846 PetscErrorCode ierr; 847 848 PetscFunctionBegin; 849 for (i=0; i<s; i++) { 850 mprk->stage_time = t + h*cf[i]; 851 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 852 /* calculate the stage value for fast and slow components respectively */ 853 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 854 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 855 856 /* slow buffer region */ 857 if (tab->np == 3 && mprk->is_medium) { 858 if (tab->rmb[i]) { 859 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 860 ierr = VecISCopy(Y[tab->rmb[i]-1],mprk->is_slowbuffer,SCATTER_REVERSE,Yslowbuffer);CHKERRQ(ierr); 861 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 862 } else { 863 PetscScalar *wm = mprk->work_medium; 864 computedstages = 0; 865 for (j=0; j<i; j++) { 866 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wsb[j]; 867 else wm[computedstages++] = wsb[j]; 868 } 869 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 870 ierr = VecMAXPY(Yslowbuffer,computedstages,wm,YdotRHS_slowbuffer);CHKERRQ(ierr); 871 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 872 } 873 } else { 874 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 875 ierr = VecMAXPY(Yslowbuffer,i,wsb,YdotRHS_slowbuffer);CHKERRQ(ierr); 876 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 877 } 878 879 /* slow region */ 880 if (mprk->is_slow) { 881 if (tab->rsb[i]) { /* repeat previous stage */ 882 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 883 ierr = VecISCopy(Y[tab->rsb[i]-1],mprk->is_slow,SCATTER_REVERSE,Yslow);CHKERRQ(ierr); 884 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 885 } else { 886 computedstages = 0; 887 for (j=0; j<i; j++) { 888 if (tab->rsb[j]) ws[tab->rsb[j]-1] += wsb[j]; 889 else ws[computedstages++] = wsb[j]; 890 } 891 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 892 ierr = VecMAXPY(Yslow,computedstages,ws,YdotRHS_slow);CHKERRQ(ierr); 893 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 894 /* only depends on the slow buffer region */ 895 ierr = TSComputeRHSFunction(mprk->subts_slow,t+h*csb[i],Y[i],YdotRHS_slow[computedstages]);CHKERRQ(ierr); 896 } 897 } 898 899 /* fast region */ 900 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 901 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 902 ierr = VecMAXPY(Yfast,i,wf,YdotRHS_fast);CHKERRQ(ierr); 903 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 904 905 if (tab->np == 3) { 906 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 907 Vec Ymedium,Ymediumbuffer; 908 const PetscReal *Amb = tab->Amb,*cmb = tab->cmb; 909 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 910 911 for (j=0; j<i; j++) wmb[j] = h*Amb[i*s+j]; 912 /* medium buffer region */ 913 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 914 ierr = VecMAXPY(Ymediumbuffer,i,wmb,YdotRHS_mediumbuffer);CHKERRQ(ierr); 915 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 916 /* medium region */ 917 if (mprk->is_medium) { 918 if (tab->rmb[i]) { /* repeat previous stage */ 919 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 920 ierr = VecISCopy(Y[tab->rmb[i]-1],mprk->is_medium,SCATTER_REVERSE,Ymedium);CHKERRQ(ierr); 921 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 922 } else { 923 computedstages = 0; 924 for (j=0; j<i; j++) { 925 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wmb[j]; 926 else wm[computedstages++] = wmb[j]; 927 928 } 929 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 930 ierr = VecMAXPY(Ymedium,computedstages,wm,YdotRHS_medium);CHKERRQ(ierr); 931 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 932 /* only depends on the medium buffer region */ 933 ierr = TSComputeRHSFunction(mprk->subts_medium,t+h*cmb[i],Y[i],YdotRHS_medium[computedstages]);CHKERRQ(ierr); 934 /* must be computed after all regions are updated in Y */ 935 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[computedstages]);CHKERRQ(ierr); 936 } 937 } 938 /* must be computed after fast region and slow region are updated in Y */ 939 ierr = TSComputeRHSFunction(mprk->subts_mediumbuffer,t+h*cmb[i],Y[i],YdotRHS_mediumbuffer[i]);CHKERRQ(ierr); 940 } 941 if (!(tab->np == 3 && mprk->is_medium)) { 942 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[i]);CHKERRQ(ierr); 943 } 944 ierr = TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i]); 945 } 946 947 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 948 ts->ptime += ts->time_step; 949 ts->time_step = next_time_step; 950 PetscFunctionReturn(0); 951 } 952 953 static PetscErrorCode TSMPRKTableauReset(TS ts) 954 { 955 TS_MPRK *mprk = (TS_MPRK*)ts->data; 956 MPRKTableau tab = mprk->tableau; 957 PetscErrorCode ierr; 958 959 PetscFunctionBegin; 960 if (!tab) PetscFunctionReturn(0); 961 ierr = PetscFree(mprk->work_fast);CHKERRQ(ierr); 962 ierr = PetscFree(mprk->work_slow);CHKERRQ(ierr); 963 ierr = PetscFree(mprk->work_slowbuffer);CHKERRQ(ierr); 964 ierr = PetscFree(mprk->work_medium);CHKERRQ(ierr); 965 ierr = PetscFree(mprk->work_mediumbuffer);CHKERRQ(ierr); 966 ierr = VecDestroyVecs(tab->s,&mprk->Y);CHKERRQ(ierr); 967 if (ts->use_splitrhsfunction) { 968 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 969 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 970 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 971 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 972 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 973 } else { 974 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 975 if (mprk->is_slow) { 976 ierr = PetscFree(mprk->YdotRHS_slow);CHKERRQ(ierr); 977 } 978 ierr = PetscFree(mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 979 if (tab->np == 3) { 980 if (mprk->is_medium) { 981 ierr = PetscFree(mprk->YdotRHS_medium);CHKERRQ(ierr); 982 } 983 ierr = PetscFree(mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 984 } 985 ierr = PetscFree(mprk->YdotRHS_fast);CHKERRQ(ierr); 986 } 987 PetscFunctionReturn(0); 988 } 989 990 static PetscErrorCode TSReset_MPRK(TS ts) 991 { 992 PetscErrorCode ierr; 993 994 PetscFunctionBegin; 995 ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr); 996 PetscFunctionReturn(0); 997 } 998 999 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 1000 { 1001 PetscFunctionBegin; 1002 PetscFunctionReturn(0); 1003 } 1004 1005 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1006 { 1007 PetscFunctionBegin; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 1012 { 1013 PetscFunctionBegin; 1014 PetscFunctionReturn(0); 1015 } 1016 1017 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1018 { 1019 PetscFunctionBegin; 1020 PetscFunctionReturn(0); 1021 } 1022 1023 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 1024 { 1025 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1026 MPRKTableau tab = mprk->tableau; 1027 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 1028 PetscErrorCode ierr; 1029 1030 PetscFunctionBegin; 1031 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y);CHKERRQ(ierr); 1032 if (mprk->is_slow) { 1033 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 1034 } 1035 ierr = PetscMalloc1(tab->s,&mprk->work_slowbuffer);CHKERRQ(ierr); 1036 if (tab->np == 3) { 1037 if (mprk->is_medium) { 1038 ierr = PetscMalloc1(tab->s,&mprk->work_medium);CHKERRQ(ierr); 1039 } 1040 ierr = PetscMalloc1(tab->s,&mprk->work_mediumbuffer);CHKERRQ(ierr); 1041 } 1042 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 1043 1044 if (ts->use_splitrhsfunction) { 1045 if (mprk->is_slow) { 1046 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1047 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1048 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1049 } 1050 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1051 ierr = VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1052 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1053 if (tab->np == 3) { 1054 if (mprk->is_medium) { 1055 ierr = VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1056 ierr = VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1057 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1058 } 1059 ierr = VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1060 ierr = VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1061 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1062 } 1063 ierr = VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1064 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1065 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1066 } else { 1067 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 1068 if (mprk->is_slow) { 1069 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1070 } 1071 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1072 if (tab->np == 3) { 1073 if (mprk->is_medium) { 1074 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1075 } 1076 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1077 } 1078 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1079 } 1080 PetscFunctionReturn(0); 1081 } 1082 1083 static PetscErrorCode TSSetUp_MPRK(TS ts) 1084 { 1085 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1086 MPRKTableau tab = mprk->tableau; 1087 DM dm; 1088 PetscErrorCode ierr; 1089 1090 PetscFunctionBegin; 1091 ierr = TSRHSSplitGetIS(ts,"slow",&mprk->is_slow);CHKERRQ(ierr); 1092 ierr = TSRHSSplitGetIS(ts,"fast",&mprk->is_fast);CHKERRQ(ierr); 1093 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); 1094 1095 if (tab->np == 3) { 1096 ierr = TSRHSSplitGetIS(ts,"medium",&mprk->is_medium);CHKERRQ(ierr); 1097 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); 1098 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 1099 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1100 mprk->is_mediumbuffer = mprk->is_medium; 1101 mprk->is_medium = NULL; 1102 } 1103 } 1104 1105 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1106 ierr = TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer);CHKERRQ(ierr); 1107 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1108 mprk->is_slowbuffer = mprk->is_slow; 1109 mprk->is_slow = NULL; 1110 } 1111 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1112 ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr); 1113 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1114 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1115 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1116 if (ts->use_splitrhsfunction) { 1117 ts->ops->step = TSStep_MPRKSPLIT; 1118 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1119 ierr = TSMPRKSetSplits(ts);CHKERRQ(ierr); 1120 } else { 1121 ts->ops->step = TSStep_MPRK; 1122 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1123 } 1124 PetscFunctionReturn(0); 1125 } 1126 1127 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1128 { 1129 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscOptionsHead(PetscOptionsObject,"PRK ODE solver options");CHKERRQ(ierr); 1134 { 1135 MPRKTableauLink link; 1136 PetscInt count,choice; 1137 PetscBool flg; 1138 const char **namelist; 1139 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1140 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1141 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1142 ierr = PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1143 if (flg) {ierr = TSMPRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1144 ierr = PetscFree(namelist);CHKERRQ(ierr); 1145 } 1146 ierr = PetscOptionsTail();CHKERRQ(ierr); 1147 PetscFunctionReturn(0); 1148 } 1149 1150 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1151 { 1152 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1153 PetscBool iascii; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1158 if (iascii) { 1159 MPRKTableau tab = mprk->tableau; 1160 TSMPRKType mprktype; 1161 char fbuf[512]; 1162 char sbuf[512]; 1163 PetscInt i; 1164 ierr = TSMPRKGetType(ts,&mprktype);CHKERRQ(ierr); 1165 ierr = PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype);CHKERRQ(ierr); 1166 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1167 1168 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf);CHKERRQ(ierr); 1169 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf);CHKERRQ(ierr); 1170 ierr = PetscViewerASCIIPrintf(viewer," Af = \n");CHKERRQ(ierr); 1171 for (i=0; i<tab->s; i++) { 1172 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s]);CHKERRQ(ierr); 1173 ierr = PetscViewerASCIIPrintf(viewer," %s\n",fbuf);CHKERRQ(ierr); 1174 } 1175 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf);CHKERRQ(ierr); 1176 ierr = PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf);CHKERRQ(ierr); 1177 1178 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb);CHKERRQ(ierr); 1179 ierr = PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf);CHKERRQ(ierr); 1180 ierr = PetscViewerASCIIPrintf(viewer," Asb = \n");CHKERRQ(ierr); 1181 for (i=0; i<tab->s; i++) { 1182 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s]);CHKERRQ(ierr); 1183 ierr = PetscViewerASCIIPrintf(viewer," %s\n",sbuf);CHKERRQ(ierr); 1184 } 1185 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb);CHKERRQ(ierr); 1186 ierr = PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf);CHKERRQ(ierr); 1187 1188 if (tab->np == 3) { 1189 char mbuf[512]; 1190 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb);CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf);CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," Amb = \n");CHKERRQ(ierr); 1193 for (i=0; i<tab->s; i++) { 1194 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s]);CHKERRQ(ierr); 1195 ierr = PetscViewerASCIIPrintf(viewer," %s\n",mbuf);CHKERRQ(ierr); 1196 } 1197 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb);CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf);CHKERRQ(ierr); 1199 } 1200 } 1201 PetscFunctionReturn(0); 1202 } 1203 1204 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1205 { 1206 PetscErrorCode ierr; 1207 TSAdapt adapt; 1208 1209 PetscFunctionBegin; 1210 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1211 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1212 PetscFunctionReturn(0); 1213 } 1214 1215 /*@C 1216 TSMPRKSetType - Set the type of MPRK scheme 1217 1218 Not collective 1219 1220 Input Parameter: 1221 + ts - timestepping context 1222 - mprktype - type of MPRK-scheme 1223 1224 Options Database: 1225 . -ts_mprk_type - <pm2,p2,p3> 1226 1227 Level: intermediate 1228 1229 .seealso: TSMPRKGetType(), TSMPRK, TSMPRKType 1230 @*/ 1231 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1232 { 1233 PetscErrorCode ierr; 1234 1235 PetscFunctionBegin; 1236 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1237 PetscValidCharPointer(mprktype,2); 1238 ierr = PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype));CHKERRQ(ierr); 1239 PetscFunctionReturn(0); 1240 } 1241 1242 /*@C 1243 TSMPRKGetType - Get the type of MPRK scheme 1244 1245 Not collective 1246 1247 Input Parameter: 1248 . ts - timestepping context 1249 1250 Output Parameter: 1251 . mprktype - type of MPRK-scheme 1252 1253 Level: intermediate 1254 1255 .seealso: TSMPRKGetType() 1256 @*/ 1257 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1258 { 1259 PetscErrorCode ierr; 1260 1261 PetscFunctionBegin; 1262 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1263 ierr = PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype));CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1268 { 1269 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1270 1271 PetscFunctionBegin; 1272 *mprktype = mprk->tableau->name; 1273 PetscFunctionReturn(0); 1274 } 1275 1276 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1277 { 1278 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1279 PetscBool match; 1280 MPRKTableauLink link; 1281 PetscErrorCode ierr; 1282 1283 PetscFunctionBegin; 1284 if (mprk->tableau) { 1285 ierr = PetscStrcmp(mprk->tableau->name,mprktype,&match);CHKERRQ(ierr); 1286 if (match) PetscFunctionReturn(0); 1287 } 1288 for (link = MPRKTableauList; link; link=link->next) { 1289 ierr = PetscStrcmp(link->tab.name,mprktype,&match);CHKERRQ(ierr); 1290 if (match) { 1291 if (ts->setupcalled) {ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr);} 1292 mprk->tableau = &link->tab; 1293 if (ts->setupcalled) {ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr);} 1294 PetscFunctionReturn(0); 1295 } 1296 } 1297 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1302 { 1303 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1304 1305 PetscFunctionBegin; 1306 *ns = mprk->tableau->s; 1307 if (Y) *Y = mprk->Y; 1308 PetscFunctionReturn(0); 1309 } 1310 1311 static PetscErrorCode TSDestroy_MPRK(TS ts) 1312 { 1313 PetscErrorCode ierr; 1314 1315 PetscFunctionBegin; 1316 ierr = TSReset_MPRK(ts);CHKERRQ(ierr); 1317 if (ts->dm) { 1318 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1319 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1320 } 1321 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1322 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL);CHKERRQ(ierr); 1323 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL);CHKERRQ(ierr); 1324 PetscFunctionReturn(0); 1325 } 1326 1327 /*MC 1328 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes 1329 1330 The user should provide the right hand side of the equation 1331 using TSSetRHSFunction(). 1332 1333 Notes: 1334 The default is TSMPRKPM2, it can be changed with TSMPRKSetType() or -ts_mprk_type 1335 1336 Level: beginner 1337 1338 .seealso: TSCreate(), TS, TSSetType(), TSMPRKSetType(), TSMPRKGetType(), TSMPRKType, TSMPRKRegister(), TSMPRKSetMultirateType() 1339 TSMPRKM2, TSMPRKM3, TSMPRKRFSMR3, TSMPRKRFSMR2 1340 1341 M*/ 1342 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1343 { 1344 TS_MPRK *mprk; 1345 PetscErrorCode ierr; 1346 1347 PetscFunctionBegin; 1348 ierr = TSMPRKInitializePackage();CHKERRQ(ierr); 1349 1350 ts->ops->reset = TSReset_MPRK; 1351 ts->ops->destroy = TSDestroy_MPRK; 1352 ts->ops->view = TSView_MPRK; 1353 ts->ops->load = TSLoad_MPRK; 1354 ts->ops->setup = TSSetUp_MPRK; 1355 ts->ops->step = TSStep_MPRK; 1356 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1357 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1358 ts->ops->getstages = TSGetStages_MPRK; 1359 1360 ierr = PetscNewLog(ts,&mprk);CHKERRQ(ierr); 1361 ts->data = (void*)mprk; 1362 1363 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK);CHKERRQ(ierr); 1364 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK);CHKERRQ(ierr); 1365 1366 ierr = TSMPRKSetType(ts,TSMPRKDefault);CHKERRQ(ierr); 1367 PetscFunctionReturn(0); 1368 } 1369