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