1 /* 2 Code for time stepping with the Partitioned Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) for nonsplit RHS multi-rate RK, 7 user should give 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) for split RHS multi-rate RK, 11 user should partioned RHS by themselves and also provide the indexes for both slow and fast components. 12 */ 13 14 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 15 #include <petscdm.h> 16 17 static TSMPRKType TSMPRKDefault = TSMPRKPM2; 18 static PetscBool TSMPRKRegisterAllCalled; 19 static PetscBool TSMPRKPackageInitialized; 20 21 typedef struct _MPRKTableau *MPRKTableau; 22 struct _MPRKTableau { 23 char *name; 24 PetscInt order; /* Classical approximation order of the method i */ 25 PetscInt s; /* Number of stages */ 26 PetscReal *Af,*bf,*cf; /* Tableau for fast components */ 27 PetscReal *As,*bs,*cs; /* Tableau for slow components */ 28 }; 29 typedef struct _MPRKTableauLink *MPRKTableauLink; 30 struct _MPRKTableauLink { 31 struct _MPRKTableau tab; 32 MPRKTableauLink next; 33 }; 34 static MPRKTableauLink MPRKTableauList; 35 36 typedef struct { 37 MPRKTableau tableau; 38 TSMPRKMultirateType mprkmtype; 39 Vec *Y; /* States computed during the step */ 40 Vec Ytmp; 41 Vec *YdotRHS_fast; /* Function evaluations by fast tableau for fast components */ 42 Vec *YdotRHS_slow; /* Function evaluations by slow tableau for slow components */ 43 PetscScalar *work_fast; /* Scalar work_fast by fast tableau */ 44 PetscScalar *work_slow; /* Scalar work_slow by slow tableau */ 45 PetscReal stage_time; 46 TSStepStatus status; 47 PetscReal ptime; 48 PetscReal time_step; 49 IS is_slow,is_fast; 50 TS subts_slow,subts_fast; 51 } TS_MPRK; 52 53 /*MC 54 TSMPRKPM2 - Second Order Partitioned Runge Kutta scheme. 55 56 This method has four stages for both slow and fast parts. 57 58 Options database: 59 . -ts_mprk_type pm2 60 61 Level: advanced 62 63 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 64 M*/ 65 /*MC 66 TSMPRKPM3 - Third Order Partitioned Runge Kutta scheme. 67 68 This method has eight stages for both slow and fast parts. 69 70 Options database: 71 . -ts_mprk_type pm3 (put here temporarily) 72 73 Level: advanced 74 75 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 76 M*/ 77 /*MC 78 TSMPRKP2 - Second Order Partitioned Runge Kutta scheme. 79 80 This method has five stages for both slow and fast parts. 81 82 Options database: 83 . -ts_mprk_type p2 84 85 Level: advanced 86 87 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 88 M*/ 89 /*MC 90 TSMPRKP3 - Third Order Partitioned Runge Kutta scheme. 91 92 This method has ten stages for both slow and fast parts. 93 94 Options database: 95 . -ts_mprk_type p3 96 97 Level: advanced 98 99 .seealso: TSMPRK, TSMPRKType, TSMPRKSetType() 100 M*/ 101 102 /*@C 103 TSMPRKRegisterAll - Registers all of the Partirioned Runge-Kutta explicit methods in TSMPRK 104 105 Not Collective, but should be called by all processes which will need the schemes to be registered 106 107 Level: advanced 108 109 .keywords: TS, TSMPRK, register, all 110 111 .seealso: TSMPRKRegisterDestroy() 112 @*/ 113 PetscErrorCode TSMPRKRegisterAll(void) 114 { 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 if (TSMPRKRegisterAllCalled) PetscFunctionReturn(0); 119 TSMPRKRegisterAllCalled = PETSC_TRUE; 120 121 #define RC PetscRealConstant 122 { 123 const PetscReal 124 As[4][4] = {{0,0,0,0}, 125 {RC(1.0),0,0,0}, 126 {0,0,0,0}, 127 {0,0,RC(1.0),0}}, 128 A[4][4] = {{0,0,0,0}, 129 {RC(0.5),0,0,0}, 130 {RC(0.25),RC(0.25),0,0}, 131 {RC(0.25),RC(0.25),RC(0.5),0}}, 132 bs[4] = {RC(0.25),RC(0.25),RC(0.25),RC(0.25)}, 133 b[4] = {RC(0.25),RC(0.25),RC(0.25),RC(0.25)}; 134 ierr = TSMPRKRegister(TSMPRKPM2,2,4,&As[0][0],bs,NULL,&A[0][0],b,NULL);CHKERRQ(ierr); 135 } 136 137 /*{ 138 const PetscReal 139 As[8][8] = {{0,0,0,0,0,0,0,0}, 140 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0}, 141 {RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0,0,0,0,0}, 142 {RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0,0}, 143 {0,0,0,0,0,0,0,0}, 144 {0,0,0,0,RC(1.0)/RC(2.0),0,0,0}, 145 {0,0,0,0,RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0}, 146 {0,0,0,0,RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0}}, 147 A[8][8] = {{0,0,0,0,0,0,0,0}, 148 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0}, 149 {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0}, 150 {RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0,0}, 151 {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}, 152 {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}, 153 {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}, 154 {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}}, 155 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)}, 156 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)}; 157 ierr = TSMPRKRegister(TSMPRKPM3,3,8,&As[0][0],bs,NULL,&A[0][0],b,NULL);CHKERRQ(ierr); 158 }*/ 159 160 { 161 const PetscReal 162 As[5][5] = {{0,0,0,0,0}, 163 {RC(1.0)/RC(2.0),0,0,0,0}, 164 {RC(1.0)/RC(2.0),0,0,0,0}, 165 {RC(1.0),0,0,0,0}, 166 {RC(1.0),0,0,0,0}}, 167 A[5][5] = {{0,0,0,0,0}, 168 {RC(1.0)/RC(2.0),0,0,0,0}, 169 {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),0,0,0}, 170 {RC(1.0)/RC(4.0),RC(1.0)/RC(4.0),RC(1.0)/RC(2.0),0,0}, 171 {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}}, 172 bs[5] = {RC(1.0)/RC(2.0),0,0,0,RC(1.0)/RC(2.0)}, 173 b[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}; 174 ierr = TSMPRKRegister(TSMPRKP2,2,5,&As[0][0],bs,NULL,&A[0][0],b,NULL);CHKERRQ(ierr); 175 } 176 177 { 178 const PetscReal 179 As[10][10] = {{0,0,0,0,0,0,0,0,0,0}, 180 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 181 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 182 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0,0,0}, 183 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0,0,0}, 184 {RC(-1.0)/RC(6.0),0,0,0,RC(2.0)/RC(3.0),0,0,0,0,0}, 185 {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}, 186 {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}, 187 {RC(1.0)/RC(3.0),0,0,0,RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0}, 188 {RC(1.0)/RC(3.0),0,0,0,RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0}}, 189 A[10][10] = {{0,0,0,0,0,0,0,0,0,0}, 190 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0,0,0}, 191 {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0,0,0}, 192 {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}, 193 {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}, 194 {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}, 195 {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}, 196 {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}, 197 {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}, 198 {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}}, 199 bs[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)}, 200 b[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}; 201 ierr = TSMPRKRegister(TSMPRKP3,3,10,&As[0][0],bs,NULL,&A[0][0],b,NULL);CHKERRQ(ierr); 202 } 203 #undef RC 204 PetscFunctionReturn(0); 205 } 206 207 /*@C 208 TSMPRKRegisterDestroy - Frees the list of schemes that were registered by TSMPRKRegister(). 209 210 Not Collective 211 212 Level: advanced 213 214 .keywords: TSMPRK, register, destroy 215 .seealso: TSMPRKRegister(), TSMPRKRegisterAll() 216 @*/ 217 PetscErrorCode TSMPRKRegisterDestroy(void) 218 { 219 PetscErrorCode ierr; 220 MPRKTableauLink link; 221 222 PetscFunctionBegin; 223 while ((link = MPRKTableauList)) { 224 MPRKTableau t = &link->tab; 225 MPRKTableauList = link->next; 226 ierr = PetscFree3(t->Af,t->bf,t->cf);CHKERRQ(ierr); 227 ierr = PetscFree3(t->As,t->bs,t->cs);CHKERRQ(ierr); 228 ierr = PetscFree (t->name);CHKERRQ(ierr); 229 ierr = PetscFree (link);CHKERRQ(ierr); 230 } 231 TSMPRKRegisterAllCalled = PETSC_FALSE; 232 PetscFunctionReturn(0); 233 } 234 235 /*@C 236 TSMPRKInitializePackage - This function initializes everything in the TSMPRK package. It is called 237 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_MPRK() 238 when using static libraries. 239 240 Level: developer 241 242 .keywords: TS, TSMPRK, initialize, package 243 .seealso: PetscInitialize() 244 @*/ 245 PetscErrorCode TSMPRKInitializePackage(void) 246 { 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 if (TSMPRKPackageInitialized) PetscFunctionReturn(0); 251 TSMPRKPackageInitialized = PETSC_TRUE; 252 ierr = TSMPRKRegisterAll();CHKERRQ(ierr); 253 ierr = PetscRegisterFinalize(TSMPRKFinalizePackage);CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 /*@C 258 TSRKFinalizePackage - This function destroys everything in the TSMPRK package. It is 259 called from PetscFinalize(). 260 261 Level: developer 262 263 .keywords: Petsc, destroy, package 264 .seealso: PetscFinalize() 265 @*/ 266 PetscErrorCode TSMPRKFinalizePackage(void) 267 { 268 PetscErrorCode ierr; 269 270 PetscFunctionBegin; 271 TSMPRKPackageInitialized = PETSC_FALSE; 272 ierr = TSMPRKRegisterDestroy();CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 /*@C 277 TSMPRKRegister - register a MPRK scheme by providing the entries in the Butcher tableau 278 279 Not Collective, but the same schemes should be registered on all processes on which they will be used 280 281 Input Parameters: 282 + name - identifier for method 283 . order - approximation order of method 284 . s - number of stages, this is the dimension of the matrices below 285 . Af - stage coefficients for fast components(dimension s*s, row-major) 286 . bf - step completion table for fast components(dimension s) 287 . cf - abscissa for fast components(dimension s) 288 . As - stage coefficients for slow components(dimension s*s, row-major) 289 . bs - step completion table for slow components(dimension s) 290 - cs - abscissa for slow components(dimension s) 291 292 Notes: 293 Several MPRK methods are provided, this function is only needed to create new methods. 294 295 Level: advanced 296 297 .keywords: TS, register 298 299 .seealso: TSMPRK 300 @*/ 301 PetscErrorCode TSMPRKRegister(TSMPRKType name,PetscInt order,PetscInt s, 302 const PetscReal As[],const PetscReal bs[],const PetscReal cs[], 303 const PetscReal Af[],const PetscReal bf[],const PetscReal cf[]) 304 { 305 MPRKTableauLink link; 306 MPRKTableau t; 307 PetscInt i,j; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 PetscValidCharPointer(name,1); 312 PetscValidRealPointer(Af,7); 313 if (bf) PetscValidRealPointer(bf,8); 314 if (cf) PetscValidRealPointer(cf,9); 315 PetscValidRealPointer(As,4); 316 if (bs) PetscValidRealPointer(bs,5); 317 if (cs) PetscValidRealPointer(cs,6); 318 319 ierr = PetscNew(&link);CHKERRQ(ierr); 320 t = &link->tab; 321 322 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 323 t->order = order; 324 t->s = s; 325 ierr = PetscMalloc3(s*s,&t->Af,s,&t->bf,s,&t->cf);CHKERRQ(ierr); 326 ierr = PetscMemcpy(t->Af,Af,s*s*sizeof(Af[0]));CHKERRQ(ierr); 327 if (bf) { 328 ierr = PetscMemcpy(t->bf,bf,s*sizeof(bf[0]));CHKERRQ(ierr); 329 } 330 else 331 for (i=0; i<s; i++) t->bf[i] = Af[(s-1)*s+i]; 332 if (cf) { 333 ierr = PetscMemcpy(t->cf,cf,s*sizeof(cf[0]));CHKERRQ(ierr); 334 } 335 else { 336 for (i=0; i<s; i++) 337 for (j=0,t->cf[i]=0; j<s; j++) 338 t->cf[i] += Af[i*s+j]; 339 } 340 ierr = PetscMalloc3(s*s,&t->As,s,&t->bs,s,&t->cs);CHKERRQ(ierr); 341 ierr = PetscMemcpy(t->As,As,s*s*sizeof(As[0]));CHKERRQ(ierr); 342 if (bs) { 343 ierr = PetscMemcpy(t->bs,bs,s*sizeof(bs[0]));CHKERRQ(ierr); 344 } 345 else 346 for (i=0; i<s; i++) t->bs[i] = As[(s-1)*s+i]; 347 if (cs) { 348 ierr = PetscMemcpy(t->cs,cs,s*sizeof(cs[0]));CHKERRQ(ierr); 349 } 350 else { 351 for (i=0; i<s; i++) 352 for (j=0,t->cs[i]=0; j<s; j++) 353 t->cs[i] += As[i*s+j]; 354 } 355 link->next = MPRKTableauList; 356 MPRKTableauList = link; 357 PetscFunctionReturn(0); 358 } 359 360 static PetscErrorCode TSMPRKSetSplits(TS ts) 361 { 362 TS_MPRK *mprk = (TS_MPRK*)ts->data; 363 DM dm,subdm,newdm; 364 PetscErrorCode ierr; 365 366 PetscFunctionBegin; 367 ierr = TSRHSSplitGetSubTS(ts,"slow",&mprk->subts_slow);CHKERRQ(ierr); 368 ierr = TSRHSSplitGetSubTS(ts,"fast",&mprk->subts_fast);CHKERRQ(ierr); 369 if (!mprk->subts_slow || !mprk->subts_fast) SETERRQ(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"); 370 371 /* Only copy */ 372 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 373 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 374 ierr = TSGetDM(mprk->subts_fast,&subdm);CHKERRQ(ierr); 375 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 376 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 377 ierr = TSSetDM(mprk->subts_fast,newdm);CHKERRQ(ierr); 378 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 379 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 380 ierr = TSGetDM(mprk->subts_slow,&subdm);CHKERRQ(ierr); 381 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 382 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 383 ierr = TSSetDM(mprk->subts_slow,newdm);CHKERRQ(ierr); 384 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 385 PetscFunctionReturn(0); 386 } 387 388 /* 389 This if for nonsplit RHS MPRK 390 The step completion formula is 391 392 x1 = x0 + h b^T YdotRHS 393 394 */ 395 static PetscErrorCode TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool *done) 396 { 397 TS_MPRK *mprk = (TS_MPRK*)ts->data; 398 MPRKTableau tab = mprk->tableau; 399 Vec Xtmp = mprk->Ytmp,Xslow,Xfast; 400 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow; 401 PetscReal h = ts->time_step; 402 PetscInt s = tab->s,j; 403 PetscErrorCode ierr; 404 405 PetscFunctionBegin; 406 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 407 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 408 for (j=0; j<s; j++) ws[j] = h*tab->bs[j]; 409 ierr = VecCopy(X,Xtmp);CHKERRQ(ierr); 410 ierr = VecMAXPY(Xtmp,s,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 411 ierr = VecGetSubVector(Xtmp,mprk->is_slow,&Xslow);CHKERRQ(ierr); 412 ierr = VecISCopy(X,mprk->is_slow,SCATTER_FORWARD,Xslow);CHKERRQ(ierr); 413 ierr = VecRestoreSubVector(Xtmp,mprk->is_slow,&Xslow);CHKERRQ(ierr); 414 415 /* Update fast part of X, note that the slow part has been changed but is simply discarded here */ 416 ierr = VecCopy(X,Xtmp);CHKERRQ(ierr); 417 ierr = VecMAXPY(Xtmp,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 418 ierr = VecGetSubVector(Xtmp,mprk->is_fast,&Xfast);CHKERRQ(ierr); 419 ierr = VecISCopy(X,mprk->is_fast,SCATTER_FORWARD,Xfast);CHKERRQ(ierr); 420 ierr = VecRestoreSubVector(Xtmp,mprk->is_fast,&Xfast);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 static PetscErrorCode TSStep_MPRK(TS ts) 425 { 426 TS_MPRK *mprk = (TS_MPRK*)ts->data; 427 Vec *Y = mprk->Y,Ytmp = mprk->Ytmp,*YdotRHS_fast = mprk->YdotRHS_fast,*YdotRHS_slow = mprk->YdotRHS_slow; 428 Vec Yfast,Yslow; 429 MPRKTableau tab = mprk->tableau; 430 const PetscInt s = tab->s; 431 const PetscReal *Af = tab->Af,*cf = tab->cf,*As = tab->As,*cs = tab->cs; 432 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow; 433 PetscInt i,j; 434 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 for (i=0; i<s; i++) { 439 mprk->stage_time = t + h*cf[i]; 440 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 441 442 /* update the satge value for all components by slow and fast tableau respectively */ 443 for (j=0; j<i; j++) ws[j] = h*As[i*s+j]; 444 ierr = VecCopy(ts->vec_sol,Ytmp);CHKERRQ(ierr); 445 ierr = VecMAXPY(Ytmp,i,ws,YdotRHS_slow);CHKERRQ(ierr); 446 ierr = VecGetSubVector(Ytmp,mprk->is_slow,&Yslow);CHKERRQ(ierr); 447 ierr = VecISCopy(Y[i],mprk->is_slow,SCATTER_FORWARD,Yslow);CHKERRQ(ierr); 448 ierr = VecRestoreSubVector(Ytmp,mprk->is_slow,&Yslow);CHKERRQ(ierr); 449 450 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 451 ierr = VecCopy(ts->vec_sol,Ytmp);CHKERRQ(ierr); 452 ierr = VecMAXPY(Ytmp,i,wf,YdotRHS_fast);CHKERRQ(ierr); 453 ierr = VecGetSubVector(Ytmp,mprk->is_fast,&Yfast);CHKERRQ(ierr); 454 ierr = VecISCopy(Y[i],mprk->is_fast,SCATTER_FORWARD,Yfast);CHKERRQ(ierr); 455 ierr = VecRestoreSubVector(Ytmp,mprk->is_fast,&Yfast);CHKERRQ(ierr); 456 457 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 458 /* compute the stage RHS by fast and slow tableau respectively */ 459 ierr = TSComputeRHSFunction(ts,t+h*cs[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 460 ierr = TSComputeRHSFunction(ts,t+h*cf[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 461 } 462 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 463 ts->ptime += ts->time_step; 464 ts->time_step = next_time_step; 465 PetscFunctionReturn(0); 466 } 467 468 /* 469 This if for partitioned RHS MPRK 470 The step completion formula is 471 472 x1 = x0 + h b^T YdotRHS 473 474 */ 475 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 476 { 477 TS_MPRK *mprk = (TS_MPRK*)ts->data; 478 MPRKTableau tab = mprk->tableau; 479 Vec Xslow,Xfast; /* subvectors for slow and fast componets in X respectively */ 480 PetscScalar *wf = mprk->work_fast,*ws = mprk->work_slow; 481 PetscReal h = ts->time_step; 482 PetscInt s = tab->s,j; 483 PetscErrorCode ierr; 484 485 PetscFunctionBegin; 486 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 487 for (j=0; j<s; j++) wf[j] = h*tab->bf[j]; 488 for (j=0; j<s; j++) ws[j] = h*tab->bs[j]; 489 ierr = VecGetSubVector(X,mprk->is_slow,&Xslow);CHKERRQ(ierr); 490 ierr = VecGetSubVector(X,mprk->is_fast,&Xfast);CHKERRQ(ierr); 491 ierr = VecMAXPY(Xslow,s,ws,mprk->YdotRHS_slow);CHKERRQ(ierr); 492 ierr = VecMAXPY(Xfast,s,wf,mprk->YdotRHS_fast);CHKERRQ(ierr); 493 ierr = VecRestoreSubVector(X,mprk->is_slow,&Xfast);CHKERRQ(ierr); 494 ierr = VecRestoreSubVector(X,mprk->is_fast,&Xslow);CHKERRQ(ierr); 495 PetscFunctionReturn(0); 496 } 497 498 static PetscErrorCode TSStep_MPRKSPLIT(TS ts) 499 { 500 TS_MPRK *mprk = (TS_MPRK*)ts->data; 501 MPRKTableau tab = mprk->tableau; 502 Vec *Y = mprk->Y,*YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow; 503 Vec Yslow,Yfast; /* subvectors for slow and fast components in Y[i] respectively */ 504 const PetscInt s = tab->s; 505 const PetscReal *Af = tab->Af,*cf = tab->cf,*As = tab->As,*cs = tab->cs; 506 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow; 507 PetscInt i,j; 508 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 509 PetscErrorCode ierr; 510 511 PetscFunctionBegin; 512 for (i=0; i<s; i++) { 513 mprk->stage_time = t + h*cf[i]; 514 ierr = TSPreStage(ts,mprk->stage_time);CHKERRQ(ierr); 515 /* calculate the stage value for fast and slow components respectively */ 516 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 517 for (j=0; j<i; j++) wf[j] = h*Af[i*s+j]; 518 for (j=0; j<i; j++) ws[j] = h*As[i*s+j]; 519 ierr = VecGetSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 520 ierr = VecGetSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 521 ierr = VecMAXPY(Yslow,i,ws,YdotRHS_slow);CHKERRQ(ierr); 522 ierr = VecMAXPY(Yfast,i,wf,YdotRHS_fast);CHKERRQ(ierr); 523 ierr = VecRestoreSubVector(Y[i],mprk->is_slow,&Yslow);CHKERRQ(ierr); 524 ierr = VecRestoreSubVector(Y[i],mprk->is_fast,&Yfast);CHKERRQ(ierr); 525 ierr = TSPostStage(ts,mprk->stage_time,i,Y); CHKERRQ(ierr); 526 /* calculate the stage RHS for slow and fast components respectively */ 527 ierr = TSComputeRHSFunction(mprk->subts_slow,t+h*cs[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 528 ierr = TSComputeRHSFunction(mprk->subts_fast,t+h*cf[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 529 } 530 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 531 ts->ptime += ts->time_step; 532 ts->time_step = next_time_step; 533 PetscFunctionReturn(0); 534 } 535 536 static PetscErrorCode TSMPRKTableauReset(TS ts) 537 { 538 TS_MPRK *mprk = (TS_MPRK*)ts->data; 539 MPRKTableau tab = mprk->tableau; 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 if (!tab) PetscFunctionReturn(0); 544 ierr = PetscFree(mprk->work_fast);CHKERRQ(ierr); 545 ierr = PetscFree(mprk->work_slow);CHKERRQ(ierr); 546 ierr = VecDestroyVecs(tab->s,&mprk->Y);CHKERRQ(ierr); 547 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 548 ierr = VecDestroy(&mprk->Ytmp);CHKERRQ(ierr); 549 } 550 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 551 ierr = VecDestroyVecs(tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 static PetscErrorCode TSReset_MPRK(TS ts) 556 { 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 564 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine,DM coarse,void *ctx) 565 { 566 PetscFunctionBegin; 567 PetscFunctionReturn(0); 568 } 569 570 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 571 { 572 PetscFunctionBegin; 573 PetscFunctionReturn(0); 574 } 575 576 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm,DM subdm,void *ctx) 577 { 578 PetscFunctionBegin; 579 PetscFunctionReturn(0); 580 } 581 582 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 583 { 584 PetscFunctionBegin; 585 PetscFunctionReturn(0); 586 } 587 588 static PetscErrorCode TSMPRKTableauSetUp(TS ts) 589 { 590 TS_MPRK *mprk = (TS_MPRK*)ts->data; 591 MPRKTableau tab = mprk->tableau; 592 Vec YdotRHS_fast,YdotRHS_slow; 593 PetscErrorCode ierr; 594 595 PetscFunctionBegin; 596 ierr = PetscMalloc1(tab->s,&mprk->work_fast);CHKERRQ(ierr); 597 ierr = PetscMalloc1(tab->s,&mprk->work_slow);CHKERRQ(ierr); 598 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->Y);CHKERRQ(ierr); 599 if (mprk->mprkmtype == TSMPRKNONSPLIT) { 600 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 601 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 602 ierr = VecDuplicate(ts->vec_sol,&mprk->Ytmp);CHKERRQ(ierr); 603 } 604 if (mprk->mprkmtype == TSMPRKSPLIT) { 605 ierr = VecGetSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 606 ierr = VecGetSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 607 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&mprk->YdotRHS_slow);CHKERRQ(ierr); 608 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&mprk->YdotRHS_fast);CHKERRQ(ierr); 609 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_slow,&YdotRHS_slow);CHKERRQ(ierr); 610 ierr = VecRestoreSubVector(ts->vec_sol,mprk->is_fast,&YdotRHS_fast);CHKERRQ(ierr); 611 } 612 PetscFunctionReturn(0); 613 } 614 615 static PetscErrorCode TSSetUp_MPRK(TS ts) 616 { 617 TS_MPRK *mprk = (TS_MPRK*)ts->data; 618 DM dm; 619 PetscErrorCode ierr; 620 621 PetscFunctionBegin; 622 ierr = TSRHSSplitGetIS(ts,"slow",&mprk->is_slow);CHKERRQ(ierr); 623 ierr = TSRHSSplitGetIS(ts,"fast",&mprk->is_fast);CHKERRQ(ierr); 624 if (!mprk->is_slow || !mprk->is_fast) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'fast' respectively in order to use -ts_type mprk"); 625 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 626 ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr); 627 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 628 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 629 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 630 ierr = PetscTryMethod(ts,"TSMPRKSetSplits_C",(TS),(ts));CHKERRQ(ierr); 631 PetscFunctionReturn(0); 632 } 633 634 /* construct a database to chose nonsplit RHS mutirate mprk method or split RHS MPRK method */ 635 const char *const TSMPRKMultirateTypes[] = {"NONSPLIT","SPLIT","TSMPRKMultirateType","TSMPRK",0}; 636 637 static PetscErrorCode TSSetFromOptions_MPRK(PetscOptionItems *PetscOptionsObject,TS ts) 638 { 639 TS_MPRK *mprk = (TS_MPRK*)ts->data; 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 ierr = PetscOptionsHead(PetscOptionsObject,"PRK ODE solver options");CHKERRQ(ierr); 644 { 645 MPRKTableauLink link; 646 PetscInt count,choice; 647 PetscBool flg; 648 const char **namelist; 649 PetscInt mprkmtype = 0; 650 for (link=MPRKTableauList,count=0; link; link=link->next,count++) ; 651 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 652 for (link=MPRKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 653 ierr = PetscOptionsEList("-ts_mprk_type","Family of MPRK method","TSMPRKSetType",(const char*const*)namelist,count,mprk->tableau->name,&choice,&flg);CHKERRQ(ierr); 654 if (flg) {ierr = TSMPRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 655 ierr = PetscFree(namelist);CHKERRQ(ierr); 656 ierr = PetscOptionsEList("-ts_mprk_multirate_type","Use Combined RHS Multirate or Partioned RHS Multirat MPRK method","TSMPRKSetMultirateType",TSMPRKMultirateTypes,2,TSMPRKMultirateTypes[0],&mprkmtype,&flg);CHKERRQ(ierr); 657 if (flg) {ierr = TSMPRKSetMultirateType(ts,mprkmtype);CHKERRQ(ierr);} 658 } 659 ierr = PetscOptionsTail();CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 static PetscErrorCode TSView_MPRK(TS ts,PetscViewer viewer) 664 { 665 TS_MPRK *mprk = (TS_MPRK*)ts->data; 666 PetscBool iascii; 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 671 if (iascii) { 672 MPRKTableau tab = mprk->tableau; 673 TSMPRKType mprktype; 674 char fbuf[512]; 675 char sbuf[512]; 676 ierr = TSMPRKGetType(ts,&mprktype);CHKERRQ(ierr); 677 ierr = PetscViewerASCIIPrintf(viewer," MPRK type %s\n",mprktype);CHKERRQ(ierr); 678 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 679 ierr = PetscFormatRealArray(fbuf,sizeof(fbuf),"% 8.6f",tab->s,tab->cf);CHKERRQ(ierr); 680 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cf = %s\n",fbuf);CHKERRQ(ierr); 681 ierr = PetscFormatRealArray(sbuf,sizeof(sbuf),"% 8.6f",tab->s,tab->cs);CHKERRQ(ierr); 682 ierr = PetscViewerASCIIPrintf(viewer," Abscissa cs = %s\n",sbuf);CHKERRQ(ierr); 683 } 684 PetscFunctionReturn(0); 685 } 686 687 static PetscErrorCode TSLoad_MPRK(TS ts,PetscViewer viewer) 688 { 689 PetscErrorCode ierr; 690 TSAdapt adapt; 691 692 PetscFunctionBegin; 693 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 694 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 /*@C 699 TSMPRKSetType - Set the type of MPRK scheme 700 701 Logically collective 702 703 Input Parameter: 704 + ts - timestepping context 705 - mprktype - type of MPRK-scheme 706 707 Options Database: 708 . -ts_mprk_type - <pm2,p2,p3> 709 710 Level: intermediate 711 712 .seealso: TSMPRKGetType(), TSMPRK, TSMPRKType 713 @*/ 714 PetscErrorCode TSMPRKSetType(TS ts,TSMPRKType mprktype) 715 { 716 PetscErrorCode ierr; 717 718 PetscFunctionBegin; 719 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 720 PetscValidCharPointer(mprktype,2); 721 ierr = PetscTryMethod(ts,"TSMPRKSetType_C",(TS,TSMPRKType),(ts,mprktype));CHKERRQ(ierr); 722 PetscFunctionReturn(0); 723 } 724 725 /*@C 726 TSMPRKGetType - Get the type of MPRK scheme 727 728 Logically collective 729 730 Input Parameter: 731 . ts - timestepping context 732 733 Output Parameter: 734 . mprktype - type of MPRK-scheme 735 736 Level: intermediate 737 738 .seealso: TSMPRKGetType() 739 @*/ 740 PetscErrorCode TSMPRKGetType(TS ts,TSMPRKType *mprktype) 741 { 742 PetscErrorCode ierr; 743 744 PetscFunctionBegin; 745 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 746 ierr = PetscUseMethod(ts,"TSMPRKGetType_C",(TS,TSMPRKType*),(ts,mprktype));CHKERRQ(ierr); 747 PetscFunctionReturn(0); 748 } 749 750 /*@C 751 TSMPRKSetMultirateType - Set the type of MPRK multirate scheme 752 753 Logically collective 754 755 Input Parameter: 756 + ts - timestepping context 757 - mprkmtype - type of the multirate configuration 758 759 Options Database: 760 . -ts_mprk_multirate_type - <nonsplit,split> 761 762 Level: intermediate 763 @*/ 764 PetscErrorCode TSMPRKSetMultirateType(TS ts, TSMPRKMultirateType mprkmtype) 765 { 766 TS_MPRK *mprk = (TS_MPRK*)ts->data; 767 PetscErrorCode ierr; 768 769 PetscFunctionBegin; 770 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 771 switch(mprkmtype){ 772 case TSMPRKNONSPLIT: 773 ts->ops->step = TSStep_MPRK; 774 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 775 break; 776 case TSMPRKSPLIT: 777 ts->ops->step = TSStep_MPRKSPLIT; 778 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT; 779 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetSplits_C",TSMPRKSetSplits);CHKERRQ(ierr); 780 break; 781 default : 782 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",mprkmtype); 783 } 784 mprk->mprkmtype = mprkmtype; 785 PetscFunctionReturn(0); 786 } 787 788 static PetscErrorCode TSMPRKGetType_MPRK(TS ts,TSMPRKType *mprktype) 789 { 790 TS_MPRK *mprk = (TS_MPRK*)ts->data; 791 792 PetscFunctionBegin; 793 *mprktype = mprk->tableau->name; 794 PetscFunctionReturn(0); 795 } 796 797 static PetscErrorCode TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype) 798 { 799 TS_MPRK *mprk = (TS_MPRK*)ts->data; 800 PetscBool match; 801 MPRKTableauLink link; 802 PetscErrorCode ierr; 803 804 PetscFunctionBegin; 805 if (mprk->tableau) { 806 ierr = PetscStrcmp(mprk->tableau->name,mprktype,&match);CHKERRQ(ierr); 807 if (match) PetscFunctionReturn(0); 808 } 809 for (link = MPRKTableauList; link; link=link->next) { 810 ierr = PetscStrcmp(link->tab.name,mprktype,&match);CHKERRQ(ierr); 811 if (match) { 812 if (ts->setupcalled) {ierr = TSMPRKTableauReset(ts);CHKERRQ(ierr);} 813 mprk->tableau = &link->tab; 814 if (ts->setupcalled) {ierr = TSMPRKTableauSetUp(ts);CHKERRQ(ierr);} 815 PetscFunctionReturn(0); 816 } 817 } 818 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mprktype); 819 PetscFunctionReturn(0); 820 } 821 822 static PetscErrorCode TSGetStages_MPRK(TS ts,PetscInt *ns,Vec **Y) 823 { 824 TS_MPRK *mprk = (TS_MPRK*)ts->data; 825 826 PetscFunctionBegin; 827 *ns = mprk->tableau->s; 828 if (Y) *Y = mprk->Y; 829 PetscFunctionReturn(0); 830 } 831 832 static PetscErrorCode TSDestroy_MPRK(TS ts) 833 { 834 PetscErrorCode ierr; 835 836 PetscFunctionBegin; 837 ierr = TSReset_MPRK(ts);CHKERRQ(ierr); 838 if (ts->dm) { 839 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSMPRK,DMRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 840 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSMPRK,DMSubDomainRestrictHook_TSMPRK,ts);CHKERRQ(ierr); 841 } 842 ierr = PetscFree(ts->data);CHKERRQ(ierr); 843 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",NULL);CHKERRQ(ierr); 844 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",NULL);CHKERRQ(ierr); 845 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",NULL);CHKERRQ(ierr); 846 PetscFunctionReturn(0); 847 } 848 849 /*MC 850 TSMPRK - ODE solver using Partitioned Runge-Kutta schemes 851 852 The user should provide the right hand side of the equation 853 using TSSetRHSFunction(). 854 855 Notes: 856 The default is TSMPRKPM2, it can be changed with TSRKSetType() or -ts_mprk_type 857 858 Level: beginner 859 860 .seealso: TSCreate(), TS, TSSetType(), TSMPRKSetType(), TSMPRKGetType(), TSMPRKType, TSMPRKRegister(), TSMPRKSetMultirateType() 861 TSMPRKM2, TSMPRKM3, TSMPRKRFSMR3, TSMPRKRFSMR2 862 863 M*/ 864 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts) 865 { 866 TS_MPRK *mprk; 867 PetscErrorCode ierr; 868 869 PetscFunctionBegin; 870 ierr = TSMPRKInitializePackage();CHKERRQ(ierr); 871 872 ts->ops->reset = TSReset_MPRK; 873 ts->ops->destroy = TSDestroy_MPRK; 874 ts->ops->view = TSView_MPRK; 875 ts->ops->load = TSLoad_MPRK; 876 ts->ops->setup = TSSetUp_MPRK; 877 ts->ops->step = TSStep_MPRK; 878 ts->ops->evaluatestep = TSEvaluateStep_MPRK; 879 ts->ops->setfromoptions = TSSetFromOptions_MPRK; 880 ts->ops->getstages = TSGetStages_MPRK; 881 882 ierr = PetscNewLog(ts,&mprk);CHKERRQ(ierr); 883 ts->data = (void*)mprk; 884 885 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKGetType_C",TSMPRKGetType_MPRK);CHKERRQ(ierr); 886 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetType_C",TSMPRKSetType_MPRK);CHKERRQ(ierr); 887 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSMPRKSetMultirateType_C",TSMPRKSetMultirateType);CHKERRQ(ierr); 888 889 ierr = TSMPRKSetType(ts,TSMPRKDefault);CHKERRQ(ierr); 890 PetscFunctionReturn(0); 891 } 892