1 /* 2 Code for time stepping with the Multirate Partitioned Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) 7 if one does not split the RHS function, but gives 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) 11 for component-wise partitioned system, 12 users should split the RHS function themselves and also provide the indexes for both slow and fast components. 13 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'. 14 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. 15 16 Reference: 17 Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007 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 Multirate 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 136 Options database: 137 . -ts_mprk_type 2a22 - select this scheme 138 139 Level: advanced 140 141 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 142 M*/ 143 /*MC 144 TSMPRK2A23 - Second Order Multirate 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 149 Options database: 150 . -ts_mprk_type 2a23 - select this scheme 151 152 Level: advanced 153 154 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 155 M*/ 156 /*MC 157 TSMPRK2A32 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A. 158 159 This method has four stages for slow and fast parts. The refinement factor of the stepsize is 3. 160 r = 3, np = 2 161 162 Options database: 163 . -ts_mprk_type 2a32 - select this scheme 164 165 Level: advanced 166 167 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 168 M*/ 169 /*MC 170 TSMPRK2A33 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A. 171 172 This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 3. 173 r = 3, np = 3 174 175 Options database: 176 . -ts_mprk_type 2a33- select this scheme 177 178 Level: advanced 179 180 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 181 M*/ 182 /*MC 183 TSMPRK3P2M - Third Order Multirate Partitioned Runge-Kutta scheme. 184 185 This method has eight stages for both slow and fast parts. 186 187 Options database: 188 . -ts_mprk_type pm3 - select this scheme 189 190 Level: advanced 191 192 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 193 M*/ 194 /*MC 195 TSMPRKP2 - Second Order Multirate Partitioned Runge-Kutta scheme. 196 197 This method has five stages for both slow and fast parts. 198 199 Options database: 200 . -ts_mprk_type p2 - select this scheme 201 202 Level: advanced 203 204 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 205 M*/ 206 /*MC 207 TSMPRKP3 - Third Order Multirate Partitioned Runge-Kutta scheme. 208 209 This method has ten stages for both slow and fast parts. 210 211 Options database: 212 . -ts_mprk_type p3 - select this scheme 213 214 Level: advanced 215 216 .seealso: `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()` 217 M*/ 218 219 /*@C 220 TSMPRKRegisterAll - Registers all of the Partirioned Runge-Kutta explicit methods in TSMPRK 221 222 Not Collective, but should be called by all processes which will need the schemes to be registered 223 224 Level: advanced 225 226 .seealso: `TSMPRKRegisterDestroy()` 227 @*/ 228 PetscErrorCode TSMPRKRegisterAll(void) 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 PetscCall(TSMPRKGenerateTableau2(2,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Af[0][0],bf)); 245 PetscCall(TSMPRKRegister(TSMPRK2A22,2,2,2,1,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL)); 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 PetscCall(TSMPRKGenerateTableau3(2,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Amb[0][0],bmb,&Af[0][0],bf)); 257 PetscCall(TSMPRKRegister(TSMPRK2A23,2,2,2,2,&Asb[0][0],bsb,NULL,rsb,&Amb[0][0],bmb,NULL,rmb,&Af[0][0],bf,NULL)); 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 PetscCall(TSMPRKGenerateTableau2(3,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Af[0][0],bf)); 269 PetscCall(TSMPRKRegister(TSMPRK2A32,2,2,3,1,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL)); 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 PetscCall(TSMPRKGenerateTableau3(3,2,&Abase[0][0],bbase,&Asb[0][0],bsb,&Amb[0][0],bmb,&Af[0][0],bf)); 281 PetscCall(TSMPRKRegister(TSMPRK2A33,2,2,3,3,&Asb[0][0],bsb,NULL,rsb,&Amb[0][0],bmb,NULL,rmb,&Af[0][0],bf,NULL)); 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 PetscCall(TSMPRKRegister(TSMPRKPM3,3,8,&As[0][0],bs,NULL,&A[0][0],b,NULL)); 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 PetscCall(TSMPRKRegister(TSMPRKP2,2,5,1,1,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL)); 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 PetscCall(TSMPRKRegister(TSMPRKP3,3,5,2,1,&Asb[0][0],bsb,NULL,rsb,NULL,NULL,NULL,NULL,&Af[0][0],bf,NULL)); 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 .seealso: `TSMPRKRegister()`, `TSMPRKRegisterAll()` 395 @*/ 396 PetscErrorCode TSMPRKRegisterDestroy(void) 397 { 398 MPRKTableauLink link; 399 400 PetscFunctionBegin; 401 while ((link = MPRKTableauList)) { 402 MPRKTableau t = &link->tab; 403 MPRKTableauList = link->next; 404 PetscCall(PetscFree3(t->Asb,t->bsb,t->csb)); 405 PetscCall(PetscFree3(t->Amb,t->bmb,t->cmb)); 406 PetscCall(PetscFree3(t->Af,t->bf,t->cf)); 407 PetscCall(PetscFree(t->rsb)); 408 PetscCall(PetscFree(t->rmb)); 409 PetscCall(PetscFree(t->name)); 410 PetscCall(PetscFree(link)); 411 } 412 TSMPRKRegisterAllCalled = PETSC_FALSE; 413 PetscFunctionReturn(0); 414 } 415 416 /*@C 417 TSMPRKInitializePackage - This function initializes everything in the TSMPRK package. It is called 418 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_MPRK() 419 when using static libraries. 420 421 Level: developer 422 423 .seealso: `PetscInitialize()` 424 @*/ 425 PetscErrorCode TSMPRKInitializePackage(void) 426 { 427 PetscFunctionBegin; 428 if (TSMPRKPackageInitialized) PetscFunctionReturn(0); 429 TSMPRKPackageInitialized = PETSC_TRUE; 430 PetscCall(TSMPRKRegisterAll()); 431 PetscCall(PetscRegisterFinalize(TSMPRKFinalizePackage)); 432 PetscFunctionReturn(0); 433 } 434 435 /*@C 436 TSMPRKFinalizePackage - This function destroys everything in the TSMPRK package. It is 437 called from PetscFinalize(). 438 439 Level: developer 440 441 .seealso: `PetscFinalize()` 442 @*/ 443 PetscErrorCode TSMPRKFinalizePackage(void) 444 { 445 PetscFunctionBegin; 446 TSMPRKPackageInitialized = PETSC_FALSE; 447 PetscCall(TSMPRKRegisterDestroy()); 448 PetscFunctionReturn(0); 449 } 450 451 /*@C 452 TSMPRKRegister - register a MPRK scheme by providing the entries in the Butcher tableau 453 454 Not Collective, but the same schemes should be registered on all processes on which they will be used 455 456 Input Parameters: 457 + name - identifier for method 458 . order - approximation order of method 459 . s - number of stages in the base methods 460 . ratio1 - stepsize ratio at 1st level (e.g. slow/medium) 461 . ratio2 - stepsize ratio at 2nd level (e.g. medium/fast) 462 . Af - stage coefficients for fast components(dimension s*s, row-major) 463 . bf - step completion table for fast components(dimension s) 464 . cf - abscissa for fast components(dimension s) 465 . As - stage coefficients for slow components(dimension s*s, row-major) 466 . bs - step completion table for slow components(dimension s) 467 - cs - abscissa for slow components(dimension s) 468 469 Notes: 470 Several MPRK methods are provided, this function is only needed to create new methods. 471 472 Level: advanced 473 474 .seealso: `TSMPRK` 475 @*/ 476 PetscErrorCode TSMPRKRegister(TSMPRKType name,PetscInt order, 477 PetscInt sbase,PetscInt ratio1,PetscInt ratio2, 478 const PetscReal Asb[],const PetscReal bsb[],const PetscReal csb[],const PetscInt rsb[], 479 const PetscReal Amb[],const PetscReal bmb[],const PetscReal cmb[],const PetscInt rmb[], 480 const PetscReal Af[],const PetscReal bf[],const PetscReal cf[]) 481 { 482 MPRKTableauLink link; 483 MPRKTableau t; 484 PetscInt s,i,j; 485 486 PetscFunctionBegin; 487 PetscValidCharPointer(name,1); 488 PetscValidRealPointer(Asb,6); 489 if (bsb) PetscValidRealPointer(bsb,7); 490 if (csb) PetscValidRealPointer(csb,8); 491 if (rsb) PetscValidIntPointer(rsb,9); 492 if (Amb) PetscValidRealPointer(Amb,10); 493 if (bmb) PetscValidRealPointer(bmb,11); 494 if (cmb) PetscValidRealPointer(cmb,12); 495 if (rmb) PetscValidIntPointer(rmb,13); 496 PetscValidRealPointer(Af,14); 497 if (bf) PetscValidRealPointer(bf,15); 498 if (cf) PetscValidRealPointer(cf,16); 499 500 PetscCall(PetscNew(&link)); 501 t = &link->tab; 502 503 PetscCall(PetscStrallocpy(name,&t->name)); 504 s = sbase*ratio1*ratio2; /* this is the dimension of the matrices below */ 505 t->order = order; 506 t->sbase = sbase; 507 t->s = s; 508 t->np = 2; 509 510 PetscCall(PetscMalloc3(s*s,&t->Af,s,&t->bf,s,&t->cf)); 511 PetscCall(PetscArraycpy(t->Af,Af,s*s)); 512 if (bf) { 513 PetscCall(PetscArraycpy(t->bf,bf,s)); 514 } else 515 for (i=0; i<s; i++) t->bf[i] = Af[(s-1)*s+i]; 516 if (cf) { 517 PetscCall(PetscArraycpy(t->cf,cf,s)); 518 } else { 519 for (i=0; i<s; i++) 520 for (j=0,t->cf[i]=0; j<s; j++) 521 t->cf[i] += Af[i*s+j]; 522 } 523 524 if (Amb) { 525 t->np = 3; 526 PetscCall(PetscMalloc3(s*s,&t->Amb,s,&t->bmb,s,&t->cmb)); 527 PetscCall(PetscArraycpy(t->Amb,Amb,s*s)); 528 if (bmb) { 529 PetscCall(PetscArraycpy(t->bmb,bmb,s)); 530 } else { 531 for (i=0; i<s; i++) t->bmb[i] = Amb[(s-1)*s+i]; 532 } 533 if (cmb) { 534 PetscCall(PetscArraycpy(t->cmb,cmb,s)); 535 } else { 536 for (i=0; i<s; i++) 537 for (j=0,t->cmb[i]=0; j<s; j++) 538 t->cmb[i] += Amb[i*s+j]; 539 } 540 if (rmb) { 541 PetscCall(PetscMalloc1(s,&t->rmb)); 542 PetscCall(PetscArraycpy(t->rmb,rmb,s)); 543 } else { 544 PetscCall(PetscCalloc1(s,&t->rmb)); 545 } 546 } 547 548 PetscCall(PetscMalloc3(s*s,&t->Asb,s,&t->bsb,s,&t->csb)); 549 PetscCall(PetscArraycpy(t->Asb,Asb,s*s)); 550 if (bsb) { 551 PetscCall(PetscArraycpy(t->bsb,bsb,s)); 552 } else 553 for (i=0; i<s; i++) t->bsb[i] = Asb[(s-1)*s+i]; 554 if (csb) { 555 PetscCall(PetscArraycpy(t->csb,csb,s)); 556 } else { 557 for (i=0; i<s; i++) 558 for (j=0,t->csb[i]=0; j<s; j++) 559 t->csb[i] += Asb[i*s+j]; 560 } 561 if (rsb) { 562 PetscCall(PetscMalloc1(s,&t->rsb)); 563 PetscCall(PetscArraycpy(t->rsb,rsb,s)); 564 } else { 565 PetscCall(PetscCalloc1(s,&t->rsb)); 566 } 567 link->next = MPRKTableauList; 568 MPRKTableauList = link; 569 PetscFunctionReturn(0); 570 } 571 572 static PetscErrorCode TSMPRKSetSplits(TS ts) 573 { 574 TS_MPRK *mprk = (TS_MPRK*)ts->data; 575 MPRKTableau tab = mprk->tableau; 576 DM dm,subdm,newdm; 577 578 PetscFunctionBegin; 579 PetscCall(TSRHSSplitGetSubTS(ts,"slow",&mprk->subts_slow)); 580 PetscCall(TSRHSSplitGetSubTS(ts,"fast",&mprk->subts_fast)); 581 PetscCheck(mprk->subts_slow && mprk->subts_fast,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"); 582 583 /* Only copy the DM */ 584 PetscCall(TSGetDM(ts,&dm)); 585 586 PetscCall(TSRHSSplitGetSubTS(ts,"slowbuffer",&mprk->subts_slowbuffer)); 587 if (!mprk->subts_slowbuffer) { 588 mprk->subts_slowbuffer = mprk->subts_slow; 589 mprk->subts_slow = NULL; 590 } 591 if (mprk->subts_slow) { 592 PetscCall(DMClone(dm,&newdm)); 593 PetscCall(TSGetDM(mprk->subts_slow,&subdm)); 594 PetscCall(DMCopyDMTS(subdm,newdm)); 595 PetscCall(DMCopyDMSNES(subdm,newdm)); 596 PetscCall(TSSetDM(mprk->subts_slow,newdm)); 597 PetscCall(DMDestroy(&newdm)); 598 } 599 PetscCall(DMClone(dm,&newdm)); 600 PetscCall(TSGetDM(mprk->subts_slowbuffer,&subdm)); 601 PetscCall(DMCopyDMTS(subdm,newdm)); 602 PetscCall(DMCopyDMSNES(subdm,newdm)); 603 PetscCall(TSSetDM(mprk->subts_slowbuffer,newdm)); 604 PetscCall(DMDestroy(&newdm)); 605 606 PetscCall(DMClone(dm,&newdm)); 607 PetscCall(TSGetDM(mprk->subts_fast,&subdm)); 608 PetscCall(DMCopyDMTS(subdm,newdm)); 609 PetscCall(DMCopyDMSNES(subdm,newdm)); 610 PetscCall(TSSetDM(mprk->subts_fast,newdm)); 611 PetscCall(DMDestroy(&newdm)); 612 613 if (tab->np == 3) { 614 PetscCall(TSRHSSplitGetSubTS(ts,"medium",&mprk->subts_medium)); 615 PetscCall(TSRHSSplitGetSubTS(ts,"mediumbuffer",&mprk->subts_mediumbuffer)); 616 if (mprk->subts_medium && !mprk->subts_mediumbuffer) { 617 mprk->subts_mediumbuffer = mprk->subts_medium; 618 mprk->subts_medium = NULL; 619 } 620 if (mprk->subts_medium) { 621 PetscCall(DMClone(dm,&newdm)); 622 PetscCall(TSGetDM(mprk->subts_medium,&subdm)); 623 PetscCall(DMCopyDMTS(subdm,newdm)); 624 PetscCall(DMCopyDMSNES(subdm,newdm)); 625 PetscCall(TSSetDM(mprk->subts_medium,newdm)); 626 PetscCall(DMDestroy(&newdm)); 627 } 628 PetscCall(DMClone(dm,&newdm)); 629 PetscCall(TSGetDM(mprk->subts_mediumbuffer,&subdm)); 630 PetscCall(DMCopyDMTS(subdm,newdm)); 631 PetscCall(DMCopyDMSNES(subdm,newdm)); 632 PetscCall(TSSetDM(mprk->subts_mediumbuffer,newdm)); 633 PetscCall(DMDestroy(&newdm)); 634 } 635 PetscFunctionReturn(0); 636 } 637 638 /* 639 This if for nonsplit RHS MPRK 640 The step completion formula is 641 642 x1 = x0 + h b^T YdotRHS 643 644 */ 645 static PetscErrorCode TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool *done) 646 { 647 TS_MPRK *mprk = (TS_MPRK*)ts->data; 648 MPRKTableau tab = mprk->tableau; 649 PetscScalar *wf = mprk->work_fast; 650 PetscReal h = ts->time_step; 651 PetscInt s = tab->s,j; 652 653 PetscFunctionBegin; 654 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 655 PetscCall(VecCopy(ts->vec_sol,X)); 656 PetscCall(VecMAXPY(X,s,wf,mprk->YdotRHS)); 657 PetscFunctionReturn(0); 658 } 659 660 static PetscErrorCode TSStep_MPRK(TS ts) 661 { 662 TS_MPRK *mprk = (TS_MPRK*)ts->data; 663 Vec *Y = mprk->Y,*YdotRHS = mprk->YdotRHS,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 664 Vec Yslow,Yslowbuffer,Yfast; 665 MPRKTableau tab = mprk->tableau; 666 const PetscInt s = tab->s; 667 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 668 PetscScalar *wf = mprk->work_fast,*wsb = mprk->work_slowbuffer; 669 PetscInt i,j; 670 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 671 672 PetscFunctionBegin; 673 for (i=0; i<s; i++) { 674 mprk->stage_time = t + h*cf[i]; 675 PetscCall(TSPreStage(ts,mprk->stage_time)); 676 PetscCall(VecCopy(ts->vec_sol,Y[i])); 677 678 /* slow buffer region */ 679 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 680 for (j=0; j<i; j++) { 681 PetscCall(VecGetSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j])); 682 } 683 PetscCall(VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 684 PetscCall(VecMAXPY(Yslowbuffer,i,wsb,mprk->YdotRHS_slowbuffer)); 685 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 686 for (j=0; j<i; j++) { 687 PetscCall(VecRestoreSubVector(YdotRHS[j],mprk->is_slowbuffer,&YdotRHS_slowbuffer[j])); 688 } 689 /* slow region */ 690 if (mprk->is_slow) { 691 for (j=0; j<i; j++) { 692 PetscCall(VecGetSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j])); 693 } 694 PetscCall(VecGetSubVector(Y[i],mprk->is_slow,&Yslow)); 695 PetscCall(VecMAXPY(Yslow,i,wsb,mprk->YdotRHS_slow)); 696 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow)); 697 for (j=0; j<i; j++) { 698 PetscCall(VecRestoreSubVector(YdotRHS[j],mprk->is_slow,&YdotRHS_slow[j])); 699 } 700 } 701 702 /* fast region */ 703 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 704 for (j=0; j<i; j++) { 705 PetscCall(VecGetSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j])); 706 } 707 PetscCall(VecGetSubVector(Y[i],mprk->is_fast,&Yfast)); 708 PetscCall(VecMAXPY(Yfast,i,wf,mprk->YdotRHS_fast)); 709 PetscCall(VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast)); 710 for (j=0; j<i; j++) { 711 PetscCall(VecRestoreSubVector(YdotRHS[j],mprk->is_fast,&YdotRHS_fast[j])); 712 } 713 if (tab->np == 3) { 714 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 715 Vec Ymedium,Ymediumbuffer; 716 PetscScalar *wmb = mprk->work_mediumbuffer; 717 718 for (j=0; j<i; j++) wmb[j] = h*tab->Amb[i*s+j]; 719 /* medium region */ 720 if (mprk->is_medium) { 721 for (j=0; j<i; j++) { 722 PetscCall(VecGetSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j])); 723 } 724 PetscCall(VecGetSubVector(Y[i],mprk->is_medium,&Ymedium)); 725 PetscCall(VecMAXPY(Ymedium,i,wmb,mprk->YdotRHS_medium)); 726 PetscCall(VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium)); 727 for (j=0; j<i; j++) { 728 PetscCall(VecRestoreSubVector(YdotRHS[j],mprk->is_medium,&YdotRHS_medium[j])); 729 } 730 } 731 /* medium buffer region */ 732 for (j=0; j<i; j++) { 733 PetscCall(VecGetSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j])); 734 } 735 PetscCall(VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer)); 736 PetscCall(VecMAXPY(Ymediumbuffer,i,wmb,mprk->YdotRHS_mediumbuffer)); 737 PetscCall(VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer)); 738 for (j=0; j<i; j++) { 739 PetscCall(VecRestoreSubVector(YdotRHS[j],mprk->is_mediumbuffer,&YdotRHS_mediumbuffer[j])); 740 } 741 } 742 PetscCall(TSPostStage(ts,mprk->stage_time,i,Y)); 743 /* compute the stage RHS by fast and slow tableau respectively */ 744 PetscCall(TSComputeRHSFunction(ts,t+h*csb[i],Y[i],YdotRHS[i])); 745 } 746 PetscCall(TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL)); 747 ts->ptime += ts->time_step; 748 ts->time_step = next_time_step; 749 PetscFunctionReturn(0); 750 } 751 752 /* 753 This if for the case when split RHS is used 754 The step completion formula is 755 x1 = x0 + h b^T YdotRHS 756 */ 757 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 758 { 759 TS_MPRK *mprk = (TS_MPRK*)ts->data; 760 MPRKTableau tab = mprk->tableau; 761 Vec Xslow,Xfast,Xslowbuffer; /* subvectors for slow and fast componets in X respectively */ 762 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 763 PetscReal h = ts->time_step; 764 PetscInt s = tab->s,j,computedstages; 765 766 PetscFunctionBegin; 767 PetscCall(VecCopy(ts->vec_sol,X)); 768 769 /* slow region */ 770 if (mprk->is_slow) { 771 computedstages = 0; 772 for (j=0; j<s; j++) { 773 if (tab->rsb[j]) ws[tab->rsb[j]-1] += h*tab->bsb[j]; 774 else ws[computedstages++] = h*tab->bsb[j]; 775 } 776 PetscCall(VecGetSubVector(X,mprk->is_slow,&Xslow)); 777 PetscCall(VecMAXPY(Xslow,computedstages,ws,mprk->YdotRHS_slow)); 778 PetscCall(VecRestoreSubVector(X,mprk->is_slow,&Xslow)); 779 } 780 781 if (tab->np == 3 && mprk->is_medium) { 782 computedstages = 0; 783 for (j=0; j<s; j++) { 784 if (tab->rmb[j]) wsb[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bsb[j]; 785 else wsb[computedstages++] = h*tab->bsb[j]; 786 } 787 PetscCall(VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer)); 788 PetscCall(VecMAXPY(Xslowbuffer,computedstages,wsb,mprk->YdotRHS_slowbuffer)); 789 PetscCall(VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer)); 790 } else { 791 /* slow buffer region */ 792 for (j=0; j<s; j++) wsb[j] = h*tab->bsb[j]; 793 PetscCall(VecGetSubVector(X,mprk->is_slowbuffer,&Xslowbuffer)); 794 PetscCall(VecMAXPY(Xslowbuffer,s,wsb,mprk->YdotRHS_slowbuffer)); 795 PetscCall(VecRestoreSubVector(X,mprk->is_slowbuffer,&Xslowbuffer)); 796 } 797 if (tab->np == 3) { 798 Vec Xmedium,Xmediumbuffer; 799 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 800 /* medium region and slow buffer region */ 801 if (mprk->is_medium) { 802 computedstages = 0; 803 for (j=0; j<s; j++) { 804 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += h*tab->bmb[j]; 805 else wm[computedstages++] = h*tab->bmb[j]; 806 } 807 PetscCall(VecGetSubVector(X,mprk->is_medium,&Xmedium)); 808 PetscCall(VecMAXPY(Xmedium,computedstages,wm,mprk->YdotRHS_medium)); 809 PetscCall(VecRestoreSubVector(X,mprk->is_medium,&Xmedium)); 810 } 811 /* medium buffer region */ 812 for (j=0; j<s; j++) wmb[j] = h*tab->bmb[j]; 813 PetscCall(VecGetSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer)); 814 PetscCall(VecMAXPY(Xmediumbuffer,s,wmb,mprk->YdotRHS_mediumbuffer)); 815 PetscCall(VecRestoreSubVector(X,mprk->is_mediumbuffer,&Xmediumbuffer)); 816 } 817 /* fast region */ 818 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 819 PetscCall(VecGetSubVector(X,mprk->is_fast,&Xfast)); 820 PetscCall(VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast)); 821 PetscCall(VecRestoreSubVector(X,mprk->is_fast,&Xfast)); 822 PetscFunctionReturn(0); 823 } 824 825 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 826 { 827 TS_MPRK *mprk = (TS_MPRK*)ts->data; 828 MPRKTableau tab = mprk->tableau; 829 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow,*YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 830 Vec Yslow,Yslowbuffer,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 831 PetscInt s = tab->s; 832 const PetscReal *Af = tab->Af,*cf = tab->cf,*Asb = tab->Asb,*csb = tab->csb; 833 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow,*wsb = mprk->work_slowbuffer; 834 PetscInt i,j,computedstages; 835 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 836 837 PetscFunctionBegin; 838 for (i=0; i<s; i++) { 839 mprk->stage_time = t + h*cf[i]; 840 PetscCall(TSPreStage(ts,mprk->stage_time)); 841 /* calculate the stage value for fast and slow components respectively */ 842 PetscCall(VecCopy(ts->vec_sol,Y[i])); 843 for (j=0; j<i; j++) wsb[j] = h*Asb[i*s+j]; 844 845 /* slow buffer region */ 846 if (tab->np == 3 && mprk->is_medium) { 847 if (tab->rmb[i]) { 848 PetscCall(VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 849 PetscCall(VecISCopy(Y[tab->rmb[i]-1],mprk->is_slowbuffer,SCATTER_REVERSE,Yslowbuffer)); 850 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 851 } else { 852 PetscScalar *wm = mprk->work_medium; 853 computedstages = 0; 854 for (j=0; j<i; j++) { 855 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wsb[j]; 856 else wm[computedstages++] = wsb[j]; 857 } 858 PetscCall(VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 859 PetscCall(VecMAXPY(Yslowbuffer,computedstages,wm,YdotRHS_slowbuffer)); 860 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 861 } 862 } else { 863 PetscCall(VecGetSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 864 PetscCall(VecMAXPY(Yslowbuffer,i,wsb,YdotRHS_slowbuffer)); 865 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slowbuffer,&Yslowbuffer)); 866 } 867 868 /* slow region */ 869 if (mprk->is_slow) { 870 if (tab->rsb[i]) { /* repeat previous stage */ 871 PetscCall(VecGetSubVector(Y[i],mprk->is_slow,&Yslow)); 872 PetscCall(VecISCopy(Y[tab->rsb[i]-1],mprk->is_slow,SCATTER_REVERSE,Yslow)); 873 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow)); 874 } else { 875 computedstages = 0; 876 for (j=0; j<i; j++) { 877 if (tab->rsb[j]) ws[tab->rsb[j]-1] += wsb[j]; 878 else ws[computedstages++] = wsb[j]; 879 } 880 PetscCall(VecGetSubVector(Y[i],mprk->is_slow,&Yslow)); 881 PetscCall(VecMAXPY(Yslow,computedstages,ws,YdotRHS_slow)); 882 PetscCall(VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow)); 883 /* only depends on the slow buffer region */ 884 PetscCall(TSComputeRHSFunction(mprk->subts_slow,t+h*csb[i],Y[i],YdotRHS_slow[computedstages])); 885 } 886 } 887 888 /* fast region */ 889 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 890 PetscCall(VecGetSubVector(Y[i],mprk->is_fast,&Yfast)); 891 PetscCall(VecMAXPY(Yfast,i,wf,YdotRHS_fast)); 892 PetscCall(VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast)); 893 894 if (tab->np == 3) { 895 Vec *YdotRHS_medium = mprk->YdotRHS_medium,*YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 896 Vec Ymedium,Ymediumbuffer; 897 const PetscReal *Amb = tab->Amb,*cmb = tab->cmb; 898 PetscScalar *wm = mprk->work_medium,*wmb = mprk->work_mediumbuffer; 899 900 for (j=0; j<i; j++) wmb[j] = h*Amb[i*s+j]; 901 /* medium buffer region */ 902 PetscCall(VecGetSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer)); 903 PetscCall(VecMAXPY(Ymediumbuffer,i,wmb,YdotRHS_mediumbuffer)); 904 PetscCall(VecRestoreSubVector(Y[i],mprk->is_mediumbuffer,&Ymediumbuffer)); 905 /* medium region */ 906 if (mprk->is_medium) { 907 if (tab->rmb[i]) { /* repeat previous stage */ 908 PetscCall(VecGetSubVector(Y[i],mprk->is_medium,&Ymedium)); 909 PetscCall(VecISCopy(Y[tab->rmb[i]-1],mprk->is_medium,SCATTER_REVERSE,Ymedium)); 910 PetscCall(VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium)); 911 } else { 912 computedstages = 0; 913 for (j=0; j<i; j++) { 914 if (tab->rmb[j]) wm[computedstages-tab->sbase+(tab->rmb[j]-1)%tab->sbase] += wmb[j]; 915 else wm[computedstages++] = wmb[j]; 916 917 } 918 PetscCall(VecGetSubVector(Y[i],mprk->is_medium,&Ymedium)); 919 PetscCall(VecMAXPY(Ymedium,computedstages,wm,YdotRHS_medium)); 920 PetscCall(VecRestoreSubVector(Y[i],mprk->is_medium,&Ymedium)); 921 /* only depends on the medium buffer region */ 922 PetscCall(TSComputeRHSFunction(mprk->subts_medium,t+h*cmb[i],Y[i],YdotRHS_medium[computedstages])); 923 /* must be computed after all regions are updated in Y */ 924 PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[computedstages])); 925 } 926 } 927 /* must be computed after fast region and slow region are updated in Y */ 928 PetscCall(TSComputeRHSFunction(mprk->subts_mediumbuffer,t+h*cmb[i],Y[i],YdotRHS_mediumbuffer[i])); 929 } 930 if (!(tab->np == 3 && mprk->is_medium)) { 931 PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer,t+h*csb[i],Y[i],YdotRHS_slowbuffer[i])); 932 } 933 PetscCall(TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i])); 934 } 935 936 PetscCall(TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL)); 937 ts->ptime += ts->time_step; 938 ts->time_step = next_time_step; 939 PetscFunctionReturn(0); 940 } 941 942 static PetscErrorCode TSMPRKTableauReset(TS ts) 943 { 944 TS_MPRK *mprk = (TS_MPRK*)ts->data; 945 MPRKTableau tab = mprk->tableau; 946 947 PetscFunctionBegin; 948 if (!tab) PetscFunctionReturn(0); 949 PetscCall(PetscFree(mprk->work_fast)); 950 PetscCall(PetscFree(mprk->work_slow)); 951 PetscCall(PetscFree(mprk->work_slowbuffer)); 952 PetscCall(PetscFree(mprk->work_medium)); 953 PetscCall(PetscFree(mprk->work_mediumbuffer)); 954 PetscCall(VecDestroyVecs(tab->s,&mprk->Y)); 955 if (ts->use_splitrhsfunction) { 956 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS_fast)); 957 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS_slow)); 958 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS_slowbuffer)); 959 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS_medium)); 960 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS_mediumbuffer)); 961 } else { 962 PetscCall(VecDestroyVecs(tab->s,&mprk->YdotRHS)); 963 if (mprk->is_slow) { 964 PetscCall(PetscFree(mprk->YdotRHS_slow)); 965 } 966 PetscCall(PetscFree(mprk->YdotRHS_slowbuffer)); 967 if (tab->np == 3) { 968 if (mprk->is_medium) { 969 PetscCall(PetscFree(mprk->YdotRHS_medium)); 970 } 971 PetscCall(PetscFree(mprk->YdotRHS_mediumbuffer)); 972 } 973 PetscCall(PetscFree(mprk->YdotRHS_fast)); 974 } 975 PetscFunctionReturn(0); 976 } 977 978 static PetscErrorCode TSReset_MPRK(TS ts) 979 { 980 PetscFunctionBegin; 981 PetscCall(TSMPRKTableauReset(ts)); 982 PetscFunctionReturn(0); 983 } 984 985 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 986 { 987 PetscFunctionBegin; 988 PetscFunctionReturn(0); 989 } 990 991 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 992 { 993 PetscFunctionBegin; 994 PetscFunctionReturn(0); 995 } 996 997 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 998 { 999 PetscFunctionBegin; 1000 PetscFunctionReturn(0); 1001 } 1002 1003 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1004 { 1005 PetscFunctionBegin; 1006 PetscFunctionReturn(0); 1007 } 1008 1009 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 1010 { 1011 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1012 MPRKTableau tab = mprk->tableau; 1013 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 1014 1015 PetscFunctionBegin; 1016 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y)); 1017 if (mprk->is_slow) { 1018 PetscCall(PetscMalloc1(tab->s,&mprk->work_slow)); 1019 } 1020 PetscCall(PetscMalloc1(tab->s,&mprk->work_slowbuffer)); 1021 if (tab->np == 3) { 1022 if (mprk->is_medium) { 1023 PetscCall(PetscMalloc1(tab->s,&mprk->work_medium)); 1024 } 1025 PetscCall(PetscMalloc1(tab->s,&mprk->work_mediumbuffer)); 1026 } 1027 PetscCall(PetscMalloc1(tab->s,&mprk->work_fast)); 1028 1029 if (ts->use_splitrhsfunction) { 1030 if (mprk->is_slow) { 1031 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow)); 1032 PetscCall(VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow)); 1033 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow)); 1034 } 1035 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer)); 1036 PetscCall(VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer)); 1037 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer)); 1038 if (tab->np == 3) { 1039 if (mprk->is_medium) { 1040 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium)); 1041 PetscCall(VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium)); 1042 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium)); 1043 } 1044 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer)); 1045 PetscCall(VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer)); 1046 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer)); 1047 } 1048 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast)); 1049 PetscCall(VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast)); 1050 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast)); 1051 } else { 1052 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS)); 1053 if (mprk->is_slow) { 1054 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_slow)); 1055 } 1056 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_slowbuffer)); 1057 if (tab->np == 3) { 1058 if (mprk->is_medium) { 1059 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_medium)); 1060 } 1061 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_mediumbuffer)); 1062 } 1063 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_fast)); 1064 } 1065 PetscFunctionReturn(0); 1066 } 1067 1068 static PetscErrorCode TSSetUp_MPRK(TS ts) 1069 { 1070 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1071 MPRKTableau tab = mprk->tableau; 1072 DM dm; 1073 1074 PetscFunctionBegin; 1075 PetscCall(TSRHSSplitGetIS(ts,"slow",&mprk->is_slow)); 1076 PetscCall(TSRHSSplitGetIS(ts,"fast",&mprk->is_fast)); 1077 PetscCheck(mprk->is_slow && mprk->is_fast,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); 1078 1079 if (tab->np == 3) { 1080 PetscCall(TSRHSSplitGetIS(ts,"medium",&mprk->is_medium)); 1081 PetscCheck(mprk->is_medium,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); 1082 PetscCall(TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer)); 1083 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1084 mprk->is_mediumbuffer = mprk->is_medium; 1085 mprk->is_medium = NULL; 1086 } 1087 } 1088 1089 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1090 PetscCall(TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer)); 1091 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1092 mprk->is_slowbuffer = mprk->is_slow; 1093 mprk->is_slow = NULL; 1094 } 1095 PetscCall(TSCheckImplicitTerm(ts)); 1096 PetscCall(TSMPRKTableauSetUp(ts)); 1097 PetscCall(TSGetDM(ts,&dm)); 1098 PetscCall(DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts)); 1099 PetscCall(DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts)); 1100 if (ts->use_splitrhsfunction) { 1101 ts->ops->step = TSStep_MPRKSPLIT; 1102 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1103 PetscCall(TSMPRKSetSplits(ts)); 1104 } else { 1105 ts->ops->step = TSStep_MPRK; 1106 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1107 } 1108 PetscFunctionReturn(0); 1109 } 1110 1111 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1112 { 1113 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1114 1115 PetscFunctionBegin; 1116 PetscOptionsHeadBegin(PetscOptionsObject,"PRK ODE solver options"); 1117 { 1118 MPRKTableauLink link; 1119 PetscInt count,choice; 1120 PetscBool flg; 1121 const char **namelist; 1122 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1123 PetscCall(PetscMalloc1(count,(char***)&namelist)); 1124 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1125 PetscCall(PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg)); 1126 if (flg) PetscCall(TSMPRKSetType(ts,namelist[choice])); 1127 PetscCall(PetscFree(namelist)); 1128 } 1129 PetscOptionsHeadEnd(); 1130 PetscFunctionReturn(0); 1131 } 1132 1133 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1134 { 1135 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1136 PetscBool iascii; 1137 1138 PetscFunctionBegin; 1139 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1140 if (iascii) { 1141 MPRKTableau tab = mprk->tableau; 1142 TSMPRKType mprktype; 1143 char fbuf[512]; 1144 char sbuf[512]; 1145 PetscInt i; 1146 PetscCall(TSMPRKGetType(ts,&mprktype)); 1147 PetscCall(PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype)); 1148 PetscCall(PetscViewerASCIIPrintf(viewer," Order: %" PetscInt_FMT "\n",tab->order)); 1149 1150 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf)); 1151 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf)); 1152 PetscCall(PetscViewerASCIIPrintf(viewer," Af = \n")); 1153 for (i=0; i<tab->s; i++) { 1154 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s])); 1155 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",fbuf)); 1156 } 1157 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf)); 1158 PetscCall(PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf)); 1159 1160 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb)); 1161 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf)); 1162 PetscCall(PetscViewerASCIIPrintf(viewer," Asb = \n")); 1163 for (i=0; i<tab->s; i++) { 1164 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s])); 1165 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",sbuf)); 1166 } 1167 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb)); 1168 PetscCall(PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf)); 1169 1170 if (tab->np == 3) { 1171 char mbuf[512]; 1172 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb)); 1173 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf)); 1174 PetscCall(PetscViewerASCIIPrintf(viewer," Amb = \n")); 1175 for (i=0; i<tab->s; i++) { 1176 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s])); 1177 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",mbuf)); 1178 } 1179 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb)); 1180 PetscCall(PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf)); 1181 } 1182 } 1183 PetscFunctionReturn(0); 1184 } 1185 1186 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1187 { 1188 TSAdapt adapt; 1189 1190 PetscFunctionBegin; 1191 PetscCall(TSGetAdapt(ts,&adapt)); 1192 PetscCall(TSAdaptLoad(adapt,viewer)); 1193 PetscFunctionReturn(0); 1194 } 1195 1196 /*@C 1197 TSMPRKSetType - Set the type of MPRK scheme 1198 1199 Not collective 1200 1201 Input Parameters: 1202 + ts - timestepping context 1203 - mprktype - type of MPRK-scheme 1204 1205 Options Database: 1206 . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme 1207 1208 Level: intermediate 1209 1210 .seealso: `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType` 1211 @*/ 1212 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1213 { 1214 PetscFunctionBegin; 1215 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1216 PetscValidCharPointer(mprktype,2); 1217 PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype)); 1218 PetscFunctionReturn(0); 1219 } 1220 1221 /*@C 1222 TSMPRKGetType - Get the type of MPRK scheme 1223 1224 Not collective 1225 1226 Input Parameter: 1227 . ts - timestepping context 1228 1229 Output Parameter: 1230 . mprktype - type of MPRK-scheme 1231 1232 Level: intermediate 1233 1234 .seealso: `TSMPRKGetType()` 1235 @*/ 1236 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1237 { 1238 PetscFunctionBegin; 1239 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1240 PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype)); 1241 PetscFunctionReturn(0); 1242 } 1243 1244 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1245 { 1246 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1247 1248 PetscFunctionBegin; 1249 *mprktype = mprk->tableau->name; 1250 PetscFunctionReturn(0); 1251 } 1252 1253 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1254 { 1255 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1256 PetscBool match; 1257 MPRKTableauLink link; 1258 1259 PetscFunctionBegin; 1260 if (mprk->tableau) { 1261 PetscCall(PetscStrcmp(mprk->tableau->name,mprktype,&match)); 1262 if (match) PetscFunctionReturn(0); 1263 } 1264 for (link = MPRKTableauList; link; link=link->next) { 1265 PetscCall(PetscStrcmp(link->tab.name,mprktype,&match)); 1266 if (match) { 1267 if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts)); 1268 mprk->tableau = &link->tab; 1269 if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts)); 1270 PetscFunctionReturn(0); 1271 } 1272 } 1273 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1274 } 1275 1276 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1277 { 1278 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1279 1280 PetscFunctionBegin; 1281 *ns = mprk->tableau->s; 1282 if (Y) *Y = mprk->Y; 1283 PetscFunctionReturn(0); 1284 } 1285 1286 static PetscErrorCode TSDestroy_MPRK(TS ts) 1287 { 1288 PetscFunctionBegin; 1289 PetscCall(TSReset_MPRK(ts)); 1290 if (ts->dm) { 1291 PetscCall(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts)); 1292 PetscCall(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts)); 1293 } 1294 PetscCall(PetscFree(ts->data)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL)); 1297 PetscFunctionReturn(0); 1298 } 1299 1300 /*MC 1301 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes 1302 1303 The user should provide the right hand side of the equation 1304 using TSSetRHSFunction(). 1305 1306 Notes: 1307 The default is TSMPRKPM2, it can be changed with TSMPRKSetType() or -ts_mprk_type 1308 1309 Level: beginner 1310 1311 .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()` 1312 `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2` 1313 1314 M*/ 1315 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1316 { 1317 TS_MPRK *mprk; 1318 1319 PetscFunctionBegin; 1320 PetscCall(TSMPRKInitializePackage()); 1321 1322 ts->ops->reset = TSReset_MPRK; 1323 ts->ops->destroy = TSDestroy_MPRK; 1324 ts->ops->view = TSView_MPRK; 1325 ts->ops->load = TSLoad_MPRK; 1326 ts->ops->setup = TSSetUp_MPRK; 1327 ts->ops->step = TSStep_MPRK; 1328 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1329 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1330 ts->ops->getstages = TSGetStages_MPRK; 1331 1332 PetscCall(PetscNewLog(ts,&mprk)); 1333 ts->data = (void*)mprk; 1334 1335 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK)); 1336 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK)); 1337 1338 PetscCall(TSMPRKSetType(ts,TSMPRKDefault)); 1339 PetscFunctionReturn(0); 1340 } 1341