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