1 /* 2 Code for time stepping with the Partitioned Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) for nonsplit RHS multi-rate RK, 7 user should give the indexes for both slow and fast components; 8 2) The general system is written as 9 Usdot = Fs(t,Us,Uf) 10 Ufdot = Ff(t,Us,Uf) for split RHS multi-rate RK, 11 user should partioned RHS by themselves and also provide the indexes for both slow and fast components. 12 3) To correct The confusing terminology in the paper, we use 'slow method', 'slow buffer method' and 'fast method' to denote the methods applied to 'slow region', 'slow buffer region' and 'fast region' respectively. The 'slow method' in the original paper actually means the 'slow buffer method'. 13 4) Why does the buffer region have to be inside the slow region? The buffer region is treated with a slow method essentially. Applying the slow method to a region with a fast characteristic time scale is apparently not a good choice. 14 15 Reference: 16 Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007 17 18 */ 19 20 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 21 #include <petscdm.h> 22 23 static TSMPRKType TSMPRKDefault = TSMPRK2A22; 24 static PetscBool TSMPRKRegisterAllCalled; 25 static PetscBool TSMPRKPackageInitialized; 26 27 typedef struct _MPRKTableau *MPRKTableau; 28 struct _MPRKTableau { 29 char *name; 30 PetscInt order; /* Classical approximation order of the method i */ 31 PetscInt sbase; /* Number of stages in the base method*/ 32 PetscInt s; /* Number of stages */ 33 PetscInt np; /* Number of partitions */ 34 PetscReal *Af,*bf,*cf; /* Tableau for fast components */ 35 PetscReal *Amb,*bmb,*cmb; /* Tableau for medium components */ 36 PetscInt *rmb; /* Array of flags for repeated stages in medium method */ 37 PetscReal *Asb,*bsb,*csb; /* Tableau for slow components */ 38 PetscInt *rsb; /* Array of flags for repeated staged in slow method*/ 39 }; 40 typedef struct _MPRKTableauLink *MPRKTableauLink; 41 struct _MPRKTableauLink { 42 struct _MPRKTableau tab; 43 MPRKTableauLink next; 44 }; 45 static MPRKTableauLink MPRKTableauList; 46 47 typedef struct { 48 MPRKTableau tableau; 49 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,2,2,1,&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,2,2,2,&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,2,3,1,&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,2,3,3,&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,1,1,&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,5,2,1,&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 in the base methods 467 . ratio1 - stepsize ratio at 1st level (e.g. slow/medium) 468 . ratio2 - stepsize ratio at 2nd level (e.g. medium/fast) 469 . Af - stage coefficients for fast components(dimension s*s, row-major) 470 . bf - step completion table for fast components(dimension s) 471 . cf - abscissa for fast components(dimension s) 472 . As - stage coefficients for slow components(dimension s*s, row-major) 473 . bs - step completion table for slow components(dimension s) 474 - cs - abscissa for slow components(dimension s) 475 476 Notes: 477 Several MPRK methods are provided, this function is only needed to create new methods. 478 479 Level: advanced 480 481 .keywords: TS, register 482 483 .seealso: TSMPRK 484 @*/ 485 PetscErrorCode TSMPRKRegister(TSMPRKType name,PetscInt order, 486 PetscInt sbase,PetscInt ratio1,PetscInt ratio2, 487 const PetscReal Asb[],const PetscReal bsb[],const PetscReal csb[],const PetscInt rsb[], 488 const PetscReal Amb[],const PetscReal bmb[],const PetscReal cmb[],const PetscInt rmb[], 489 const PetscReal Af[],const PetscReal bf[],const PetscReal cf[]) 490 { 491 MPRKTableauLink link; 492 MPRKTableau t; 493 PetscInt s,i,j; 494 PetscErrorCode ierr; 495 496 PetscFunctionBegin; 497 PetscValidCharPointer(name,1); 498 PetscValidRealPointer(Asb,4); 499 if (bsb) PetscValidRealPointer(bsb,7); 500 if (csb) PetscValidRealPointer(csb,8); 501 if (rsb) PetscValidRealPointer(rsb,9); 502 if (Amb) PetscValidRealPointer(Amb,10); 503 if (bmb) PetscValidRealPointer(bmb,11); 504 if (cmb) PetscValidRealPointer(cmb,12); 505 if (rmb) PetscValidRealPointer(rmb,13); 506 PetscValidRealPointer(Af,14); 507 if (bf) PetscValidRealPointer(bf,15); 508 if (cf) PetscValidRealPointer(cf,16); 509 510 ierr = PetscNew(&link);CHKERRQ(ierr); 511 t = &link->tab; 512 513 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 514 s = sbase*ratio1*ratio2; /* this is the dimension of the matrices below */ 515 t->order = order; 516 t->sbase = sbase; 517 t->s = s; 518 t->np = 2; 519 520 ierr = PetscMalloc3(s*s,&t->Af,s,&t->bf,s,&t->cf);CHKERRQ(ierr); 521 ierr = PetscMemcpy(t->Af,Af,s*s*sizeof(Af[0]));CHKERRQ(ierr); 522 if (bf) { 523 ierr = PetscMemcpy(t->bf,bf,s*sizeof(bf[0]));CHKERRQ(ierr); 524 } else 525 for (i=0; i<s; i++) t->bf[i] = Af[(s-1)*s+i]; 526 if (cf) { 527 ierr = PetscMemcpy(t->cf,cf,s*sizeof(cf[0]));CHKERRQ(ierr); 528 } else { 529 for (i=0; i<s; i++) 530 for (j=0,t->cf[i]=0; j<s; j++) 531 t->cf[i] += Af[i*s+j]; 532 } 533 534 if (Amb) { 535 t->np = 3; 536 ierr = PetscMalloc3(s*s,&t->Amb,s,&t->bmb,s,&t->cmb);CHKERRQ(ierr); 537 ierr = PetscCalloc1(s,&t->rmb);CHKERRQ(ierr); 538 ierr = PetscMemcpy(t->Amb,Amb,s*s*sizeof(Amb[0]));CHKERRQ(ierr); 539 if (bmb) { 540 ierr = PetscMemcpy(t->bmb,bmb,s*sizeof(bmb[0]));CHKERRQ(ierr); 541 } else { 542 for (i=0; i<s; i++) t->bmb[i] = Amb[(s-1)*s+i]; 543 } 544 if (cmb) { 545 ierr = PetscMemcpy(t->cmb,cmb,s*sizeof(cmb[0]));CHKERRQ(ierr); 546 } else { 547 for (i=0; i<s; i++) 548 for (j=0,t->cmb[i]=0; j<s; j++) 549 t->cmb[i] += Amb[i*s+j]; 550 } 551 if (rmb) { 552 ierr = PetscMemcpy(t->rmb,rmb,s*sizeof(rmb[0]));CHKERRQ(ierr); 553 } 554 } 555 556 ierr = PetscMalloc3(s*s,&t->Asb,s,&t->bsb,s,&t->csb);CHKERRQ(ierr); 557 ierr = PetscMemcpy(t->Asb,Asb,s*s*sizeof(Asb[0]));CHKERRQ(ierr); 558 ierr = PetscCalloc1(s,&t->rsb);CHKERRQ(ierr); 559 if (bsb) { 560 ierr = PetscMemcpy(t->bsb,bsb,s*sizeof(bsb[0]));CHKERRQ(ierr); 561 } else 562 for (i=0; i<s; i++) t->bsb[i] = Asb[(s-1)*s+i]; 563 if (csb) { 564 ierr = PetscMemcpy(t->csb,csb,s*sizeof(csb[0]));CHKERRQ(ierr); 565 } else { 566 for (i=0; i<s; i++) 567 for (j=0,t->csb[i]=0; j<s; j++) 568 t->csb[i] += Asb[i*s+j]; 569 } 570 if (rsb) { 571 ierr = PetscMemcpy(t->rsb,rsb,s*sizeof(rsb[0]));CHKERRQ(ierr); 572 } 573 link->next = MPRKTableauList; 574 MPRKTableauList = link; 575 PetscFunctionReturn(0); 576 } 577 578 static PetscErrorCode TSMPRKSetSplits(TS ts) 579 { 580 TS_MPRK *mprk = (TS_MPRK*)ts->data; 581 MPRKTableau tab = mprk->tableau; 582 DM dm,subdm,newdm; 583 PetscErrorCode ierr; 584 585 PetscFunctionBegin; 586 ierr = TSRHSSplitGetSubTS(ts,"slow",&mprk->subts_slow);CHKERRQ(ierr); 587 ierr = TSRHSSplitGetSubTS(ts,"fast",&mprk->subts_fast);CHKERRQ(ierr); 588 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"); 589 590 /* Only copy the DM */ 591 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 592 593 ierr = TSRHSSplitGetSubTS(ts,"slowbuffer",&mprk->subts_slowbuffer);CHKERRQ(ierr); 594 if (!mprk->subts_slowbuffer) { 595 mprk->subts_slowbuffer = mprk->subts_slow; 596 mprk->subts_slow = NULL; 597 } 598 if (mprk->subts_slow) { 599 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 600 ierr = TSGetDM(mprk->subts_slow,&subdm);CHKERRQ(ierr); 601 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 602 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 603 ierr = TSSetDM(mprk->subts_slow,newdm);CHKERRQ(ierr); 604 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 605 } 606 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 607 ierr = TSGetDM(mprk->subts_slowbuffer,&subdm);CHKERRQ(ierr); 608 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 609 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 610 ierr = TSSetDM(mprk->subts_slowbuffer,newdm);CHKERRQ(ierr); 611 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 612 613 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 614 ierr = TSGetDM(mprk->subts_fast,&subdm);CHKERRQ(ierr); 615 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 616 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 617 ierr = TSSetDM(mprk->subts_fast,newdm);CHKERRQ(ierr); 618 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 619 620 if (tab->np == 3) { 621 ierr = TSRHSSplitGetSubTS(ts,"medium",&mprk->subts_medium);CHKERRQ(ierr); 622 ierr = TSRHSSplitGetSubTS(ts,"mediumbuffer",&mprk->subts_mediumbuffer);CHKERRQ(ierr); 623 if (mprk->subts_medium && !mprk->subts_mediumbuffer) { 624 mprk->subts_mediumbuffer = mprk->subts_medium; 625 mprk->subts_medium = NULL; 626 } 627 if (mprk->subts_medium) { 628 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 629 ierr = TSGetDM(mprk->subts_medium,&subdm);CHKERRQ(ierr); 630 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 631 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 632 ierr = TSSetDM(mprk->subts_medium,newdm);CHKERRQ(ierr); 633 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 634 } 635 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 636 ierr = TSGetDM(mprk->subts_mediumbuffer,&subdm);CHKERRQ(ierr); 637 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 638 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 639 ierr = TSSetDM(mprk->subts_mediumbuffer,newdm);CHKERRQ(ierr); 640 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 641 } 642 PetscFunctionReturn(0); 643 } 644 645 /* 646 This if for nonsplit RHS MPRK 647 The step completion formula is 648 649 x1 = x0 + h b^T YdotRHS 650 651 */ 652 static PetscErrorCode TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool *done) 653 { 654 TS_MPRK *mprk = (TS_MPRK*)ts->data; 655 MPRKTableau tab = mprk->tableau; 656 PetscScalar *wf = mprk->work_fast; 657 PetscReal h = ts->time_step; 658 PetscInt s = tab->s,j; 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 663 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 664 ierr = VecMAXPY(X,s,wf,mprk->YdotRHS);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 static PetscErrorCode TSStep_MPRK(TS ts) 669 { 670 TS_MPRK *mprk = (TS_MPRK*)ts->data; 671 Vec *Y = mprk->Y,*YdotRHS = mprk->YdotRHS,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 672 Vec Yslow,Yslowbuffer,Yfast; 673 MPRKTableau tab = mprk->tableau; 674 const PetscInt s = tab->s; 675 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 676 PetscScalar *wf = mprk->work_fast,*wsb = mprk->work_slowbuffer; 677 PetscInt i,j; 678 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 679 PetscErrorCode ierr; 680 681 PetscFunctionBegin; 682 for (i=0; i<s; i++) { 683 mprk->stage_time = t + h*cf[i]; 684 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 685 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 686 687 /* slow buffer region */ 688 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 689 for (j=0; j<i; j++) { 690 ierr = VecGetSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j]);CHKERRQ(ierr); 691 } 692 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 693 ierr = VecMAXPY(Yslowbuffer,i,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 694 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 695 for (j=0; j<i; j++) { 696 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j]);CHKERRQ(ierr); 697 } 698 /* slow region */ 699 if (mprk->is_slow) { 700 for (j=0; j<i; 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,wsb,mprk->YdotRHS_slow);CHKERRQ(ierr); 705 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 706 for (j=0; j<i; 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<i; 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<i; 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 *wmb = mprk->work_mediumbuffer; 726 727 for (j=0; j<i; j++) wmb[j] = h*tab->Amb[i*s+j]; 728 /* medium region */ 729 if (mprk->is_medium) { 730 for (j=0; j<i; j++) { 731 ierr = VecGetSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 732 } 733 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 734 ierr = VecMAXPY(Ymedium,i,wmb,mprk->YdotRHS_medium);CHKERRQ(ierr); 735 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 736 for (j=0; j<i; j++) { 737 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j]);CHKERRQ(ierr); 738 } 739 } 740 /* medium buffer region */ 741 for (j=0; j<i; j++) { 742 ierr = VecGetSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 743 } 744 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 745 ierr = VecMAXPY(Ymediumbuffer,i,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 746 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 747 for (j=0; j<i; j++) { 748 ierr = VecRestoreSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j]);CHKERRQ(ierr); 749 } 750 } 751 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 752 /* compute the stage RHS by fast and slow tableau respectively */ 753 ierr = TSComputeRHSFunction(ts,t+h*csb[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 754 } 755 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 756 ts->ptime += ts->time_step; 757 ts->time_step = next_time_step; 758 PetscFunctionReturn(0); 759 } 760 761 /* 762 This if for partitioned RHS MPRK 763 The step completion formula is 764 765 x1 = x0 + h b^T YdotRHS 766 767 */ 768 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 769 { 770 TS_MPRK *mprk = (TS_MPRK*)ts->data; 771 MPRKTableau tab = mprk->tableau; 772 Vec Xslow,Xfast,Xslowbuffer; /* subvectors for slow and fast componets in X respectively */ 773 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 774 PetscReal h = ts->time_step; 775 PetscInt s = tab->s,j,computedstages; 776 PetscErrorCode ierr; 777 778 PetscFunctionBegin; 779 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 780 781 /* slow region */ 782 if (mprk->is_slow) { 783 computedstages = 0; 784 for (j=0; j<s; j++) { 785 if (tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->bsb[j]; 786 else ws[computedstages++] = h*tab->bsb[j]; 787 } 788 ierr = VecGetSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 789 ierr = VecMAXPY(Xslow,computedstages,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 790 ierr = VecRestoreSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 791 } 792 793 if (tab->np == 3 && mprk->is_medium) { 794 computedstages = 0; 795 for (j=0; j<s; j++) { 796 if (tab->rmb[j]) wsb[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bsb[j]; 797 else wsb[computedstages++] = h*tab->bsb[j]; 798 } 799 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 800 ierr = VecMAXPY(Xslowbuffer,computedstages,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 801 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 802 } else { 803 /* slow buffer region */ 804 for (j=0; j<s; j++) wsb[j] = h*tab->bsb[j]; 805 ierr = VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 806 ierr = VecMAXPY(Xslowbuffer,s,wsb,mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 807 ierr = VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer);CHKERRQ(ierr); 808 } 809 if (tab->np == 3) { 810 Vec Xmedium,Xmediumbuffer; 811 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 812 /* medium region and slow buffer region */ 813 if (mprk->is_medium) { 814 computedstages = 0; 815 for (j=0; j<s; j++) { 816 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bmb[j]; 817 else wm[computedstages++] = h*tab->bmb[j]; 818 } 819 ierr = VecGetSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 820 ierr = VecMAXPY(Xmedium,computedstages,wm,mprk->YdotRHS_medium);CHKERRQ(ierr); 821 ierr = VecRestoreSubVector(X,mprk->is_medium,&Xmedium);CHKERRQ(ierr); 822 } 823 /* medium buffer region */ 824 for (j=0; j<s; j++) wmb[j] = h*tab->bmb[j]; 825 ierr = VecGetSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 826 ierr = VecMAXPY(Xmediumbuffer,s,wmb,mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 827 ierr = VecRestoreSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer);CHKERRQ(ierr); 828 } 829 /* fast region */ 830 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 831 ierr = VecGetSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 832 ierr = VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 833 ierr = VecRestoreSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 838 { 839 TS_MPRK *mprk = (TS_MPRK*)ts->data; 840 MPRKTableau tab = mprk->tableau; 841 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 842 Vec Yslow,Yslowbuffer,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 843 PetscInt s = tab->s; 844 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 845 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 846 PetscInt i,j,computedstages; 847 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 848 PetscErrorCode ierr; 849 850 PetscFunctionBegin; 851 for (i=0; i<s; i++) { 852 mprk->stage_time = t + h*cf[i]; 853 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 854 /* calculate the stage value for fast and slow components respectively */ 855 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 856 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 857 858 /* slow buffer region */ 859 if (tab->np == 3 && mprk->is_medium) { 860 if (tab->rmb[i]) { 861 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 862 ierr = VecISCopy(Y[tab->rmb[i]-1],mprk->is_slowbuffer,SCATTER_REVERSE,Yslowbuffer);CHKERRQ(ierr); 863 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 864 } else { 865 PetscScalar *wm = mprk->work_medium; 866 computedstages = 0; 867 for (j=0; j<i; j++) { 868 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wsb[j]; 869 else wm[computedstages++] = wsb[j]; 870 } 871 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 872 ierr = VecMAXPY(Yslowbuffer,computedstages,wm,YdotRHS_slowbuffer);CHKERRQ(ierr); 873 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 874 } 875 } else { 876 ierr = VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 877 ierr = VecMAXPY(Yslowbuffer,i,wsb,YdotRHS_slowbuffer);CHKERRQ(ierr); 878 ierr = VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer);CHKERRQ(ierr); 879 } 880 881 /* slow region */ 882 if (mprk->is_slow) { 883 if (tab->rsb[i]) { /* repeat previous stage */ 884 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 885 ierr = VecISCopy(Y[tab->rsb[i]-1],mprk->is_slow,SCATTER_REVERSE,Yslow);CHKERRQ(ierr); 886 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 887 } else { 888 computedstages = 0; 889 for (j=0; j<i; j++) { 890 if (tab->rsb[j]) ws[tab->rsb[j]-1] += wsb[j]; 891 else ws[computedstages++] = wsb[j]; 892 } 893 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 894 ierr = VecMAXPY(Yslow,computedstages,ws,YdotRHS_slow);CHKERRQ(ierr); 895 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 896 /* only depends on the slow buffer region */ 897 ierr = TSComputeRHSFunction(mprk->subts_slow,t+h*csb[i],Y[i],YdotRHS_slow[computedstages]);CHKERRQ(ierr); 898 } 899 } 900 901 /* fast region */ 902 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 903 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 904 ierr = VecMAXPY(Yfast,i,wf,YdotRHS_fast);CHKERRQ(ierr); 905 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 906 907 if (tab->np == 3) { 908 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 909 Vec Ymedium,Ymediumbuffer; 910 const PetscReal *Amb = tab->Amb,*cmb = tab->cmb; 911 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 912 913 for (j=0; j<i; j++) wmb[j] = h*Amb[i*s+j]; 914 /* medium buffer region */ 915 ierr = VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 916 ierr = VecMAXPY(Ymediumbuffer,i,wmb,YdotRHS_mediumbuffer);CHKERRQ(ierr); 917 ierr = VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer);CHKERRQ(ierr); 918 /* medium region */ 919 if (mprk->is_medium) { 920 if (tab->rmb[i]) { /* repeat previous stage */ 921 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 922 ierr = VecISCopy(Y[tab->rmb[i]-1],mprk->is_medium,SCATTER_REVERSE,Ymedium);CHKERRQ(ierr); 923 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 924 } else { 925 computedstages = 0; 926 for (j=0; j<i; j++) { 927 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wmb[j]; 928 else wm[computedstages++] = wmb[j]; 929 930 } 931 ierr = VecGetSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 932 ierr = VecMAXPY(Ymedium,computedstages,wm,YdotRHS_medium);CHKERRQ(ierr); 933 ierr = VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium);CHKERRQ(ierr); 934 /* only depends on the medium buffer region */ 935 ierr = TSComputeRHSFunction(mprk->subts_medium,t+h*cmb[i],Y[i],YdotRHS_medium[computedstages]);CHKERRQ(ierr); 936 /* must be computed after all regions are updated in Y */ 937 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[computedstages]);CHKERRQ(ierr); 938 } 939 } 940 /* must be computed after fast region and slow region are updated in Y */ 941 ierr = TSComputeRHSFunction(mprk->subts_mediumbuffer,t+h*cmb[i],Y[i],YdotRHS_mediumbuffer[i]);CHKERRQ(ierr); 942 } 943 if (!(tab->np == 3 && mprk->is_medium)) { 944 ierr = TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[i]);CHKERRQ(ierr); 945 } 946 ierr = TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i]); 947 } 948 949 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 950 ts->ptime += ts->time_step; 951 ts->time_step = next_time_step; 952 PetscFunctionReturn(0); 953 } 954 955 static PetscErrorCode TSMPRKTableauReset(TS ts) 956 { 957 TS_MPRK *mprk = (TS_MPRK*)ts->data; 958 MPRKTableau tab = mprk->tableau; 959 PetscErrorCode ierr; 960 961 PetscFunctionBegin; 962 if (!tab) PetscFunctionReturn(0); 963 ierr = PetscFree(mprk->work_fast);CHKERRQ(ierr); 964 ierr = PetscFree(mprk->work_slow);CHKERRQ(ierr); 965 ierr = PetscFree(mprk->work_slowbuffer);CHKERRQ(ierr); 966 ierr = PetscFree(mprk->work_medium);CHKERRQ(ierr); 967 ierr = PetscFree(mprk->work_mediumbuffer);CHKERRQ(ierr); 968 ierr = VecDestroyVecs(tab->s,&mprk->Y);CHKERRQ(ierr); 969 if (ts->use_splitrhsfunction) { 970 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 971 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 972 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 973 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 974 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 975 } else { 976 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 977 if (mprk->is_slow) { 978 ierr = PetscFree(mprk->YdotRHS_slow);CHKERRQ(ierr); 979 } 980 ierr = PetscFree(mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 981 if (tab->np == 3) { 982 if (mprk->is_medium) { 983 ierr = PetscFree(mprk->YdotRHS_medium);CHKERRQ(ierr); 984 } 985 ierr = PetscFree(mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 986 } 987 ierr = PetscFree(mprk->YdotRHS_fast);CHKERRQ(ierr); 988 } 989 PetscFunctionReturn(0); 990 } 991 992 static PetscErrorCode TSReset_MPRK(TS ts) 993 { 994 PetscErrorCode ierr; 995 996 PetscFunctionBegin; 997 ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr); 998 PetscFunctionReturn(0); 999 } 1000 1001 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 1002 { 1003 PetscFunctionBegin; 1004 PetscFunctionReturn(0); 1005 } 1006 1007 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1008 { 1009 PetscFunctionBegin; 1010 PetscFunctionReturn(0); 1011 } 1012 1013 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 1014 { 1015 PetscFunctionBegin; 1016 PetscFunctionReturn(0); 1017 } 1018 1019 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1020 { 1021 PetscFunctionBegin; 1022 PetscFunctionReturn(0); 1023 } 1024 1025 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 1026 { 1027 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1028 MPRKTableau tab = mprk->tableau; 1029 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y);CHKERRQ(ierr); 1034 if (mprk->is_slow) { 1035 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 1036 } 1037 ierr = PetscMalloc1(tab->s,&mprk->work_slowbuffer);CHKERRQ(ierr); 1038 if (tab->np == 3) { 1039 if (mprk->is_medium) { 1040 ierr = PetscMalloc1(tab->s,&mprk->work_medium);CHKERRQ(ierr); 1041 } 1042 ierr = PetscMalloc1(tab->s,&mprk->work_mediumbuffer);CHKERRQ(ierr); 1043 } 1044 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 1045 1046 if (ts->use_splitrhsfunction) { 1047 if (mprk->is_slow) { 1048 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1049 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1050 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 1051 } 1052 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1053 ierr = VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1054 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer);CHKERRQ(ierr); 1055 if (tab->np == 3) { 1056 if (mprk->is_medium) { 1057 ierr = VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1058 ierr = VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1059 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium);CHKERRQ(ierr); 1060 } 1061 ierr = VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1062 ierr = VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1063 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer);CHKERRQ(ierr); 1064 } 1065 ierr = VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1066 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1067 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 1068 } else { 1069 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS);CHKERRQ(ierr); 1070 if (mprk->is_slow) { 1071 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 1072 } 1073 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_slowbuffer);CHKERRQ(ierr); 1074 if (tab->np == 3) { 1075 if (mprk->is_medium) { 1076 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_medium);CHKERRQ(ierr); 1077 } 1078 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_mediumbuffer);CHKERRQ(ierr); 1079 } 1080 ierr = PetscMalloc1(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 1081 } 1082 PetscFunctionReturn(0); 1083 } 1084 1085 static PetscErrorCode TSSetUp_MPRK(TS ts) 1086 { 1087 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1088 MPRKTableau tab = mprk->tableau; 1089 DM dm; 1090 PetscErrorCode ierr; 1091 1092 PetscFunctionBegin; 1093 ierr = TSRHSSplitGetIS(ts,"slow",&mprk->is_slow);CHKERRQ(ierr); 1094 ierr = TSRHSSplitGetIS(ts,"fast",&mprk->is_fast);CHKERRQ(ierr); 1095 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); 1096 1097 if (tab->np == 3) { 1098 ierr = TSRHSSplitGetIS(ts,"medium",&mprk->is_medium);CHKERRQ(ierr); 1099 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); 1100 ierr = TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer);CHKERRQ(ierr); 1101 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1102 mprk->is_mediumbuffer = mprk->is_medium; 1103 mprk->is_medium = NULL; 1104 } 1105 } 1106 1107 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1108 ierr = TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer);CHKERRQ(ierr); 1109 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1110 mprk->is_slowbuffer = mprk->is_slow; 1111 mprk->is_slow = NULL; 1112 } 1113 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1114 ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr); 1115 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1116 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1117 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1118 if (ts->use_splitrhsfunction) { 1119 ts->ops->step = TSStep_MPRKSPLIT; 1120 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1121 ierr = TSMPRKSetSplits(ts);CHKERRQ(ierr); 1122 } else { 1123 ts->ops->step = TSStep_MPRK; 1124 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1125 } 1126 PetscFunctionReturn(0); 1127 } 1128 1129 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1130 { 1131 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1132 PetscErrorCode ierr; 1133 1134 PetscFunctionBegin; 1135 ierr = PetscOptionsHead(PetscOptionsObject,"PRK ODE solver options");CHKERRQ(ierr); 1136 { 1137 MPRKTableauLink link; 1138 PetscInt count,choice; 1139 PetscBool flg; 1140 const char **namelist; 1141 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1142 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1143 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1144 ierr = PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1145 if (flg) {ierr = TSMPRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1146 ierr = PetscFree(namelist);CHKERRQ(ierr); 1147 } 1148 ierr = PetscOptionsTail();CHKERRQ(ierr); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1153 { 1154 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1155 PetscBool iascii; 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1160 if (iascii) { 1161 MPRKTableau tab = mprk->tableau; 1162 TSMPRKType mprktype; 1163 char fbuf[512]; 1164 char sbuf[512]; 1165 PetscInt i; 1166 ierr = TSMPRKGetType(ts,&mprktype);CHKERRQ(ierr); 1167 ierr = PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype);CHKERRQ(ierr); 1168 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1169 1170 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf);CHKERRQ(ierr); 1171 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf);CHKERRQ(ierr); 1172 ierr = PetscViewerASCIIPrintf(viewer," Af = \n");CHKERRQ(ierr); 1173 for (i=0; i<tab->s; i++) { 1174 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s]);CHKERRQ(ierr); 1175 ierr = PetscViewerASCIIPrintf(viewer," %s\n",fbuf);CHKERRQ(ierr); 1176 } 1177 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf);CHKERRQ(ierr); 1178 ierr = PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf);CHKERRQ(ierr); 1179 1180 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb);CHKERRQ(ierr); 1181 ierr = PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf);CHKERRQ(ierr); 1182 ierr = PetscViewerASCIIPrintf(viewer," Asb = \n");CHKERRQ(ierr); 1183 for (i=0; i<tab->s; i++) { 1184 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s]);CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," %s\n",sbuf);CHKERRQ(ierr); 1186 } 1187 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb);CHKERRQ(ierr); 1188 ierr = PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf);CHKERRQ(ierr); 1189 1190 if (tab->np == 3) { 1191 char mbuf[512]; 1192 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb);CHKERRQ(ierr); 1193 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf);CHKERRQ(ierr); 1194 ierr = PetscViewerASCIIPrintf(viewer," Amb = \n");CHKERRQ(ierr); 1195 for (i=0; i<tab->s; i++) { 1196 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s]);CHKERRQ(ierr); 1197 ierr = PetscViewerASCIIPrintf(viewer," %s\n",mbuf);CHKERRQ(ierr); 1198 } 1199 ierr = PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb);CHKERRQ(ierr); 1200 ierr = PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf);CHKERRQ(ierr); 1201 } 1202 } 1203 PetscFunctionReturn(0); 1204 } 1205 1206 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1207 { 1208 PetscErrorCode ierr; 1209 TSAdapt adapt; 1210 1211 PetscFunctionBegin; 1212 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1213 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /*@C 1218 TSMPRKSetType - Set the type of MPRK scheme 1219 1220 Not collective 1221 1222 Input Parameter: 1223 + ts - timestepping context 1224 - mprktype - type of MPRK-scheme 1225 1226 Options Database: 1227 . -ts_mprk_type - <pm2,p2,p3> 1228 1229 Level: intermediate 1230 1231 .seealso: TSMPRKGetType(), TSMPRK, TSMPRKType 1232 @*/ 1233 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1234 { 1235 PetscErrorCode ierr; 1236 1237 PetscFunctionBegin; 1238 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1239 PetscValidCharPointer(mprktype,2); 1240 ierr = PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype));CHKERRQ(ierr); 1241 PetscFunctionReturn(0); 1242 } 1243 1244 /*@C 1245 TSMPRKGetType - Get the type of MPRK scheme 1246 1247 Not collective 1248 1249 Input Parameter: 1250 . ts - timestepping context 1251 1252 Output Parameter: 1253 . mprktype - type of MPRK-scheme 1254 1255 Level: intermediate 1256 1257 .seealso: TSMPRKGetType() 1258 @*/ 1259 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1260 { 1261 PetscErrorCode ierr; 1262 1263 PetscFunctionBegin; 1264 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1265 ierr = PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype));CHKERRQ(ierr); 1266 PetscFunctionReturn(0); 1267 } 1268 1269 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1270 { 1271 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1272 1273 PetscFunctionBegin; 1274 *mprktype = mprk->tableau->name; 1275 PetscFunctionReturn(0); 1276 } 1277 1278 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1279 { 1280 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1281 PetscBool match; 1282 MPRKTableauLink link; 1283 PetscErrorCode ierr; 1284 1285 PetscFunctionBegin; 1286 if (mprk->tableau) { 1287 ierr = PetscStrcmp(mprk->tableau->name,mprktype,&match);CHKERRQ(ierr); 1288 if (match) PetscFunctionReturn(0); 1289 } 1290 for (link = MPRKTableauList; link; link=link->next) { 1291 ierr = PetscStrcmp(link->tab.name,mprktype,&match);CHKERRQ(ierr); 1292 if (match) { 1293 if (ts->setupcalled) {ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr);} 1294 mprk->tableau = &link->tab; 1295 if (ts->setupcalled) {ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr);} 1296 PetscFunctionReturn(0); 1297 } 1298 } 1299 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1300 PetscFunctionReturn(0); 1301 } 1302 1303 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1304 { 1305 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1306 1307 PetscFunctionBegin; 1308 *ns = mprk->tableau->s; 1309 if (Y) *Y = mprk->Y; 1310 PetscFunctionReturn(0); 1311 } 1312 1313 static PetscErrorCode TSDestroy_MPRK(TS ts) 1314 { 1315 PetscErrorCode ierr; 1316 1317 PetscFunctionBegin; 1318 ierr = TSReset_MPRK(ts);CHKERRQ(ierr); 1319 if (ts->dm) { 1320 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1321 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 1322 } 1323 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1324 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL);CHKERRQ(ierr); 1325 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL);CHKERRQ(ierr); 1326 PetscFunctionReturn(0); 1327 } 1328 1329 /*MC 1330 TSMPRK - ODE solver using Partitioned Runge-Kutta schemes 1331 1332 The user should provide the right hand side of the equation 1333 using TSSetRHSFunction(). 1334 1335 Notes: 1336 The default is TSMPRKPM2, it can be changed with TSRKSetType() or -ts_mprk_type 1337 1338 Level: beginner 1339 1340 .seealso: TSCreate(), TS, TSSetType(), TSMPRKSetType(), TSMPRKGetType(), TSMPRKType, TSMPRKRegister(), TSMPRKSetMultirateType() 1341 TSMPRKM2, TSMPRKM3, TSMPRKRFSMR3, TSMPRKRFSMR2 1342 1343 M*/ 1344 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1345 { 1346 TS_MPRK *mprk; 1347 PetscErrorCode ierr; 1348 1349 PetscFunctionBegin; 1350 ierr = TSMPRKInitializePackage();CHKERRQ(ierr); 1351 1352 ts->ops->reset = TSReset_MPRK; 1353 ts->ops->destroy = TSDestroy_MPRK; 1354 ts->ops->view = TSView_MPRK; 1355 ts->ops->load = TSLoad_MPRK; 1356 ts->ops->setup = TSSetUp_MPRK; 1357 ts->ops->step = TSStep_MPRK; 1358 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1359 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1360 ts->ops->getstages = TSGetStages_MPRK; 1361 1362 ierr = PetscNewLog(ts,&mprk);CHKERRQ(ierr); 1363 ts->data = (void*)mprk; 1364 1365 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK);CHKERRQ(ierr); 1366 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK);CHKERRQ(ierr); 1367 1368 ierr = TSMPRKSetType(ts,TSMPRKDefault);CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371