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,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 678 PetscInt i,j,computedstages; 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 computedstages = 0; 702 for (j=0; j<i; j++) { 703 if (tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->Asb[i*s+j]; 704 else ws[computedstages++] = h*tab->Asb[i*s+j]; 705 } 706 for (j=0; j<computedstages; j++) { 707 ierr = VecGetSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j]);CHKERRQ(ierr); 708 } 709 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 710 ierr = VecMAXPY(Yslow,computedstages,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 711 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 712 for (j=0; j<computedstages; j++) { 713 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j]);CHKERRQ(ierr); 714 } 715 } 716 717 /* fast region */ 718 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 719 for (j=0; j<i; j++) { 720 ierr = VecGetSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j]);CHKERRQ(ierr); 721 } 722 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 723 ierr = VecMAXPY(Yfast,i,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 724 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 725 for (j=0; j<i; j++) { 726 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j]);CHKERRQ(ierr); 727 } 728 if (tab->np == 3) { 729 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 730 Vec Ymedium,Ymediumbuffer; 731 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 732 /* medium region */ 733 if (mprk->is_medium) { 734 computedstages = 0; 735 for (j=0; j<i; j++) { 736 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->Amb[i*s+j]; 737 else wm[computedstages++] = h*tab->Amb[i*s+j]; 738 } 739 for (j=0; j<computedstages; j++) { 740 ierr = VecGetSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 741 } 742 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 743 ierr = VecMAXPY(Ymedium,computedstages,wm,mprk->YdotRHS_medium);CHKERRQ(ierr); 744 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 745 for (j=0; j<computedstages; j++) { 746 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 747 } 748 } 749 /* medium buffer region */ 750 for (j=0; j<i; j++) wmb[j] = h*tab->Amb[i*s+j]; 751 for (j=0; j<i; j++) { 752 ierr = VecGetSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 753 } 754 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 755 ierr = VecMAXPY(Ymediumbuffer,i,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 756 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 757 for (j=0; j<i; j++) { 758 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 759 } 760 } 761 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 762 /* compute the stage RHS by fast and slow tableau respectively */ 763 ierr = TSComputeRHSFunction(ts,t+h*csb[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 764 } 765 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 766 ts->ptime += ts->time_step; 767 ts->time_step = next_time_step; 768 PetscFunctionReturn(0); 769 } 770 771 /* 772 This if for partitioned RHS MPRK 773 The step completion formula is 774 775 x1 = x0 + h b^T YdotRHS 776 777 */ 778 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 779 { 780 TS_MPRK *mprk = (TS_MPRK*)ts->data; 781 MPRKTableau tab = mprk->tableau; 782 Vec Xslow,Xfast,Xslowbuffer; /* subvectors for slow and fast componets in X respectively */ 783 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 784 PetscReal h = ts->time_step; 785 PetscInt s = tab->s,j,computedstages; 786 PetscErrorCode ierr; 787 788 PetscFunctionBegin; 789 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 790 791 /* slow region */ 792 if (mprk->is_slow) { 793 computedstages = 0; 794 for (j=0; j<s; j++) { 795 if (tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->bsb[j]; 796 else ws[computedstages++] = h*tab->bsb[j]; 797 } 798 ierr = VecGetSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 799 ierr = VecMAXPY(Xslow,computedstages,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 800 ierr = VecRestoreSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 801 } 802 803 /* slow buffer region */ 804 for (j=0; j<s; j++) wsb[j] = h*tab->bsb[j]; 805 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 806 ierr = VecMAXPY(Xslowbuffer,s,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 807 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 808 809 if (tab->np == 3) { 810 Vec Xmedium,Xmediumbuffer; 811 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 812 /* medium region */ 813 if (mprk->is_medium) { 814 computedstages = 0; 815 for (j=0; j<s; j++) { 816 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bmb[j]; 817 else wm[computedstages++] = h*tab->bmb[j]; 818 } 819 ierr = VecGetSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 820 ierr = VecMAXPY(Xmedium,computedstages,wm,mprk->YdotRHS_medium);CHKERRQ(ierr); 821 ierr = VecRestoreSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 822 } 823 /* medium buffer region */ 824 for (j=0; j<s; j++) wmb[j] = h*tab->bmb[j]; 825 ierr = VecGetSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 826 ierr = VecMAXPY(Xmediumbuffer,s,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 827 ierr = VecRestoreSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 828 } 829 /* fast region */ 830 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 831 ierr = VecGetSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 832 ierr = VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 833 ierr = VecRestoreSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 838 { 839 TS_MPRK *mprk = (TS_MPRK*)ts->data; 840 MPRKTableau tab = mprk->tableau; 841 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 842 Vec Yslow,Yslowbuffer,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 843 PetscInt s = tab->s; 844 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 845 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 846 PetscInt i,j,computedstages; 847 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 848 PetscErrorCode ierr; 849 850 PetscFunctionBegin; 851 for (i=0; i<s; i++) { 852 mprk->stage_time = t + h*cf[i]; 853 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 854 /* calculate the stage value for fast and slow components respectively */ 855 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 856 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 857 858 /* slow buffer region */ 859 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 860 ierr = VecMAXPY(Yslowbuffer,i,wsb,YdotRHS_slowbuffer);CHKERRQ(ierr); 861 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 862 863 /* slow region */ 864 if (mprk->is_slow) { 865 if (tab->rsb[i]) { /* repeat previous stage */ 866 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 867 ierr = VecISCopy(Y[tab->rsb[i]-1],mprk->is_slow,SCATTER_REVERSE,Yslow);CHKERRQ(ierr); 868 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 869 } else { 870 computedstages = 0; 871 for (j=0; j<i; j++) { 872 if (tab->rsb[j]) ws[tab->rsb[j]-1] += wsb[j]; 873 else ws[computedstages++] = wsb[j]; 874 } 875 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 876 ierr = VecMAXPY(Yslow,computedstages,ws,YdotRHS_slow);CHKERRQ(ierr); 877 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 878 /* only depends on the slow buffer region */ 879 ierr = TSComputeRHSFunction(mprk->subts_slow,t+h*csb[i],Y[i],YdotRHS_slow[computedstages]);CHKERRQ(ierr); 880 } 881 } 882 883 /* fast region */ 884 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 885 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 886 ierr = VecMAXPY(Yfast,i,wf,YdotRHS_fast);CHKERRQ(ierr); 887 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 888 889 if (tab->np == 3) { 890 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 891 Vec Ymedium,Ymediumbuffer; 892 const PetscReal *Amb = tab->Amb,*cmb = tab->cmb; 893 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 894 895 for (j=0; j<i; j++) wmb[j] = h*Amb[i*s+j]; 896 /* medium buffer region */ 897 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 898 ierr = VecMAXPY(Ymediumbuffer,i,wmb,YdotRHS_mediumbuffer);CHKERRQ(ierr); 899 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 900 /* medium region */ 901 if (mprk->is_medium) { 902 if (tab->rmb[i]) { /* repeat previous stage */ 903 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 904 ierr = VecISCopy(Y[tab->rmb[i]-1],mprk->is_medium,SCATTER_REVERSE,Ymedium);CHKERRQ(ierr); 905 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 906 } else { 907 computedstages = 0; 908 for (j=0; j<i; j++) { 909 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wmb[j]; 910 else wm[computedstages++] = wmb[j]; 911 912 } 913 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 914 ierr = VecMAXPY(Ymedium,computedstages,wm,YdotRHS_medium);CHKERRQ(ierr); 915 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 916 /* only depends on the medium buffer region */ 917 ierr = TSComputeRHSFunction(mprk->subts_medium,t+h*cmb[i],Y[i],YdotRHS_medium[computedstages]);CHKERRQ(ierr); 918 } 919 } 920 /* must be computed after fast region and slow region are updated in Y */ 921 ierr = TSComputeRHSFunction(mprk->subts_mediumbuffer,t+h*cmb[i],Y[i],YdotRHS_mediumbuffer[i]);CHKERRQ(ierr); 922 } 923 /* must be computed after all regions are updated in Y */ 924 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[i]);CHKERRQ(ierr); 925 ierr = TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i]); 926 } 927 928 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 929 ts->ptime += ts->time_step; 930 ts->time_step = next_time_step; 931 PetscFunctionReturn(0); 932 } 933 934 static PetscErrorCode TSMPRKTableauReset(TS ts) 935 { 936 TS_MPRK *mprk = (TS_MPRK*)ts->data; 937 MPRKTableau tab = mprk->tableau; 938 PetscErrorCode ierr; 939 940 PetscFunctionBegin; 941 if (!tab) PetscFunctionReturn(0); 942 ierr = PetscFree(mprk->work_fast);CHKERRQ(ierr); 943 ierr = PetscFree(mprk->work_slow);CHKERRQ(ierr); 944 ierr = PetscFree(mprk->work_slowbuffer);CHKERRQ(ierr); 945 ierr = PetscFree(mprk->work_medium);CHKERRQ(ierr); 946 ierr = PetscFree(mprk->work_mediumbuffer);CHKERRQ(ierr); 947 ierr = VecDestroyVecs(tab->s,&mprk->Y);CHKERRQ(ierr); 948 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 949 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 950 if (mprk->is_slow) { 951 ierr = PetscFree(mprk->YdotRHS_slow);CHKERRQ(ierr); 952 } 953 ierr = PetscFree(mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 954 if (tab->np == 3) { 955 if (mprk->is_medium) { 956 ierr = PetscFree(mprk->YdotRHS_medium);CHKERRQ(ierr); 957 } 958 ierr = PetscFree(mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 959 } 960 ierr = PetscFree(mprk->YdotRHS_fast);CHKERRQ(ierr); 961 962 } else { 963 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 964 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 965 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 966 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 967 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 968 } 969 PetscFunctionReturn(0); 970 } 971 972 static PetscErrorCode TSReset_MPRK(TS ts) 973 { 974 PetscErrorCode ierr; 975 976 PetscFunctionBegin; 977 ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr); 978 PetscFunctionReturn(0); 979 } 980 981 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 982 { 983 PetscFunctionBegin; 984 PetscFunctionReturn(0); 985 } 986 987 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 988 { 989 PetscFunctionBegin; 990 PetscFunctionReturn(0); 991 } 992 993 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 994 { 995 PetscFunctionBegin; 996 PetscFunctionReturn(0); 997 } 998 999 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1000 { 1001 PetscFunctionBegin; 1002 PetscFunctionReturn(0); 1003 } 1004 1005 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 1006 { 1007 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1008 MPRKTableau tab = mprk->tableau; 1009 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y);CHKERRQ(ierr); 1014 if (mprk->is_slow) { 1015 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 1016 } 1017 ierr = PetscMalloc1(tab->s,&mprk->work_slowbuffer);CHKERRQ(ierr); 1018 if (tab->np == 3) { 1019 if (mprk->is_medium) { 1020 ierr = PetscMalloc1(tab->s,&mprk->work_medium);CHKERRQ(ierr); 1021 } 1022 ierr = PetscMalloc1(tab->s,&mprk->work_mediumbuffer);CHKERRQ(ierr); 1023 } 1024 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 1025 1026 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 1027 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 1028 if (mprk->is_slow) { 1029 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1030 } 1031 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1032 if (tab->np == 3) { 1033 if (mprk->is_medium) { 1034 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1035 } 1036 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1037 } 1038 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1039 } else { 1040 if (mprk->is_slow) { 1041 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1042 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1043 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1044 } 1045 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1046 ierr = VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1047 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1048 if (tab->np == 3) { 1049 if (mprk->is_medium) { 1050 ierr = VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1051 ierr = VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1052 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1053 } 1054 ierr = VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1055 ierr = VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1056 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1057 } 1058 ierr = VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1059 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1060 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1061 } 1062 PetscFunctionReturn(0); 1063 } 1064 1065 static PetscErrorCode TSSetUp_MPRK(TS ts) 1066 { 1067 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1068 MPRKTableau tab = mprk->tableau; 1069 DM dm; 1070 PetscErrorCode ierr; 1071 1072 PetscFunctionBegin; 1073 ierr = TSRHSSplitGetIS(ts,"slow",&mprk->is_slow);CHKERRQ(ierr); 1074 ierr = TSRHSSplitGetIS(ts,"fast",&mprk->is_fast);CHKERRQ(ierr); 1075 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); 1076 1077 if (tab->np == 3) { 1078 ierr = TSRHSSplitGetIS(ts,"medium",&mprk->is_medium);CHKERRQ(ierr); 1079 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); 1080 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 1081 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1082 mprk->is_mediumbuffer = mprk->is_medium; 1083 mprk->is_medium = NULL; 1084 } 1085 } 1086 1087 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1088 ierr = TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer);CHKERRQ(ierr); 1089 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1090 mprk->is_slowbuffer = mprk->is_slow; 1091 mprk->is_slow = NULL; 1092 } 1093 /* 1094 if (!mprk->is_medium) { 1095 mprk->is_medium = mprk->is_fast; 1096 mprk->is_fast = NULL; 1097 } else { 1098 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 1099 } 1100 */ 1101 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1102 ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr); 1103 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1104 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1105 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1106 ierr = PetscTryMethod(ts,"TSMPRKSetSplits_C",(TS),(ts));CHKERRQ(ierr); 1107 PetscFunctionReturn(0); 1108 } 1109 1110 /* construct a database to chose nonsplit RHS mutirate mprk method or split RHS MPRK method */ 1111 const char *const TSMPRKMultirateTypes[] = {"NONSPLIT","SPLIT","TSMPRKMultirateType","TSMPRK",0}; 1112 1113 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1114 { 1115 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1116 PetscErrorCode ierr; 1117 1118 PetscFunctionBegin; 1119 ierr = PetscOptionsHead(PetscOptionsObject,"PRK ODE solver options");CHKERRQ(ierr); 1120 { 1121 MPRKTableauLink link; 1122 PetscInt count,choice; 1123 PetscBool flg; 1124 const char **namelist; 1125 PetscInt mprkmtype = 0; 1126 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1127 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1128 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1129 ierr = PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1130 if (flg) {ierr = TSMPRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1131 ierr = PetscFree(namelist);CHKERRQ(ierr); 1132 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); 1133 if (flg) {ierr = TSMPRKSetMultirateType(ts,mprkmtype);CHKERRQ(ierr);} 1134 } 1135 ierr = PetscOptionsTail();CHKERRQ(ierr); 1136 PetscFunctionReturn(0); 1137 } 1138 1139 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1140 { 1141 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1142 PetscBool iascii; 1143 PetscErrorCode ierr; 1144 1145 PetscFunctionBegin; 1146 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1147 if (iascii) { 1148 MPRKTableau tab = mprk->tableau; 1149 TSMPRKType mprktype; 1150 char fbuf[512]; 1151 char sbuf[512]; 1152 PetscInt i; 1153 ierr = TSMPRKGetType(ts,&mprktype);CHKERRQ(ierr); 1154 ierr = PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype);CHKERRQ(ierr); 1155 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1156 1157 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf);CHKERRQ(ierr); 1158 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf);CHKERRQ(ierr); 1159 ierr = PetscViewerASCIIPrintf(viewer," Af = \n");CHKERRQ(ierr); 1160 for (i=0; i<tab->s; i++) { 1161 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s]);CHKERRQ(ierr); 1162 ierr = PetscViewerASCIIPrintf(viewer," %s\n",fbuf);CHKERRQ(ierr); 1163 } 1164 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf);CHKERRQ(ierr); 1165 ierr = PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf);CHKERRQ(ierr); 1166 1167 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb);CHKERRQ(ierr); 1168 ierr = PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf);CHKERRQ(ierr); 1169 ierr = PetscViewerASCIIPrintf(viewer," Asb = \n");CHKERRQ(ierr); 1170 for (i=0; i<tab->s; i++) { 1171 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s]);CHKERRQ(ierr); 1172 ierr = PetscViewerASCIIPrintf(viewer," %s\n",sbuf);CHKERRQ(ierr); 1173 } 1174 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb);CHKERRQ(ierr); 1175 ierr = PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf);CHKERRQ(ierr); 1176 1177 if (tab->np == 3) { 1178 char mbuf[512]; 1179 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb);CHKERRQ(ierr); 1180 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf);CHKERRQ(ierr); 1181 ierr = PetscViewerASCIIPrintf(viewer," Amb = \n");CHKERRQ(ierr); 1182 for (i=0; i<tab->s; i++) { 1183 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s]);CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer," %s\n",mbuf);CHKERRQ(ierr); 1185 } 1186 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb);CHKERRQ(ierr); 1187 ierr = PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf);CHKERRQ(ierr); 1188 } 1189 } 1190 PetscFunctionReturn(0); 1191 } 1192 1193 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1194 { 1195 PetscErrorCode ierr; 1196 TSAdapt adapt; 1197 1198 PetscFunctionBegin; 1199 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1200 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1201 PetscFunctionReturn(0); 1202 } 1203 1204 /*@C 1205 TSMPRKSetType - Set the type of MPRK scheme 1206 1207 Logically collective 1208 1209 Input Parameter: 1210 + ts - timestepping context 1211 - mprktype - type of MPRK-scheme 1212 1213 Options Database: 1214 . -ts_mprk_type - <pm2,p2,p3> 1215 1216 Level: intermediate 1217 1218 .seealso: TSMPRKGetType(), TSMPRK, TSMPRKType 1219 @*/ 1220 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1221 { 1222 PetscErrorCode ierr; 1223 1224 PetscFunctionBegin; 1225 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1226 PetscValidCharPointer(mprktype,2); 1227 ierr = PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype));CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 /*@C 1232 TSMPRKGetType - Get the type of MPRK scheme 1233 1234 Logically collective 1235 1236 Input Parameter: 1237 . ts - timestepping context 1238 1239 Output Parameter: 1240 . mprktype - type of MPRK-scheme 1241 1242 Level: intermediate 1243 1244 .seealso: TSMPRKGetType() 1245 @*/ 1246 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1247 { 1248 PetscErrorCode ierr; 1249 1250 PetscFunctionBegin; 1251 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1252 ierr = PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype));CHKERRQ(ierr); 1253 PetscFunctionReturn(0); 1254 } 1255 1256 /*@C 1257 TSMPRKSetMultirateType - Set the type of MPRK multirate scheme 1258 1259 Logically collective 1260 1261 Input Parameter: 1262 + ts - timestepping context 1263 - mprkmtype - type of the multirate configuration 1264 1265 Options Database: 1266 . -ts_mprk_multirate_type - <nonsplit,split> 1267 1268 Level: intermediate 1269 @*/ 1270 PetscErrorCode TSMPRKSetMultirateType(TS ts, TSMPRKMultirateType mprkmtype) 1271 { 1272 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1273 PetscErrorCode ierr; 1274 1275 PetscFunctionBegin; 1276 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1277 switch(mprkmtype){ 1278 case TSMPRKNONSPLIT: 1279 ts->ops->step = TSStep_MPRK; 1280 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1281 break; 1282 case TSMPRKSPLIT: 1283 ts->ops->step = TSStep_MPRKSPLIT; 1284 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1285 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetSplits_C",TSMPRKSetSplits);CHKERRQ(ierr); 1286 break; 1287 default : 1288 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",mprkmtype); 1289 } 1290 mprk->mprkmtype = mprkmtype; 1291 PetscFunctionReturn(0); 1292 } 1293 1294 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1295 { 1296 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1297 1298 PetscFunctionBegin; 1299 *mprktype = mprk->tableau->name; 1300 PetscFunctionReturn(0); 1301 } 1302 1303 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1304 { 1305 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1306 PetscBool match; 1307 MPRKTableauLink link; 1308 PetscErrorCode ierr; 1309 1310 PetscFunctionBegin; 1311 if (mprk->tableau) { 1312 ierr = PetscStrcmp(mprk->tableau->name,mprktype,&match);CHKERRQ(ierr); 1313 if (match) PetscFunctionReturn(0); 1314 } 1315 for (link = MPRKTableauList; link; link=link->next) { 1316 ierr = PetscStrcmp(link->tab.name,mprktype,&match);CHKERRQ(ierr); 1317 if (match) { 1318 if (ts->setupcalled) {ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr);} 1319 mprk->tableau = &link->tab; 1320 if (ts->setupcalled) {ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr);} 1321 PetscFunctionReturn(0); 1322 } 1323 } 1324 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1329 { 1330 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1331 1332 PetscFunctionBegin; 1333 *ns = mprk->tableau->s; 1334 if (Y) *Y = mprk->Y; 1335 PetscFunctionReturn(0); 1336 } 1337 1338 static PetscErrorCode TSDestroy_MPRK(TS ts) 1339 { 1340 PetscErrorCode ierr; 1341 1342 PetscFunctionBegin; 1343 ierr = TSReset_MPRK(ts);CHKERRQ(ierr); 1344 if (ts->dm) { 1345 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1346 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1347 } 1348 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1349 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL);CHKERRQ(ierr); 1350 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL);CHKERRQ(ierr); 1351 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",NULL);CHKERRQ(ierr); 1352 PetscFunctionReturn(0); 1353 } 1354 1355 /*MC 1356 TSMPRK - ODE solver using Partitioned Runge-Kutta schemes 1357 1358 The user should provide the right hand side of the equation 1359 using TSSetRHSFunction(). 1360 1361 Notes: 1362 The default is TSMPRKPM2, it can be changed with TSRKSetType() or -ts_mprk_type 1363 1364 Level: beginner 1365 1366 .seealso: TSCreate(), TS, TSSetType(), TSMPRKSetType(), TSMPRKGetType(), TSMPRKType, TSMPRKRegister(), TSMPRKSetMultirateType() 1367 TSMPRKM2, TSMPRKM3, TSMPRKRFSMR3, TSMPRKRFSMR2 1368 1369 M*/ 1370 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1371 { 1372 TS_MPRK *mprk; 1373 PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 ierr = TSMPRKInitializePackage();CHKERRQ(ierr); 1377 1378 ts->ops->reset = TSReset_MPRK; 1379 ts->ops->destroy = TSDestroy_MPRK; 1380 ts->ops->view = TSView_MPRK; 1381 ts->ops->load = TSLoad_MPRK; 1382 ts->ops->setup = TSSetUp_MPRK; 1383 ts->ops->step = TSStep_MPRK; 1384 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1385 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1386 ts->ops->getstages = TSGetStages_MPRK; 1387 1388 ierr = PetscNewLog(ts,&mprk);CHKERRQ(ierr); 1389 ts->data = (void*)mprk; 1390 1391 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK);CHKERRQ(ierr); 1392 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK);CHKERRQ(ierr); 1393 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",TSMPRKSetMultirateType);CHKERRQ(ierr); 1394 1395 ierr = TSMPRKSetType(ts,TSMPRKDefault);CHKERRQ(ierr); 1396 PetscFunctionReturn(0); 1397 } 1398