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