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