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