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 components 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) PetscCall(PetscFree(mprk->YdotRHS_slow)); 964 PetscCall(PetscFree(mprk->YdotRHS_slowbuffer)); 965 if (tab->np == 3) { 966 if (mprk->is_medium) PetscCall(PetscFree(mprk->YdotRHS_medium)); 967 PetscCall(PetscFree(mprk->YdotRHS_mediumbuffer)); 968 } 969 PetscCall(PetscFree(mprk->YdotRHS_fast)); 970 } 971 PetscFunctionReturn(0); 972 } 973 974 static PetscErrorCode TSReset_MPRK(TS ts) 975 { 976 PetscFunctionBegin; 977 PetscCall(TSMPRKTableauReset(ts)); 978 PetscFunctionReturn(0); 979 } 980 981 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 982 { 983 PetscFunctionBegin; 984 PetscFunctionReturn(0); 985 } 986 987 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 988 { 989 PetscFunctionBegin; 990 PetscFunctionReturn(0); 991 } 992 993 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 994 { 995 PetscFunctionBegin; 996 PetscFunctionReturn(0); 997 } 998 999 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1000 { 1001 PetscFunctionBegin; 1002 PetscFunctionReturn(0); 1003 } 1004 1005 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 1006 { 1007 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1008 MPRKTableau tab = mprk->tableau; 1009 Vec YdotRHS_slow,YdotRHS_slowbuffer,YdotRHS_medium,YdotRHS_mediumbuffer,YdotRHS_fast; 1010 1011 PetscFunctionBegin; 1012 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y)); 1013 if (mprk->is_slow) { 1014 PetscCall(PetscMalloc1(tab->s,&mprk->work_slow)); 1015 } 1016 PetscCall(PetscMalloc1(tab->s,&mprk->work_slowbuffer)); 1017 if (tab->np == 3) { 1018 if (mprk->is_medium) { 1019 PetscCall(PetscMalloc1(tab->s,&mprk->work_medium)); 1020 } 1021 PetscCall(PetscMalloc1(tab->s,&mprk->work_mediumbuffer)); 1022 } 1023 PetscCall(PetscMalloc1(tab->s,&mprk->work_fast)); 1024 1025 if (ts->use_splitrhsfunction) { 1026 if (mprk->is_slow) { 1027 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow)); 1028 PetscCall(VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow)); 1029 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow)); 1030 } 1031 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer)); 1032 PetscCall(VecDuplicateVecs(YdotRHS_slowbuffer,tab->s,&mprk->YdotRHS_slowbuffer)); 1033 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_slowbuffer,&YdotRHS_slowbuffer)); 1034 if (tab->np == 3) { 1035 if (mprk->is_medium) { 1036 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium)); 1037 PetscCall(VecDuplicateVecs(YdotRHS_medium,tab->s,&mprk->YdotRHS_medium)); 1038 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_medium,&YdotRHS_medium)); 1039 } 1040 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer)); 1041 PetscCall(VecDuplicateVecs(YdotRHS_mediumbuffer,tab->s,&mprk->YdotRHS_mediumbuffer)); 1042 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_mediumbuffer,&YdotRHS_mediumbuffer)); 1043 } 1044 PetscCall(VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast)); 1045 PetscCall(VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast)); 1046 PetscCall(VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast)); 1047 } else { 1048 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS)); 1049 if (mprk->is_slow) { 1050 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_slow)); 1051 } 1052 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_slowbuffer)); 1053 if (tab->np == 3) { 1054 if (mprk->is_medium) { 1055 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_medium)); 1056 } 1057 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_mediumbuffer)); 1058 } 1059 PetscCall(PetscMalloc1(tab->s,&mprk->YdotRHS_fast)); 1060 } 1061 PetscFunctionReturn(0); 1062 } 1063 1064 static PetscErrorCode TSSetUp_MPRK(TS ts) 1065 { 1066 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1067 MPRKTableau tab = mprk->tableau; 1068 DM dm; 1069 1070 PetscFunctionBegin; 1071 PetscCall(TSRHSSplitGetIS(ts,"slow",&mprk->is_slow)); 1072 PetscCall(TSRHSSplitGetIS(ts,"fast",&mprk->is_fast)); 1073 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); 1074 1075 if (tab->np == 3) { 1076 PetscCall(TSRHSSplitGetIS(ts,"medium",&mprk->is_medium)); 1077 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); 1078 PetscCall(TSRHSSplitGetIS(ts,"mediumbuffer",&mprk->is_mediumbuffer)); 1079 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1080 mprk->is_mediumbuffer = mprk->is_medium; 1081 mprk->is_medium = NULL; 1082 } 1083 } 1084 1085 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1086 PetscCall(TSRHSSplitGetIS(ts,"slowbuffer",&mprk->is_slowbuffer)); 1087 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1088 mprk->is_slowbuffer = mprk->is_slow; 1089 mprk->is_slow = NULL; 1090 } 1091 PetscCall(TSCheckImplicitTerm(ts)); 1092 PetscCall(TSMPRKTableauSetUp(ts)); 1093 PetscCall(TSGetDM(ts,&dm)); 1094 PetscCall(DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts)); 1095 PetscCall(DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts)); 1096 if (ts->use_splitrhsfunction) { 1097 ts->ops->step = TSStep_MPRKSPLIT; 1098 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1099 PetscCall(TSMPRKSetSplits(ts)); 1100 } else { 1101 ts->ops->step = TSStep_MPRK; 1102 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1103 } 1104 PetscFunctionReturn(0); 1105 } 1106 1107 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 1108 { 1109 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1110 1111 PetscFunctionBegin; 1112 PetscOptionsHeadBegin(PetscOptionsObject,"PRK ODE solver options"); 1113 { 1114 MPRKTableauLink link; 1115 PetscInt count,choice; 1116 PetscBool flg; 1117 const char **namelist; 1118 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 1119 PetscCall(PetscMalloc1(count,(char***)&namelist)); 1120 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1121 PetscCall(PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg)); 1122 if (flg) PetscCall(TSMPRKSetType(ts,namelist[choice])); 1123 PetscCall(PetscFree(namelist)); 1124 } 1125 PetscOptionsHeadEnd(); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 1130 { 1131 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1132 PetscBool iascii; 1133 1134 PetscFunctionBegin; 1135 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1136 if (iascii) { 1137 MPRKTableau tab = mprk->tableau; 1138 TSMPRKType mprktype; 1139 char fbuf[512]; 1140 char sbuf[512]; 1141 PetscInt i; 1142 PetscCall(TSMPRKGetType(ts,&mprktype)); 1143 PetscCall(PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype)); 1144 PetscCall(PetscViewerASCIIPrintf(viewer," Order: %" PetscInt_FMT "\n",tab->order)); 1145 1146 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf)); 1147 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf)); 1148 PetscCall(PetscViewerASCIIPrintf(viewer," Af = \n")); 1149 for (i=0; i<tab->s; i++) { 1150 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,&tab->Af[i*tab->s])); 1151 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",fbuf)); 1152 } 1153 PetscCall(PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->bf)); 1154 PetscCall(PetscViewerASCIIPrintf(viewer," bf = %s\n",fbuf)); 1155 1156 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->csb)); 1157 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa csb = %s\n",sbuf)); 1158 PetscCall(PetscViewerASCIIPrintf(viewer," Asb = \n")); 1159 for (i=0; i<tab->s; i++) { 1160 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,&tab->Asb[i*tab->s])); 1161 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",sbuf)); 1162 } 1163 PetscCall(PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->bsb)); 1164 PetscCall(PetscViewerASCIIPrintf(viewer," bsb = %s\n",sbuf)); 1165 1166 if (tab->np == 3) { 1167 char mbuf[512]; 1168 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->cmb)); 1169 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa cmb = %s\n",mbuf)); 1170 PetscCall(PetscViewerASCIIPrintf(viewer," Amb = \n")); 1171 for (i=0; i<tab->s; i++) { 1172 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,&tab->Amb[i*tab->s])); 1173 PetscCall(PetscViewerASCIIPrintf(viewer," %s\n",mbuf)); 1174 } 1175 PetscCall(PetscFormatRealArray(mbuf,sizeof(mbuf),"% 8.6f",tab->s,tab->bmb)); 1176 PetscCall(PetscViewerASCIIPrintf(viewer," bmb = %s\n",mbuf)); 1177 } 1178 } 1179 PetscFunctionReturn(0); 1180 } 1181 1182 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 1183 { 1184 TSAdapt adapt; 1185 1186 PetscFunctionBegin; 1187 PetscCall(TSGetAdapt(ts,&adapt)); 1188 PetscCall(TSAdaptLoad(adapt,viewer)); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 /*@C 1193 TSMPRKSetType - Set the type of MPRK scheme 1194 1195 Not collective 1196 1197 Input Parameters: 1198 + ts - timestepping context 1199 - mprktype - type of MPRK-scheme 1200 1201 Options Database: 1202 . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme 1203 1204 Level: intermediate 1205 1206 .seealso: `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType` 1207 @*/ 1208 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 1209 { 1210 PetscFunctionBegin; 1211 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1212 PetscValidCharPointer(mprktype,2); 1213 PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype)); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /*@C 1218 TSMPRKGetType - Get the type of MPRK scheme 1219 1220 Not collective 1221 1222 Input Parameter: 1223 . ts - timestepping context 1224 1225 Output Parameter: 1226 . mprktype - type of MPRK-scheme 1227 1228 Level: intermediate 1229 1230 .seealso: `TSMPRKGetType()` 1231 @*/ 1232 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 1233 { 1234 PetscFunctionBegin; 1235 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1236 PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype)); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 1241 { 1242 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1243 1244 PetscFunctionBegin; 1245 *mprktype = mprk->tableau->name; 1246 PetscFunctionReturn(0); 1247 } 1248 1249 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 1250 { 1251 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1252 PetscBool match; 1253 MPRKTableauLink link; 1254 1255 PetscFunctionBegin; 1256 if (mprk->tableau) { 1257 PetscCall(PetscStrcmp(mprk->tableau->name,mprktype,&match)); 1258 if (match) PetscFunctionReturn(0); 1259 } 1260 for (link = MPRKTableauList; link; link=link->next) { 1261 PetscCall(PetscStrcmp(link->tab.name,mprktype,&match)); 1262 if (match) { 1263 if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts)); 1264 mprk->tableau = &link->tab; 1265 if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts)); 1266 PetscFunctionReturn(0); 1267 } 1268 } 1269 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 1270 } 1271 1272 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 1273 { 1274 TS_MPRK *mprk = (TS_MPRK*)ts->data; 1275 1276 PetscFunctionBegin; 1277 *ns = mprk->tableau->s; 1278 if (Y) *Y = mprk->Y; 1279 PetscFunctionReturn(0); 1280 } 1281 1282 static PetscErrorCode TSDestroy_MPRK(TS ts) 1283 { 1284 PetscFunctionBegin; 1285 PetscCall(TSReset_MPRK(ts)); 1286 if (ts->dm) { 1287 PetscCall(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts)); 1288 PetscCall(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts)); 1289 } 1290 PetscCall(PetscFree(ts->data)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL)); 1293 PetscFunctionReturn(0); 1294 } 1295 1296 /*MC 1297 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes 1298 1299 The user should provide the right hand side of the equation 1300 using TSSetRHSFunction(). 1301 1302 Notes: 1303 The default is TSMPRKPM2, it can be changed with TSMPRKSetType() or -ts_mprk_type 1304 1305 Level: beginner 1306 1307 .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()` 1308 `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2` 1309 1310 M*/ 1311 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1312 { 1313 TS_MPRK *mprk; 1314 1315 PetscFunctionBegin; 1316 PetscCall(TSMPRKInitializePackage()); 1317 1318 ts->ops->reset = TSReset_MPRK; 1319 ts->ops->destroy = TSDestroy_MPRK; 1320 ts->ops->view = TSView_MPRK; 1321 ts->ops->load = TSLoad_MPRK; 1322 ts->ops->setup = TSSetUp_MPRK; 1323 ts->ops->step = TSStep_MPRK; 1324 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1325 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1326 ts->ops->getstages = TSGetStages_MPRK; 1327 1328 PetscCall(PetscNewLog(ts,&mprk)); 1329 ts->data = (void*)mprk; 1330 1331 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK)); 1332 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK)); 1333 1334 PetscCall(TSMPRKSetType(ts,TSMPRKDefault)); 1335 PetscFunctionReturn(0); 1336 } 1337