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, No Fortran Support 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, PetscCtx 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, PetscCtx ctx) 946 { 947 PetscFunctionBegin; 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm, DM subdm, PetscCtx ctx) 952 { 953 PetscFunctionBegin; 954 PetscFunctionReturn(PETSC_SUCCESS); 955 } 956 957 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx 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 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1070 for (link = MPRKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1071 PetscCall(PetscOptionsEList("-ts_mprk_type", "Family of MPRK method", "TSMPRKSetType", (const char *const *)namelist, count, mprk->tableau->name, &choice, &flg)); 1072 if (flg) PetscCall(TSMPRKSetType(ts, namelist[choice])); 1073 PetscCall(PetscFree(namelist)); 1074 } 1075 PetscOptionsHeadEnd(); 1076 PetscFunctionReturn(PETSC_SUCCESS); 1077 } 1078 1079 static PetscErrorCode TSView_MPRK(TS ts, PetscViewer viewer) 1080 { 1081 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1082 PetscBool isascii; 1083 1084 PetscFunctionBegin; 1085 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 1086 if (isascii) { 1087 MPRKTableau tab = mprk->tableau; 1088 TSMPRKType mprktype; 1089 char fbuf[512]; 1090 char sbuf[512]; 1091 PetscInt i; 1092 PetscCall(TSMPRKGetType(ts, &mprktype)); 1093 PetscCall(PetscViewerASCIIPrintf(viewer, " MPRK type %s\n", mprktype)); 1094 PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order)); 1095 1096 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->cf)); 1097 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cf = %s\n", fbuf)); 1098 PetscCall(PetscViewerASCIIPrintf(viewer, " Af = \n")); 1099 for (i = 0; i < tab->s; i++) { 1100 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, &tab->Af[i * tab->s])); 1101 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", fbuf)); 1102 } 1103 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->bf)); 1104 PetscCall(PetscViewerASCIIPrintf(viewer, " bf = %s\n", fbuf)); 1105 1106 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->csb)); 1107 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa csb = %s\n", sbuf)); 1108 PetscCall(PetscViewerASCIIPrintf(viewer, " Asb = \n")); 1109 for (i = 0; i < tab->s; i++) { 1110 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, &tab->Asb[i * tab->s])); 1111 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", sbuf)); 1112 } 1113 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->bsb)); 1114 PetscCall(PetscViewerASCIIPrintf(viewer, " bsb = %s\n", sbuf)); 1115 1116 if (tab->np == 3) { 1117 char mbuf[512]; 1118 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->cmb)); 1119 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cmb = %s\n", mbuf)); 1120 PetscCall(PetscViewerASCIIPrintf(viewer, " Amb = \n")); 1121 for (i = 0; i < tab->s; i++) { 1122 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, &tab->Amb[i * tab->s])); 1123 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", mbuf)); 1124 } 1125 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->bmb)); 1126 PetscCall(PetscViewerASCIIPrintf(viewer, " bmb = %s\n", mbuf)); 1127 } 1128 } 1129 PetscFunctionReturn(PETSC_SUCCESS); 1130 } 1131 1132 static PetscErrorCode TSLoad_MPRK(TS ts, PetscViewer viewer) 1133 { 1134 TSAdapt adapt; 1135 1136 PetscFunctionBegin; 1137 PetscCall(TSGetAdapt(ts, &adapt)); 1138 PetscCall(TSAdaptLoad(adapt, viewer)); 1139 PetscFunctionReturn(PETSC_SUCCESS); 1140 } 1141 1142 /*@ 1143 TSMPRKSetType - Set the type of `TSMPRK` scheme 1144 1145 Not Collective 1146 1147 Input Parameters: 1148 + ts - timestepping context 1149 - mprktype - type of `TSMPRK` scheme 1150 1151 Options Database Key: 1152 . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme 1153 1154 Level: intermediate 1155 1156 .seealso: [](ch_ts), `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType` 1157 @*/ 1158 PetscErrorCode TSMPRKSetType(TS ts, TSMPRKType mprktype) 1159 { 1160 PetscFunctionBegin; 1161 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1162 PetscAssertPointer(mprktype, 2); 1163 PetscTryMethod(ts, "TSMPRKSetType_C", (TS, TSMPRKType), (ts, mprktype)); 1164 PetscFunctionReturn(PETSC_SUCCESS); 1165 } 1166 1167 /*@ 1168 TSMPRKGetType - Get the type of `TSMPRK` scheme 1169 1170 Not Collective 1171 1172 Input Parameter: 1173 . ts - timestepping context 1174 1175 Output Parameter: 1176 . mprktype - type of `TSMPRK` scheme 1177 1178 Level: intermediate 1179 1180 .seealso: [](ch_ts), `TSMPRK` 1181 @*/ 1182 PetscErrorCode TSMPRKGetType(TS ts, TSMPRKType *mprktype) 1183 { 1184 PetscFunctionBegin; 1185 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1186 PetscUseMethod(ts, "TSMPRKGetType_C", (TS, TSMPRKType *), (ts, mprktype)); 1187 PetscFunctionReturn(PETSC_SUCCESS); 1188 } 1189 1190 static PetscErrorCode TSMPRKGetType_MPRK(TS ts, TSMPRKType *mprktype) 1191 { 1192 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1193 1194 PetscFunctionBegin; 1195 *mprktype = mprk->tableau->name; 1196 PetscFunctionReturn(PETSC_SUCCESS); 1197 } 1198 1199 static PetscErrorCode TSMPRKSetType_MPRK(TS ts, TSMPRKType mprktype) 1200 { 1201 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1202 PetscBool match; 1203 MPRKTableauLink link; 1204 1205 PetscFunctionBegin; 1206 if (mprk->tableau) { 1207 PetscCall(PetscStrcmp(mprk->tableau->name, mprktype, &match)); 1208 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1209 } 1210 for (link = MPRKTableauList; link; link = link->next) { 1211 PetscCall(PetscStrcmp(link->tab.name, mprktype, &match)); 1212 if (match) { 1213 if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts)); 1214 mprk->tableau = &link->tab; 1215 if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts)); 1216 PetscFunctionReturn(PETSC_SUCCESS); 1217 } 1218 } 1219 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mprktype); 1220 } 1221 1222 static PetscErrorCode TSGetStages_MPRK(TS ts, PetscInt *ns, Vec **Y) 1223 { 1224 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1225 1226 PetscFunctionBegin; 1227 *ns = mprk->tableau->s; 1228 if (Y) *Y = mprk->Y; 1229 PetscFunctionReturn(PETSC_SUCCESS); 1230 } 1231 1232 static PetscErrorCode TSDestroy_MPRK(TS ts) 1233 { 1234 PetscFunctionBegin; 1235 PetscCall(TSReset_MPRK(ts)); 1236 if (ts->dm) { 1237 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts)); 1238 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts)); 1239 } 1240 PetscCall(PetscFree(ts->data)); 1241 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", NULL)); 1242 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", NULL)); 1243 PetscFunctionReturn(PETSC_SUCCESS); 1244 } 1245 1246 /*MC 1247 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes 1248 1249 The user should provide the right-hand side of the equation using `TSSetRHSFunction()`. 1250 1251 Level: beginner 1252 1253 Note: 1254 The default is `TSMPRKPM2`, it can be changed with `TSMPRKSetType()` or -ts_mprk_type 1255 1256 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()` 1257 `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2`, `TSType` 1258 M*/ 1259 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1260 { 1261 TS_MPRK *mprk; 1262 1263 PetscFunctionBegin; 1264 PetscCall(TSMPRKInitializePackage()); 1265 1266 ts->ops->reset = TSReset_MPRK; 1267 ts->ops->destroy = TSDestroy_MPRK; 1268 ts->ops->view = TSView_MPRK; 1269 ts->ops->load = TSLoad_MPRK; 1270 ts->ops->setup = TSSetUp_MPRK; 1271 ts->ops->step = TSStep_MPRK; 1272 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1273 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1274 ts->ops->getstages = TSGetStages_MPRK; 1275 1276 PetscCall(PetscNew(&mprk)); 1277 ts->data = (void *)mprk; 1278 1279 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", TSMPRKGetType_MPRK)); 1280 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", TSMPRKSetType_MPRK)); 1281 1282 PetscCall(TSMPRKSetType(ts, TSMPRKDefault)); 1283 PetscFunctionReturn(PETSC_SUCCESS); 1284 } 1285