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