1 /* 2 Code for time stepping with the Multirate Partitioned Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) 7 if one does not split the RHS function, but gives the indexes for both slow and fast components; 8 2) The general system is written as 9 Usdot = Fs(t,Us,Uf) 10 Ufdot = Ff(t,Us,Uf) 11 for component-wise partitioned system, 12 users should split the RHS function themselves and also provide the indexes for both slow and fast components. 13 3) To correct The confusing terminology in the paper, we use 'slow method', 'slow buffer method' and 'fast method' to denote the methods applied to 'slow region', 'slow buffer region' and 'fast region' respectively. The 'slow method' in the original paper actually means the 'slow buffer method'. 14 4) Why does the buffer region have to be inside the slow region? The buffer region is treated with a slow method essentially. Applying the slow method to a region with a fast characteristic time scale is apparently not a good choice. 15 16 Reference: 17 Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007 18 */ 19 20 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 21 #include <petscdm.h> 22 23 static TSMPRKType TSMPRKDefault = TSMPRK2A22; 24 static PetscBool TSMPRKRegisterAllCalled; 25 static PetscBool TSMPRKPackageInitialized; 26 27 typedef struct _MPRKTableau *MPRKTableau; 28 struct _MPRKTableau { 29 char *name; 30 PetscInt order; /* Classical approximation order of the method i */ 31 PetscInt sbase; /* Number of stages in the base method*/ 32 PetscInt s; /* Number of stages */ 33 PetscInt np; /* Number of partitions */ 34 PetscReal *Af, *bf, *cf; /* Tableau for fast components */ 35 PetscReal *Amb, *bmb, *cmb; /* Tableau for medium components */ 36 PetscInt *rmb; /* Array of flags for repeated stages in medium method */ 37 PetscReal *Asb, *bsb, *csb; /* Tableau for slow components */ 38 PetscInt *rsb; /* Array of flags for repeated staged in slow method*/ 39 }; 40 typedef struct _MPRKTableauLink *MPRKTableauLink; 41 struct _MPRKTableauLink { 42 struct _MPRKTableau tab; 43 MPRKTableauLink next; 44 }; 45 static MPRKTableauLink MPRKTableauList; 46 47 typedef struct { 48 MPRKTableau tableau; 49 Vec *Y; /* States computed during the step */ 50 Vec *YdotRHS; 51 Vec *YdotRHS_slow; /* Function evaluations by slow tableau for slow components */ 52 Vec *YdotRHS_slowbuffer; /* Function evaluations by slow tableau for slow components */ 53 Vec *YdotRHS_medium; /* Function evaluations by slow tableau for slow components */ 54 Vec *YdotRHS_mediumbuffer; /* Function evaluations by slow tableau for slow components */ 55 Vec *YdotRHS_fast; /* Function evaluations by fast tableau for fast components */ 56 PetscScalar *work_slow; /* Scalar work_slow by slow tableau */ 57 PetscScalar *work_slowbuffer; /* Scalar work_slow by slow tableau */ 58 PetscScalar *work_medium; /* Scalar work_slow by medium tableau */ 59 PetscScalar *work_mediumbuffer; /* Scalar work_slow by medium tableau */ 60 PetscScalar *work_fast; /* Scalar work_fast by fast tableau */ 61 PetscReal stage_time; 62 TSStepStatus status; 63 PetscReal ptime; 64 PetscReal time_step; 65 IS is_slow, is_slowbuffer, is_medium, is_mediumbuffer, is_fast; 66 TS subts_slow, subts_slowbuffer, subts_medium, subts_mediumbuffer, subts_fast; 67 } TS_MPRK; 68 69 static PetscErrorCode TSMPRKGenerateTableau2(PetscInt ratio, PetscInt s, const PetscReal Abase[], const PetscReal bbase[], PetscReal A1[], PetscReal b1[], PetscReal A2[], PetscReal b2[]) 70 { 71 PetscInt i, j, k, l; 72 73 PetscFunctionBegin; 74 for (k = 0; k < ratio; k++) { 75 /* diagonal blocks */ 76 for (i = 0; i < s; i++) 77 for (j = 0; j < s; j++) { 78 A1[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j]; 79 A2[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j] / ratio; 80 } 81 /* off diagonal blocks */ 82 for (l = 0; l < k; l++) 83 for (i = 0; i < s; i++) 84 for (j = 0; j < s; j++) 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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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: [](chapter_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 . s - 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 . Af - stage coefficients for fast components(dimension s*s, row-major) 446 . bf - step completion table for fast components(dimension s) 447 . cf - abscissa for fast components(dimension s) 448 . As - stage coefficients for slow components(dimension s*s, row-major) 449 . bs - step completion table for slow components(dimension s) 450 - cs - abscissa for slow components(dimension s) 451 452 Level: advanced 453 454 Note: 455 Several `TSMPRK` methods are provided, this function is only needed to create new methods. 456 457 .seealso: [](chapter_ts), `TSMPRK` 458 @*/ 459 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[]) 460 { 461 MPRKTableauLink link; 462 MPRKTableau t; 463 PetscInt s, i, j; 464 465 PetscFunctionBegin; 466 PetscValidCharPointer(name, 1); 467 PetscValidRealPointer(Asb, 6); 468 if (bsb) PetscValidRealPointer(bsb, 7); 469 if (csb) PetscValidRealPointer(csb, 8); 470 if (rsb) PetscValidIntPointer(rsb, 9); 471 if (Amb) PetscValidRealPointer(Amb, 10); 472 if (bmb) PetscValidRealPointer(bmb, 11); 473 if (cmb) PetscValidRealPointer(cmb, 12); 474 if (rmb) PetscValidIntPointer(rmb, 13); 475 PetscValidRealPointer(Af, 14); 476 if (bf) PetscValidRealPointer(bf, 15); 477 if (cf) PetscValidRealPointer(cf, 16); 478 479 PetscCall(PetscNew(&link)); 480 t = &link->tab; 481 482 PetscCall(PetscStrallocpy(name, &t->name)); 483 s = sbase * ratio1 * ratio2; /* this is the dimension of the matrices below */ 484 t->order = order; 485 t->sbase = sbase; 486 t->s = s; 487 t->np = 2; 488 489 PetscCall(PetscMalloc3(s * s, &t->Af, s, &t->bf, s, &t->cf)); 490 PetscCall(PetscArraycpy(t->Af, Af, s * s)); 491 if (bf) { 492 PetscCall(PetscArraycpy(t->bf, bf, s)); 493 } else 494 for (i = 0; i < s; i++) t->bf[i] = Af[(s - 1) * s + i]; 495 if (cf) { 496 PetscCall(PetscArraycpy(t->cf, cf, s)); 497 } else { 498 for (i = 0; i < s; i++) 499 for (j = 0, t->cf[i] = 0; j < s; j++) t->cf[i] += Af[i * s + j]; 500 } 501 502 if (Amb) { 503 t->np = 3; 504 PetscCall(PetscMalloc3(s * s, &t->Amb, s, &t->bmb, s, &t->cmb)); 505 PetscCall(PetscArraycpy(t->Amb, Amb, s * s)); 506 if (bmb) { 507 PetscCall(PetscArraycpy(t->bmb, bmb, s)); 508 } else { 509 for (i = 0; i < s; i++) t->bmb[i] = Amb[(s - 1) * s + i]; 510 } 511 if (cmb) { 512 PetscCall(PetscArraycpy(t->cmb, cmb, s)); 513 } else { 514 for (i = 0; i < s; i++) 515 for (j = 0, t->cmb[i] = 0; j < s; j++) t->cmb[i] += Amb[i * s + j]; 516 } 517 if (rmb) { 518 PetscCall(PetscMalloc1(s, &t->rmb)); 519 PetscCall(PetscArraycpy(t->rmb, rmb, s)); 520 } else { 521 PetscCall(PetscCalloc1(s, &t->rmb)); 522 } 523 } 524 525 PetscCall(PetscMalloc3(s * s, &t->Asb, s, &t->bsb, s, &t->csb)); 526 PetscCall(PetscArraycpy(t->Asb, Asb, s * s)); 527 if (bsb) { 528 PetscCall(PetscArraycpy(t->bsb, bsb, s)); 529 } else 530 for (i = 0; i < s; i++) t->bsb[i] = Asb[(s - 1) * s + i]; 531 if (csb) { 532 PetscCall(PetscArraycpy(t->csb, csb, s)); 533 } else { 534 for (i = 0; i < s; i++) 535 for (j = 0, t->csb[i] = 0; j < s; j++) t->csb[i] += Asb[i * s + j]; 536 } 537 if (rsb) { 538 PetscCall(PetscMalloc1(s, &t->rsb)); 539 PetscCall(PetscArraycpy(t->rsb, rsb, s)); 540 } else { 541 PetscCall(PetscCalloc1(s, &t->rsb)); 542 } 543 link->next = MPRKTableauList; 544 MPRKTableauList = link; 545 PetscFunctionReturn(PETSC_SUCCESS); 546 } 547 548 static PetscErrorCode TSMPRKSetSplits(TS ts) 549 { 550 TS_MPRK *mprk = (TS_MPRK *)ts->data; 551 MPRKTableau tab = mprk->tableau; 552 DM dm, subdm, newdm; 553 554 PetscFunctionBegin; 555 PetscCall(TSRHSSplitGetSubTS(ts, "slow", &mprk->subts_slow)); 556 PetscCall(TSRHSSplitGetSubTS(ts, "fast", &mprk->subts_fast)); 557 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"); 558 559 /* Only copy the DM */ 560 PetscCall(TSGetDM(ts, &dm)); 561 562 PetscCall(TSRHSSplitGetSubTS(ts, "slowbuffer", &mprk->subts_slowbuffer)); 563 if (!mprk->subts_slowbuffer) { 564 mprk->subts_slowbuffer = mprk->subts_slow; 565 mprk->subts_slow = NULL; 566 } 567 if (mprk->subts_slow) { 568 PetscCall(DMClone(dm, &newdm)); 569 PetscCall(TSGetDM(mprk->subts_slow, &subdm)); 570 PetscCall(DMCopyDMTS(subdm, newdm)); 571 PetscCall(DMCopyDMSNES(subdm, newdm)); 572 PetscCall(TSSetDM(mprk->subts_slow, newdm)); 573 PetscCall(DMDestroy(&newdm)); 574 } 575 PetscCall(DMClone(dm, &newdm)); 576 PetscCall(TSGetDM(mprk->subts_slowbuffer, &subdm)); 577 PetscCall(DMCopyDMTS(subdm, newdm)); 578 PetscCall(DMCopyDMSNES(subdm, newdm)); 579 PetscCall(TSSetDM(mprk->subts_slowbuffer, newdm)); 580 PetscCall(DMDestroy(&newdm)); 581 582 PetscCall(DMClone(dm, &newdm)); 583 PetscCall(TSGetDM(mprk->subts_fast, &subdm)); 584 PetscCall(DMCopyDMTS(subdm, newdm)); 585 PetscCall(DMCopyDMSNES(subdm, newdm)); 586 PetscCall(TSSetDM(mprk->subts_fast, newdm)); 587 PetscCall(DMDestroy(&newdm)); 588 589 if (tab->np == 3) { 590 PetscCall(TSRHSSplitGetSubTS(ts, "medium", &mprk->subts_medium)); 591 PetscCall(TSRHSSplitGetSubTS(ts, "mediumbuffer", &mprk->subts_mediumbuffer)); 592 if (mprk->subts_medium && !mprk->subts_mediumbuffer) { 593 mprk->subts_mediumbuffer = mprk->subts_medium; 594 mprk->subts_medium = NULL; 595 } 596 if (mprk->subts_medium) { 597 PetscCall(DMClone(dm, &newdm)); 598 PetscCall(TSGetDM(mprk->subts_medium, &subdm)); 599 PetscCall(DMCopyDMTS(subdm, newdm)); 600 PetscCall(DMCopyDMSNES(subdm, newdm)); 601 PetscCall(TSSetDM(mprk->subts_medium, newdm)); 602 PetscCall(DMDestroy(&newdm)); 603 } 604 PetscCall(DMClone(dm, &newdm)); 605 PetscCall(TSGetDM(mprk->subts_mediumbuffer, &subdm)); 606 PetscCall(DMCopyDMTS(subdm, newdm)); 607 PetscCall(DMCopyDMSNES(subdm, newdm)); 608 PetscCall(TSSetDM(mprk->subts_mediumbuffer, newdm)); 609 PetscCall(DMDestroy(&newdm)); 610 } 611 PetscFunctionReturn(PETSC_SUCCESS); 612 } 613 614 /* 615 This if for nonsplit RHS MPRK 616 The step completion formula is 617 618 x1 = x0 + h b^T YdotRHS 619 620 */ 621 static PetscErrorCode TSEvaluateStep_MPRK(TS ts, PetscInt order, Vec X, PetscBool *done) 622 { 623 TS_MPRK *mprk = (TS_MPRK *)ts->data; 624 MPRKTableau tab = mprk->tableau; 625 PetscScalar *wf = mprk->work_fast; 626 PetscReal h = ts->time_step; 627 PetscInt s = tab->s, j; 628 629 PetscFunctionBegin; 630 for (j = 0; j < s; j++) wf[j] = h * tab->bf[j]; 631 PetscCall(VecCopy(ts->vec_sol, X)); 632 PetscCall(VecMAXPY(X, s, wf, mprk->YdotRHS)); 633 PetscFunctionReturn(PETSC_SUCCESS); 634 } 635 636 static PetscErrorCode TSStep_MPRK(TS ts) 637 { 638 TS_MPRK *mprk = (TS_MPRK *)ts->data; 639 Vec *Y = mprk->Y, *YdotRHS = mprk->YdotRHS, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 640 Vec Yslow, Yslowbuffer, Yfast; 641 MPRKTableau tab = mprk->tableau; 642 const PetscInt s = tab->s; 643 const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb; 644 PetscScalar *wf = mprk->work_fast, *wsb = mprk->work_slowbuffer; 645 PetscInt i, j; 646 PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step; 647 648 PetscFunctionBegin; 649 for (i = 0; i < s; i++) { 650 mprk->stage_time = t + h * cf[i]; 651 PetscCall(TSPreStage(ts, mprk->stage_time)); 652 PetscCall(VecCopy(ts->vec_sol, Y[i])); 653 654 /* slow buffer region */ 655 for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j]; 656 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j])); 657 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 658 PetscCall(VecMAXPY(Yslowbuffer, i, wsb, mprk->YdotRHS_slowbuffer)); 659 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 660 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j])); 661 /* slow region */ 662 if (mprk->is_slow) { 663 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j])); 664 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow)); 665 PetscCall(VecMAXPY(Yslow, i, wsb, mprk->YdotRHS_slow)); 666 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow)); 667 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j])); 668 } 669 670 /* fast region */ 671 for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j]; 672 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j])); 673 PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast)); 674 PetscCall(VecMAXPY(Yfast, i, wf, mprk->YdotRHS_fast)); 675 PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast)); 676 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j])); 677 if (tab->np == 3) { 678 Vec *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 679 Vec Ymedium, Ymediumbuffer; 680 PetscScalar *wmb = mprk->work_mediumbuffer; 681 682 for (j = 0; j < i; j++) wmb[j] = h * tab->Amb[i * s + j]; 683 /* medium region */ 684 if (mprk->is_medium) { 685 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j])); 686 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium)); 687 PetscCall(VecMAXPY(Ymedium, i, wmb, mprk->YdotRHS_medium)); 688 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium)); 689 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j])); 690 } 691 /* medium buffer region */ 692 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j])); 693 PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer)); 694 PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, mprk->YdotRHS_mediumbuffer)); 695 PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer)); 696 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j])); 697 } 698 PetscCall(TSPostStage(ts, mprk->stage_time, i, Y)); 699 /* compute the stage RHS by fast and slow tableau respectively */ 700 PetscCall(TSComputeRHSFunction(ts, t + h * csb[i], Y[i], YdotRHS[i])); 701 } 702 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL)); 703 ts->ptime += ts->time_step; 704 ts->time_step = next_time_step; 705 PetscFunctionReturn(PETSC_SUCCESS); 706 } 707 708 /* 709 This if for the case when split RHS is used 710 The step completion formula is 711 x1 = x0 + h b^T YdotRHS 712 */ 713 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts, PetscInt order, Vec X, PetscBool *done) 714 { 715 TS_MPRK *mprk = (TS_MPRK *)ts->data; 716 MPRKTableau tab = mprk->tableau; 717 Vec Xslow, Xfast, Xslowbuffer; /* subvectors for slow and fast components in X respectively */ 718 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer; 719 PetscReal h = ts->time_step; 720 PetscInt s = tab->s, j, computedstages; 721 722 PetscFunctionBegin; 723 PetscCall(VecCopy(ts->vec_sol, X)); 724 725 /* slow region */ 726 if (mprk->is_slow) { 727 computedstages = 0; 728 for (j = 0; j < s; j++) { 729 if (tab->rsb[j]) ws[tab->rsb[j] - 1] += h * tab->bsb[j]; 730 else ws[computedstages++] = h * tab->bsb[j]; 731 } 732 PetscCall(VecGetSubVector(X, mprk->is_slow, &Xslow)); 733 PetscCall(VecMAXPY(Xslow, computedstages, ws, mprk->YdotRHS_slow)); 734 PetscCall(VecRestoreSubVector(X, mprk->is_slow, &Xslow)); 735 } 736 737 if (tab->np == 3 && mprk->is_medium) { 738 computedstages = 0; 739 for (j = 0; j < s; j++) { 740 if (tab->rmb[j]) wsb[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bsb[j]; 741 else wsb[computedstages++] = h * tab->bsb[j]; 742 } 743 PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer)); 744 PetscCall(VecMAXPY(Xslowbuffer, computedstages, wsb, mprk->YdotRHS_slowbuffer)); 745 PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer)); 746 } else { 747 /* slow buffer region */ 748 for (j = 0; j < s; j++) wsb[j] = h * tab->bsb[j]; 749 PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer)); 750 PetscCall(VecMAXPY(Xslowbuffer, s, wsb, mprk->YdotRHS_slowbuffer)); 751 PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer)); 752 } 753 if (tab->np == 3) { 754 Vec Xmedium, Xmediumbuffer; 755 PetscScalar *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer; 756 /* medium region and slow buffer region */ 757 if (mprk->is_medium) { 758 computedstages = 0; 759 for (j = 0; j < s; j++) { 760 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bmb[j]; 761 else wm[computedstages++] = h * tab->bmb[j]; 762 } 763 PetscCall(VecGetSubVector(X, mprk->is_medium, &Xmedium)); 764 PetscCall(VecMAXPY(Xmedium, computedstages, wm, mprk->YdotRHS_medium)); 765 PetscCall(VecRestoreSubVector(X, mprk->is_medium, &Xmedium)); 766 } 767 /* medium buffer region */ 768 for (j = 0; j < s; j++) wmb[j] = h * tab->bmb[j]; 769 PetscCall(VecGetSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer)); 770 PetscCall(VecMAXPY(Xmediumbuffer, s, wmb, mprk->YdotRHS_mediumbuffer)); 771 PetscCall(VecRestoreSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer)); 772 } 773 /* fast region */ 774 for (j = 0; j < s; j++) wf[j] = h * tab->bf[j]; 775 PetscCall(VecGetSubVector(X, mprk->is_fast, &Xfast)); 776 PetscCall(VecMAXPY(Xfast, s, wf, mprk->YdotRHS_fast)); 777 PetscCall(VecRestoreSubVector(X, mprk->is_fast, &Xfast)); 778 PetscFunctionReturn(PETSC_SUCCESS); 779 } 780 781 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 782 { 783 TS_MPRK *mprk = (TS_MPRK *)ts->data; 784 MPRKTableau tab = mprk->tableau; 785 Vec *Y = mprk->Y, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer; 786 Vec Yslow, Yslowbuffer, Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 787 PetscInt s = tab->s; 788 const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb; 789 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer; 790 PetscInt i, j, computedstages; 791 PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step; 792 793 PetscFunctionBegin; 794 for (i = 0; i < s; i++) { 795 mprk->stage_time = t + h * cf[i]; 796 PetscCall(TSPreStage(ts, mprk->stage_time)); 797 /* calculate the stage value for fast and slow components respectively */ 798 PetscCall(VecCopy(ts->vec_sol, Y[i])); 799 for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j]; 800 801 /* slow buffer region */ 802 if (tab->np == 3 && mprk->is_medium) { 803 if (tab->rmb[i]) { 804 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 805 PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_slowbuffer, SCATTER_REVERSE, Yslowbuffer)); 806 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 807 } else { 808 PetscScalar *wm = mprk->work_medium; 809 computedstages = 0; 810 for (j = 0; j < i; j++) { 811 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wsb[j]; 812 else wm[computedstages++] = wsb[j]; 813 } 814 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 815 PetscCall(VecMAXPY(Yslowbuffer, computedstages, wm, YdotRHS_slowbuffer)); 816 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 817 } 818 } else { 819 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 820 PetscCall(VecMAXPY(Yslowbuffer, i, wsb, YdotRHS_slowbuffer)); 821 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer)); 822 } 823 824 /* slow region */ 825 if (mprk->is_slow) { 826 if (tab->rsb[i]) { /* repeat previous stage */ 827 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow)); 828 PetscCall(VecISCopy(Y[tab->rsb[i] - 1], mprk->is_slow, SCATTER_REVERSE, Yslow)); 829 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow)); 830 } else { 831 computedstages = 0; 832 for (j = 0; j < i; j++) { 833 if (tab->rsb[j]) ws[tab->rsb[j] - 1] += wsb[j]; 834 else ws[computedstages++] = wsb[j]; 835 } 836 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow)); 837 PetscCall(VecMAXPY(Yslow, computedstages, ws, YdotRHS_slow)); 838 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow)); 839 /* only depends on the slow buffer region */ 840 PetscCall(TSComputeRHSFunction(mprk->subts_slow, t + h * csb[i], Y[i], YdotRHS_slow[computedstages])); 841 } 842 } 843 844 /* fast region */ 845 for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j]; 846 PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast)); 847 PetscCall(VecMAXPY(Yfast, i, wf, YdotRHS_fast)); 848 PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast)); 849 850 if (tab->np == 3) { 851 Vec *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer; 852 Vec Ymedium, Ymediumbuffer; 853 const PetscReal *Amb = tab->Amb, *cmb = tab->cmb; 854 PetscScalar *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer; 855 856 for (j = 0; j < i; j++) wmb[j] = h * Amb[i * s + j]; 857 /* medium buffer region */ 858 PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer)); 859 PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, YdotRHS_mediumbuffer)); 860 PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer)); 861 /* medium region */ 862 if (mprk->is_medium) { 863 if (tab->rmb[i]) { /* repeat previous stage */ 864 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium)); 865 PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_medium, SCATTER_REVERSE, Ymedium)); 866 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium)); 867 } else { 868 computedstages = 0; 869 for (j = 0; j < i; j++) { 870 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wmb[j]; 871 else wm[computedstages++] = wmb[j]; 872 } 873 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium)); 874 PetscCall(VecMAXPY(Ymedium, computedstages, wm, YdotRHS_medium)); 875 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium)); 876 /* only depends on the medium buffer region */ 877 PetscCall(TSComputeRHSFunction(mprk->subts_medium, t + h * cmb[i], Y[i], YdotRHS_medium[computedstages])); 878 /* must be computed after all regions are updated in Y */ 879 PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[computedstages])); 880 } 881 } 882 /* must be computed after fast region and slow region are updated in Y */ 883 PetscCall(TSComputeRHSFunction(mprk->subts_mediumbuffer, t + h * cmb[i], Y[i], YdotRHS_mediumbuffer[i])); 884 } 885 if (!(tab->np == 3 && mprk->is_medium)) PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[i])); 886 PetscCall(TSComputeRHSFunction(mprk->subts_fast, t + h * cf[i], Y[i], YdotRHS_fast[i])); 887 } 888 889 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL)); 890 ts->ptime += ts->time_step; 891 ts->time_step = next_time_step; 892 PetscFunctionReturn(PETSC_SUCCESS); 893 } 894 895 static PetscErrorCode TSMPRKTableauReset(TS ts) 896 { 897 TS_MPRK *mprk = (TS_MPRK *)ts->data; 898 MPRKTableau tab = mprk->tableau; 899 900 PetscFunctionBegin; 901 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 902 PetscCall(PetscFree(mprk->work_fast)); 903 PetscCall(PetscFree(mprk->work_slow)); 904 PetscCall(PetscFree(mprk->work_slowbuffer)); 905 PetscCall(PetscFree(mprk->work_medium)); 906 PetscCall(PetscFree(mprk->work_mediumbuffer)); 907 PetscCall(VecDestroyVecs(tab->s, &mprk->Y)); 908 if (ts->use_splitrhsfunction) { 909 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_fast)); 910 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slow)); 911 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slowbuffer)); 912 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_medium)); 913 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_mediumbuffer)); 914 } else { 915 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS)); 916 if (mprk->is_slow) PetscCall(PetscFree(mprk->YdotRHS_slow)); 917 PetscCall(PetscFree(mprk->YdotRHS_slowbuffer)); 918 if (tab->np == 3) { 919 if (mprk->is_medium) PetscCall(PetscFree(mprk->YdotRHS_medium)); 920 PetscCall(PetscFree(mprk->YdotRHS_mediumbuffer)); 921 } 922 PetscCall(PetscFree(mprk->YdotRHS_fast)); 923 } 924 PetscFunctionReturn(PETSC_SUCCESS); 925 } 926 927 static PetscErrorCode TSReset_MPRK(TS ts) 928 { 929 PetscFunctionBegin; 930 PetscCall(TSMPRKTableauReset(ts)); 931 PetscFunctionReturn(PETSC_SUCCESS); 932 } 933 934 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine, DM coarse, void *ctx) 935 { 936 PetscFunctionBegin; 937 PetscFunctionReturn(PETSC_SUCCESS); 938 } 939 940 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 941 { 942 PetscFunctionBegin; 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 946 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm, DM subdm, void *ctx) 947 { 948 PetscFunctionBegin; 949 PetscFunctionReturn(PETSC_SUCCESS); 950 } 951 952 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 953 { 954 PetscFunctionBegin; 955 PetscFunctionReturn(PETSC_SUCCESS); 956 } 957 958 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 959 { 960 TS_MPRK *mprk = (TS_MPRK *)ts->data; 961 MPRKTableau tab = mprk->tableau; 962 Vec YdotRHS_slow, YdotRHS_slowbuffer, YdotRHS_medium, YdotRHS_mediumbuffer, YdotRHS_fast; 963 964 PetscFunctionBegin; 965 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->Y)); 966 if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->work_slow)); 967 PetscCall(PetscMalloc1(tab->s, &mprk->work_slowbuffer)); 968 if (tab->np == 3) { 969 if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->work_medium)); 970 PetscCall(PetscMalloc1(tab->s, &mprk->work_mediumbuffer)); 971 } 972 PetscCall(PetscMalloc1(tab->s, &mprk->work_fast)); 973 974 if (ts->use_splitrhsfunction) { 975 if (mprk->is_slow) { 976 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow)); 977 PetscCall(VecDuplicateVecs(YdotRHS_slow, tab->s, &mprk->YdotRHS_slow)); 978 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow)); 979 } 980 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer)); 981 PetscCall(VecDuplicateVecs(YdotRHS_slowbuffer, tab->s, &mprk->YdotRHS_slowbuffer)); 982 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer)); 983 if (tab->np == 3) { 984 if (mprk->is_medium) { 985 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium)); 986 PetscCall(VecDuplicateVecs(YdotRHS_medium, tab->s, &mprk->YdotRHS_medium)); 987 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium)); 988 } 989 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer)); 990 PetscCall(VecDuplicateVecs(YdotRHS_mediumbuffer, tab->s, &mprk->YdotRHS_mediumbuffer)); 991 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer)); 992 } 993 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast)); 994 PetscCall(VecDuplicateVecs(YdotRHS_fast, tab->s, &mprk->YdotRHS_fast)); 995 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast)); 996 } else { 997 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->YdotRHS)); 998 if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slow)); 999 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slowbuffer)); 1000 if (tab->np == 3) { 1001 if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_medium)); 1002 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_mediumbuffer)); 1003 } 1004 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_fast)); 1005 } 1006 PetscFunctionReturn(PETSC_SUCCESS); 1007 } 1008 1009 static PetscErrorCode TSSetUp_MPRK(TS ts) 1010 { 1011 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1012 MPRKTableau tab = mprk->tableau; 1013 DM dm; 1014 1015 PetscFunctionBegin; 1016 PetscCall(TSRHSSplitGetIS(ts, "slow", &mprk->is_slow)); 1017 PetscCall(TSRHSSplitGetIS(ts, "fast", &mprk->is_fast)); 1018 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); 1019 1020 if (tab->np == 3) { 1021 PetscCall(TSRHSSplitGetIS(ts, "medium", &mprk->is_medium)); 1022 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); 1023 PetscCall(TSRHSSplitGetIS(ts, "mediumbuffer", &mprk->is_mediumbuffer)); 1024 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */ 1025 mprk->is_mediumbuffer = mprk->is_medium; 1026 mprk->is_medium = NULL; 1027 } 1028 } 1029 1030 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */ 1031 PetscCall(TSRHSSplitGetIS(ts, "slowbuffer", &mprk->is_slowbuffer)); 1032 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */ 1033 mprk->is_slowbuffer = mprk->is_slow; 1034 mprk->is_slow = NULL; 1035 } 1036 PetscCall(TSCheckImplicitTerm(ts)); 1037 PetscCall(TSMPRKTableauSetUp(ts)); 1038 PetscCall(TSGetDM(ts, &dm)); 1039 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts)); 1040 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts)); 1041 if (ts->use_splitrhsfunction) { 1042 ts->ops->step = TSStep_MPRKSPLIT; 1043 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 1044 PetscCall(TSMPRKSetSplits(ts)); 1045 } else { 1046 ts->ops->step = TSStep_MPRK; 1047 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1048 } 1049 PetscFunctionReturn(PETSC_SUCCESS); 1050 } 1051 1052 static PetscErrorCode TSSetFromOptions_MPRK(TS ts, PetscOptionItems *PetscOptionsObject) 1053 { 1054 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1055 1056 PetscFunctionBegin; 1057 PetscOptionsHeadBegin(PetscOptionsObject, "PRK ODE solver options"); 1058 { 1059 MPRKTableauLink link; 1060 PetscInt count, choice; 1061 PetscBool flg; 1062 const char **namelist; 1063 for (link = MPRKTableauList, count = 0; link; link = link->next, count++) 1064 ; 1065 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1066 for (link = MPRKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1067 PetscCall(PetscOptionsEList("-ts_mprk_type", "Family of MPRK method", "TSMPRKSetType", (const char *const *)namelist, count, mprk->tableau->name, &choice, &flg)); 1068 if (flg) PetscCall(TSMPRKSetType(ts, namelist[choice])); 1069 PetscCall(PetscFree(namelist)); 1070 } 1071 PetscOptionsHeadEnd(); 1072 PetscFunctionReturn(PETSC_SUCCESS); 1073 } 1074 1075 static PetscErrorCode TSView_MPRK(TS ts, PetscViewer viewer) 1076 { 1077 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1078 PetscBool iascii; 1079 1080 PetscFunctionBegin; 1081 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1082 if (iascii) { 1083 MPRKTableau tab = mprk->tableau; 1084 TSMPRKType mprktype; 1085 char fbuf[512]; 1086 char sbuf[512]; 1087 PetscInt i; 1088 PetscCall(TSMPRKGetType(ts, &mprktype)); 1089 PetscCall(PetscViewerASCIIPrintf(viewer, " MPRK type %s\n", mprktype)); 1090 PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order)); 1091 1092 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->cf)); 1093 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cf = %s\n", fbuf)); 1094 PetscCall(PetscViewerASCIIPrintf(viewer, " Af = \n")); 1095 for (i = 0; i < tab->s; i++) { 1096 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, &tab->Af[i * tab->s])); 1097 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", fbuf)); 1098 } 1099 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->bf)); 1100 PetscCall(PetscViewerASCIIPrintf(viewer, " bf = %s\n", fbuf)); 1101 1102 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->csb)); 1103 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa csb = %s\n", sbuf)); 1104 PetscCall(PetscViewerASCIIPrintf(viewer, " Asb = \n")); 1105 for (i = 0; i < tab->s; i++) { 1106 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, &tab->Asb[i * tab->s])); 1107 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", sbuf)); 1108 } 1109 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->bsb)); 1110 PetscCall(PetscViewerASCIIPrintf(viewer, " bsb = %s\n", sbuf)); 1111 1112 if (tab->np == 3) { 1113 char mbuf[512]; 1114 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->cmb)); 1115 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cmb = %s\n", mbuf)); 1116 PetscCall(PetscViewerASCIIPrintf(viewer, " Amb = \n")); 1117 for (i = 0; i < tab->s; i++) { 1118 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, &tab->Amb[i * tab->s])); 1119 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", mbuf)); 1120 } 1121 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->bmb)); 1122 PetscCall(PetscViewerASCIIPrintf(viewer, " bmb = %s\n", mbuf)); 1123 } 1124 } 1125 PetscFunctionReturn(PETSC_SUCCESS); 1126 } 1127 1128 static PetscErrorCode TSLoad_MPRK(TS ts, PetscViewer viewer) 1129 { 1130 TSAdapt adapt; 1131 1132 PetscFunctionBegin; 1133 PetscCall(TSGetAdapt(ts, &adapt)); 1134 PetscCall(TSAdaptLoad(adapt, viewer)); 1135 PetscFunctionReturn(PETSC_SUCCESS); 1136 } 1137 1138 /*@C 1139 TSMPRKSetType - Set the type of `TSMPRK` scheme 1140 1141 Not Collective 1142 1143 Input Parameters: 1144 + ts - timestepping context 1145 - mprktype - type of `TSMPRK` scheme 1146 1147 Options Database Key: 1148 . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme 1149 1150 Level: intermediate 1151 1152 .seealso: [](chapter_ts), `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType` 1153 @*/ 1154 PetscErrorCode TSMPRKSetType(TS ts, TSMPRKType mprktype) 1155 { 1156 PetscFunctionBegin; 1157 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1158 PetscValidCharPointer(mprktype, 2); 1159 PetscTryMethod(ts, "TSMPRKSetType_C", (TS, TSMPRKType), (ts, mprktype)); 1160 PetscFunctionReturn(PETSC_SUCCESS); 1161 } 1162 1163 /*@C 1164 TSMPRKGetType - Get the type of `TSMPRK` scheme 1165 1166 Not Collective 1167 1168 Input Parameter: 1169 . ts - timestepping context 1170 1171 Output Parameter: 1172 . mprktype - type of `TSMPRK` scheme 1173 1174 Level: intermediate 1175 1176 .seealso: [](chapter_ts), `TSMPRKGetType()`, `TSMPRK` 1177 @*/ 1178 PetscErrorCode TSMPRKGetType(TS ts, TSMPRKType *mprktype) 1179 { 1180 PetscFunctionBegin; 1181 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1182 PetscUseMethod(ts, "TSMPRKGetType_C", (TS, TSMPRKType *), (ts, mprktype)); 1183 PetscFunctionReturn(PETSC_SUCCESS); 1184 } 1185 1186 static PetscErrorCode TSMPRKGetType_MPRK(TS ts, TSMPRKType *mprktype) 1187 { 1188 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1189 1190 PetscFunctionBegin; 1191 *mprktype = mprk->tableau->name; 1192 PetscFunctionReturn(PETSC_SUCCESS); 1193 } 1194 1195 static PetscErrorCode TSMPRKSetType_MPRK(TS ts, TSMPRKType mprktype) 1196 { 1197 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1198 PetscBool match; 1199 MPRKTableauLink link; 1200 1201 PetscFunctionBegin; 1202 if (mprk->tableau) { 1203 PetscCall(PetscStrcmp(mprk->tableau->name, mprktype, &match)); 1204 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1205 } 1206 for (link = MPRKTableauList; link; link = link->next) { 1207 PetscCall(PetscStrcmp(link->tab.name, mprktype, &match)); 1208 if (match) { 1209 if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts)); 1210 mprk->tableau = &link->tab; 1211 if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts)); 1212 PetscFunctionReturn(PETSC_SUCCESS); 1213 } 1214 } 1215 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mprktype); 1216 } 1217 1218 static PetscErrorCode TSGetStages_MPRK(TS ts, PetscInt *ns, Vec **Y) 1219 { 1220 TS_MPRK *mprk = (TS_MPRK *)ts->data; 1221 1222 PetscFunctionBegin; 1223 *ns = mprk->tableau->s; 1224 if (Y) *Y = mprk->Y; 1225 PetscFunctionReturn(PETSC_SUCCESS); 1226 } 1227 1228 static PetscErrorCode TSDestroy_MPRK(TS ts) 1229 { 1230 PetscFunctionBegin; 1231 PetscCall(TSReset_MPRK(ts)); 1232 if (ts->dm) { 1233 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts)); 1234 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts)); 1235 } 1236 PetscCall(PetscFree(ts->data)); 1237 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", NULL)); 1238 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", NULL)); 1239 PetscFunctionReturn(PETSC_SUCCESS); 1240 } 1241 1242 /*MC 1243 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes 1244 1245 The user should provide the right hand side of the equation using `TSSetRHSFunction()`. 1246 1247 Level: beginner 1248 1249 Note: 1250 The default is `TSMPRKPM2`, it can be changed with `TSMPRKSetType()` or -ts_mprk_type 1251 1252 .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()` 1253 `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2`, `TSType` 1254 M*/ 1255 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 1256 { 1257 TS_MPRK *mprk; 1258 1259 PetscFunctionBegin; 1260 PetscCall(TSMPRKInitializePackage()); 1261 1262 ts->ops->reset = TSReset_MPRK; 1263 ts->ops->destroy = TSDestroy_MPRK; 1264 ts->ops->view = TSView_MPRK; 1265 ts->ops->load = TSLoad_MPRK; 1266 ts->ops->setup = TSSetUp_MPRK; 1267 ts->ops->step = TSStep_MPRK; 1268 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 1269 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 1270 ts->ops->getstages = TSGetStages_MPRK; 1271 1272 PetscCall(PetscNew(&mprk)); 1273 ts->data = (void *)mprk; 1274 1275 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", TSMPRKGetType_MPRK)); 1276 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", TSMPRKSetType_MPRK)); 1277 1278 PetscCall(TSMPRKSetType(ts, TSMPRKDefault)); 1279 PetscFunctionReturn(PETSC_SUCCESS); 1280 } 1281