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