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