1 /* 2 Code for timestepping with Runge-Kutta method 3 4 Notes: 5 The general system is written as 6 7 F(t,U,Udot) = G(t,U) 8 9 where F represents the stiff part of the physics and G represents the non-stiff part. 10 11 */ 12 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 13 #include <petscdm.h> 14 15 static TSRKType TSRKDefault = TSRK3; 16 static PetscBool TSRKRegisterAllCalled; 17 static PetscBool TSRKPackageInitialized; 18 static PetscInt explicit_stage_time_id; 19 20 typedef struct _RKTableau *RKTableau; 21 struct _RKTableau { 22 char *name; 23 PetscInt order; /* Classical approximation order of the method i */ 24 PetscInt s; /* Number of stages */ 25 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 26 PetscInt pinterp; /* Interpolation order */ 27 PetscReal *A,*b,*c; /* Tableau */ 28 PetscReal *bembed; /* Embedded formula of order one less (order-1) */ 29 PetscReal *binterp; /* Dense output formula */ 30 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 31 }; 32 typedef struct _RKTableauLink *RKTableauLink; 33 struct _RKTableauLink { 34 struct _RKTableau tab; 35 RKTableauLink next; 36 }; 37 static RKTableauLink RKTableauList; 38 39 typedef struct { 40 RKTableau tableau; 41 Vec *Y; /* States computed during the step */ 42 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 43 PetscScalar *work; /* Scalar work */ 44 PetscReal stage_time; 45 TSStepStatus status; 46 } TS_RK; 47 48 /*MC 49 TSRK1 - First order forward Euler scheme. 50 51 This method has one stage. 52 53 Level: advanced 54 55 .seealso: TSRK 56 M*/ 57 /*MC 58 TSRK2 - Second order RK scheme. 59 60 This method has two stages. 61 62 Level: advanced 63 64 .seealso: TSRK 65 M*/ 66 /*MC 67 TSRK3 - Third order RK scheme. 68 69 This method has three stages. 70 71 Level: advanced 72 73 .seealso: TSRK 74 M*/ 75 /*MC 76 TSRK4 - Fourth order RK scheme. 77 78 This method has four stages. 79 80 Level: advanced 81 82 .seealso: TSRK 83 M*/ 84 /*MC 85 TSRK5F - Fifth order Fehlberg RK scheme with 4th order embedded method. 86 87 This method has six stages. 88 89 Level: advanced 90 91 .seealso: TSRK 92 M*/ 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "TSRKRegisterAll" 96 /*@C 97 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 98 99 Not Collective, but should be called by all processes which will need the schemes to be registered 100 101 Level: advanced 102 103 .keywords: TS, TSRK, register, all 104 105 .seealso: TSRKRegisterDestroy() 106 @*/ 107 PetscErrorCode TSRKRegisterAll(void) 108 { 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 113 TSRKRegisterAllCalled = PETSC_TRUE; 114 115 { 116 const PetscReal 117 A[1][1] = {{0.0}}, 118 b[1] = {1.0}; 119 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);CHKERRQ(ierr); 120 } 121 { 122 const PetscReal 123 A[2][2] = {{0.0,0.0}, 124 {1.0,0.0}}, 125 b[2] = {0.5,0.5}, 126 bembed[2] = {1.0,0}; 127 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);CHKERRQ(ierr); 128 } 129 { 130 const PetscReal 131 A[3][3] = {{0,0,0}, 132 {2.0/3.0,0,0}, 133 {-1.0/3.0,1.0,0}}, 134 b[3] = {0.25,0.5,0.25}; 135 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);CHKERRQ(ierr); 136 } 137 { 138 const PetscReal 139 A[4][4] = {{0,0,0,0}, 140 {1.0/2.0,0,0,0}, 141 {0,3.0/4.0,0,0}, 142 {2.0/9.0,1.0/3.0,4.0/9.0,0}}, 143 b[4] = {2.0/9.0,1.0/3.0,4.0/9.0,0}, 144 bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0}; 145 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);CHKERRQ(ierr); 146 } 147 { 148 const PetscReal 149 A[4][4] = {{0,0,0,0}, 150 {0.5,0,0,0}, 151 {0,0.5,0,0}, 152 {0,0,1.0,0}}, 153 b[4] = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0}; 154 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);CHKERRQ(ierr); 155 } 156 { 157 const PetscReal 158 A[6][6] = {{0,0,0,0,0,0}, 159 {0.25,0,0,0,0,0}, 160 {3.0/32.0,9.0/32.0,0,0,0,0}, 161 {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0}, 162 {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0}, 163 {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}}, 164 b[6] = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0}, 165 bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0}; 166 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 167 } 168 { 169 const PetscReal 170 A[7][7] = {{0,0,0,0,0,0,0}, 171 {1.0/5.0,0,0,0,0,0,0}, 172 {3.0/40.0,9.0/40.0,0,0,0,0,0}, 173 {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0}, 174 {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0}, 175 {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0}, 176 {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}}, 177 b[7] = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}, 178 bembed[7] = {5179.0/57600.0,0,7571.0/16695.0,393.0/640.0,-92097.0/339200.0,187.0/2100.0,1.0/40.0}; 179 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);CHKERRQ(ierr); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "TSRKRegisterDestroy" 186 /*@C 187 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 188 189 Not Collective 190 191 Level: advanced 192 193 .keywords: TSRK, register, destroy 194 .seealso: TSRKRegister(), TSRKRegisterAll() 195 @*/ 196 PetscErrorCode TSRKRegisterDestroy(void) 197 { 198 PetscErrorCode ierr; 199 RKTableauLink link; 200 201 PetscFunctionBegin; 202 while ((link = RKTableauList)) { 203 RKTableau t = &link->tab; 204 RKTableauList = link->next; 205 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 206 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 207 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 208 ierr = PetscFree (t->name); CHKERRQ(ierr); 209 ierr = PetscFree (link); CHKERRQ(ierr); 210 } 211 TSRKRegisterAllCalled = PETSC_FALSE; 212 PetscFunctionReturn(0); 213 } 214 215 #undef __FUNCT__ 216 #define __FUNCT__ "TSRKInitializePackage" 217 /*@C 218 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 219 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK() 220 when using static libraries. 221 222 Level: developer 223 224 .keywords: TS, TSRK, initialize, package 225 .seealso: PetscInitialize() 226 @*/ 227 PetscErrorCode TSRKInitializePackage(void) 228 { 229 PetscErrorCode ierr; 230 231 PetscFunctionBegin; 232 if (TSRKPackageInitialized) PetscFunctionReturn(0); 233 TSRKPackageInitialized = PETSC_TRUE; 234 ierr = TSRKRegisterAll();CHKERRQ(ierr); 235 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 236 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "TSRKFinalizePackage" 242 /*@C 243 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 244 called from PetscFinalize(). 245 246 Level: developer 247 248 .keywords: Petsc, destroy, package 249 .seealso: PetscFinalize() 250 @*/ 251 PetscErrorCode TSRKFinalizePackage(void) 252 { 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 TSRKPackageInitialized = PETSC_FALSE; 257 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "TSRKRegister" 263 /*@C 264 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 265 266 Not Collective, but the same schemes should be registered on all processes on which they will be used 267 268 Input Parameters: 269 + name - identifier for method 270 . order - approximation order of method 271 . s - number of stages, this is the dimension of the matrices below 272 . A - stage coefficients (dimension s*s, row-major) 273 . b - step completion table (dimension s; NULL to use last row of A) 274 . c - abscissa (dimension s; NULL to use row sums of A) 275 . bembed - completion table for embedded method (dimension s; NULL if not available) 276 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterp 277 - binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt) 278 279 Notes: 280 Several RK methods are provided, this function is only needed to create new methods. 281 282 Level: advanced 283 284 .keywords: TS, register 285 286 .seealso: TSRK 287 @*/ 288 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 289 const PetscReal A[],const PetscReal b[],const PetscReal c[], 290 const PetscReal bembed[], 291 PetscInt pinterp,const PetscReal binterp[]) 292 { 293 PetscErrorCode ierr; 294 RKTableauLink link; 295 RKTableau t; 296 PetscInt i,j; 297 298 PetscFunctionBegin; 299 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 300 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 301 t = &link->tab; 302 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 303 t->order = order; 304 t->s = s; 305 ierr = PetscMalloc3(s*s,PetscReal,&t->A,s,PetscReal,&t->b,s,PetscReal,&t->c);CHKERRQ(ierr); 306 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 307 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 308 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 309 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 310 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 311 t->FSAL = PETSC_TRUE; 312 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 313 if (bembed) { 314 ierr = PetscMalloc(s*sizeof(PetscReal),&t->bembed);CHKERRQ(ierr); 315 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 316 } 317 318 t->pinterp = pinterp; 319 ierr = PetscMalloc(s*pinterp*sizeof(PetscReal),&t->binterp);CHKERRQ(ierr); 320 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 321 link->next = RKTableauList; 322 RKTableauList = link; 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "TSEvaluateStep_RK" 328 /* 329 The step completion formula is 330 331 x1 = x0 + h b^T YdotRHS 332 333 This function can be called before or after ts->vec_sol has been updated. 334 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 335 We can write 336 337 x1e = x0 + h be^T YdotRHS 338 = x1 - h b^T YdotRHS + h be^T YdotRHS 339 = x1 + h (be - b)^T YdotRHS 340 341 so we can evaluate the method with different order even after the step has been optimistically completed. 342 */ 343 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 344 { 345 TS_RK *rk = (TS_RK*)ts->data; 346 RKTableau tab = rk->tableau; 347 PetscScalar *w = rk->work; 348 PetscReal h; 349 PetscInt s = tab->s,j; 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 switch (rk->status) { 354 case TS_STEP_INCOMPLETE: 355 case TS_STEP_PENDING: 356 h = ts->time_step; break; 357 case TS_STEP_COMPLETE: 358 h = ts->time_step_prev; break; 359 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 360 } 361 if (order == tab->order) { 362 if (rk->status == TS_STEP_INCOMPLETE) { 363 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 364 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 365 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 366 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 367 if (done) *done = PETSC_TRUE; 368 PetscFunctionReturn(0); 369 } else if (order == tab->order-1) { 370 if (!tab->bembed) goto unavailable; 371 if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */ 372 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 373 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 374 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 375 } else { /* Rollback and re-complete using (be-b) */ 376 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 377 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 378 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 379 } 380 if (done) *done = PETSC_TRUE; 381 PetscFunctionReturn(0); 382 } 383 unavailable: 384 if (done) *done = PETSC_FALSE; 385 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 386 PetscFunctionReturn(0); 387 } 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "TSStep_RK" 391 static PetscErrorCode TSStep_RK(TS ts) 392 { 393 TS_RK *rk = (TS_RK*)ts->data; 394 RKTableau tab = rk->tableau; 395 const PetscInt s = tab->s; 396 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 397 PetscScalar *w = rk->work; 398 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 399 TSAdapt adapt; 400 PetscInt i,j,reject,next_scheme; 401 PetscReal next_time_step; 402 PetscReal t; 403 PetscBool accept; 404 PetscErrorCode ierr; 405 406 PetscFunctionBegin; 407 408 next_time_step = ts->time_step; 409 t = ts->ptime; 410 accept = PETSC_TRUE; 411 rk->status = TS_STEP_INCOMPLETE; 412 413 414 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 415 PetscReal h = ts->time_step; 416 ierr = TSPreStep(ts);CHKERRQ(ierr); 417 for (i=0; i<s; i++) { 418 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 419 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 420 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 421 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 422 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 423 if (!accept) goto reject_step; 424 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 425 } 426 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 427 rk->status = TS_STEP_PENDING; 428 429 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 430 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 431 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 432 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 433 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 434 if (accept) { 435 /* ignore next_scheme for now */ 436 ts->ptime += ts->time_step; 437 ts->time_step = next_time_step; 438 ts->steps++; 439 rk->status = TS_STEP_COMPLETE; 440 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 441 442 break; 443 } else { /* Roll back the current step */ 444 for (j=0; j<s; j++) w[j] = -h*b[j]; 445 ierr = VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);CHKERRQ(ierr); 446 ts->time_step = next_time_step; 447 rk->status = TS_STEP_INCOMPLETE; 448 } 449 reject_step: continue; 450 } 451 if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 452 PetscFunctionReturn(0); 453 } 454 455 #undef __FUNCT__ 456 #define __FUNCT__ "TSInterpolate_RK" 457 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 458 { 459 TS_RK *rk = (TS_RK*)ts->data; 460 PetscInt s = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j; 461 PetscReal h; 462 PetscReal tt,t; 463 PetscScalar *b; 464 const PetscReal *B = rk->tableau->binterp; 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 469 switch (rk->status) { 470 case TS_STEP_INCOMPLETE: 471 case TS_STEP_PENDING: 472 h = ts->time_step; 473 t = (itime - ts->ptime)/h; 474 break; 475 case TS_STEP_COMPLETE: 476 h = ts->time_step_prev; 477 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 478 break; 479 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 480 } 481 ierr = PetscMalloc(s*sizeof(PetscScalar),&b);CHKERRQ(ierr); 482 for (i=0; i<s; i++) b[i] = 0; 483 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 484 for (i=0; i<s; i++) { 485 b[i] += h * B[i*pinterp+j] * tt; 486 } 487 } 488 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 489 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 490 ierr = PetscFree(b);CHKERRQ(ierr); 491 PetscFunctionReturn(0); 492 } 493 494 /*------------------------------------------------------------*/ 495 #undef __FUNCT__ 496 #define __FUNCT__ "TSReset_RK" 497 static PetscErrorCode TSReset_RK(TS ts) 498 { 499 TS_RK *rk = (TS_RK*)ts->data; 500 PetscInt s; 501 PetscErrorCode ierr; 502 503 PetscFunctionBegin; 504 if (!rk->tableau) PetscFunctionReturn(0); 505 s = rk->tableau->s; 506 ierr = VecDestroyVecs(s,&rk->Y);CHKERRQ(ierr); 507 ierr = VecDestroyVecs(s,&rk->YdotRHS);CHKERRQ(ierr); 508 ierr = PetscFree(rk->work);CHKERRQ(ierr); 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "TSDestroy_RK" 514 static PetscErrorCode TSDestroy_RK(TS ts) 515 { 516 PetscErrorCode ierr; 517 518 PetscFunctionBegin; 519 ierr = TSReset_RK(ts);CHKERRQ(ierr); 520 ierr = PetscFree(ts->data);CHKERRQ(ierr); 521 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 522 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 523 PetscFunctionReturn(0); 524 } 525 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "DMCoarsenHook_TSRK" 529 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 530 { 531 PetscFunctionBegin; 532 PetscFunctionReturn(0); 533 } 534 535 #undef __FUNCT__ 536 #define __FUNCT__ "DMRestrictHook_TSRK" 537 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 538 { 539 PetscFunctionBegin; 540 PetscFunctionReturn(0); 541 } 542 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "DMSubDomainHook_TSRK" 546 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 547 { 548 PetscFunctionBegin; 549 PetscFunctionReturn(0); 550 } 551 552 #undef __FUNCT__ 553 #define __FUNCT__ "DMSubDomainRestrictHook_TSRK" 554 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 555 { 556 557 PetscFunctionBegin; 558 PetscFunctionReturn(0); 559 } 560 561 #undef __FUNCT__ 562 #define __FUNCT__ "TSSetUp_RK" 563 static PetscErrorCode TSSetUp_RK(TS ts) 564 { 565 TS_RK *rk = (TS_RK*)ts->data; 566 RKTableau tab; 567 PetscInt s; 568 PetscErrorCode ierr; 569 DM dm; 570 571 PetscFunctionBegin; 572 if (!rk->tableau) { 573 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 574 } 575 tab = rk->tableau; 576 s = tab->s; 577 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->Y);CHKERRQ(ierr); 578 ierr = VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);CHKERRQ(ierr); 579 ierr = PetscMalloc(s*sizeof(rk->work[0]),&rk->work);CHKERRQ(ierr); 580 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 581 if (dm) { 582 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 583 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 584 } 585 PetscFunctionReturn(0); 586 } 587 /*------------------------------------------------------------*/ 588 589 #undef __FUNCT__ 590 #define __FUNCT__ "TSSetFromOptions_RK" 591 static PetscErrorCode TSSetFromOptions_RK(TS ts) 592 { 593 PetscErrorCode ierr; 594 char rktype[256]; 595 596 PetscFunctionBegin; 597 ierr = PetscOptionsHead("RK ODE solver options");CHKERRQ(ierr); 598 { 599 RKTableauLink link; 600 PetscInt count,choice; 601 PetscBool flg; 602 const char **namelist; 603 ierr = PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));CHKERRQ(ierr); 604 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 605 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 606 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 607 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);CHKERRQ(ierr); 608 ierr = TSRKSetType(ts,flg ? namelist[choice] : rktype);CHKERRQ(ierr); 609 ierr = PetscFree(namelist);CHKERRQ(ierr); 610 } 611 ierr = PetscOptionsTail();CHKERRQ(ierr); 612 PetscFunctionReturn(0); 613 } 614 615 #undef __FUNCT__ 616 #define __FUNCT__ "PetscFormatRealArray" 617 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 618 { 619 PetscErrorCode ierr; 620 PetscInt i; 621 size_t left,count; 622 char *p; 623 624 PetscFunctionBegin; 625 for (i=0,p=buf,left=len; i<n; i++) { 626 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 627 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 628 left -= count; 629 p += count; 630 *p++ = ' '; 631 } 632 p[i ? 0 : -1] = 0; 633 PetscFunctionReturn(0); 634 } 635 636 #undef __FUNCT__ 637 #define __FUNCT__ "TSView_RK" 638 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 639 { 640 TS_RK *rk = (TS_RK*)ts->data; 641 RKTableau tab = rk->tableau; 642 PetscBool iascii; 643 PetscErrorCode ierr; 644 TSAdapt adapt; 645 646 PetscFunctionBegin; 647 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 648 if (iascii) { 649 TSRKType rktype; 650 char buf[512]; 651 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 652 ierr = PetscViewerASCIIPrintf(viewer," RK %s\n",rktype);CHKERRQ(ierr); 653 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 654 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 655 ierr = PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 656 } 657 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 658 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "TSLoad_RK" 664 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 665 { 666 PetscErrorCode ierr; 667 TSAdapt tsadapt; 668 669 PetscFunctionBegin; 670 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 671 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "TSRKSetType" 677 /*@C 678 TSRKSetType - Set the type of RK scheme 679 680 Logically collective 681 682 Input Parameter: 683 + ts - timestepping context 684 - rktype - type of RK-scheme 685 686 Level: intermediate 687 688 .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5 689 @*/ 690 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 691 { 692 PetscErrorCode ierr; 693 694 PetscFunctionBegin; 695 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 696 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 700 #undef __FUNCT__ 701 #define __FUNCT__ "TSRKGetType" 702 /*@C 703 TSRKGetType - Get the type of RK scheme 704 705 Logically collective 706 707 Input Parameter: 708 . ts - timestepping context 709 710 Output Parameter: 711 . rktype - type of RK-scheme 712 713 Level: intermediate 714 715 .seealso: TSRKGetType() 716 @*/ 717 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 718 { 719 PetscErrorCode ierr; 720 721 PetscFunctionBegin; 722 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 723 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 724 PetscFunctionReturn(0); 725 } 726 727 #undef __FUNCT__ 728 #define __FUNCT__ "TSRKGetType_RK" 729 PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 730 { 731 TS_RK *rk = (TS_RK*)ts->data; 732 PetscErrorCode ierr; 733 734 PetscFunctionBegin; 735 if (!rk->tableau) { 736 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 737 } 738 *rktype = rk->tableau->name; 739 PetscFunctionReturn(0); 740 } 741 #undef __FUNCT__ 742 #define __FUNCT__ "TSRKSetType_RK" 743 PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 744 { 745 TS_RK *rk = (TS_RK*)ts->data; 746 PetscErrorCode ierr; 747 PetscBool match; 748 RKTableauLink link; 749 750 PetscFunctionBegin; 751 if (rk->tableau) { 752 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 753 if (match) PetscFunctionReturn(0); 754 } 755 for (link = RKTableauList; link; link=link->next) { 756 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 757 if (match) { 758 ierr = TSReset_RK(ts);CHKERRQ(ierr); 759 rk->tableau = &link->tab; 760 PetscFunctionReturn(0); 761 } 762 } 763 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 764 PetscFunctionReturn(0); 765 } 766 767 /* ------------------------------------------------------------ */ 768 /*MC 769 TSRK - ODE and DAE solver using Runge-Kutta schemes 770 771 The user should provide the right hand side of the equation 772 using TSSetRHSFunction(). 773 774 Notes: 775 The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type 776 777 Level: beginner 778 779 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 780 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister() 781 782 M*/ 783 #undef __FUNCT__ 784 #define __FUNCT__ "TSCreate_RK" 785 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 786 { 787 TS_RK *th; 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 792 ierr = TSRKInitializePackage();CHKERRQ(ierr); 793 #endif 794 795 ts->ops->reset = TSReset_RK; 796 ts->ops->destroy = TSDestroy_RK; 797 ts->ops->view = TSView_RK; 798 ts->ops->load = TSLoad_RK; 799 ts->ops->setup = TSSetUp_RK; 800 ts->ops->step = TSStep_RK; 801 ts->ops->interpolate = TSInterpolate_RK; 802 ts->ops->evaluatestep = TSEvaluateStep_RK; 803 ts->ops->setfromoptions = TSSetFromOptions_RK; 804 805 ierr = PetscNewLog(ts,TS_RK,&th);CHKERRQ(ierr); 806 ts->data = (void*)th; 807 808 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 809 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 810 PetscFunctionReturn(0); 811 } 812