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