1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 4 PetscClassId TSADAPT_CLASSID; 5 6 static PetscFunctionList TSAdaptList; 7 static PetscBool TSAdaptPackageInitialized; 8 static PetscBool TSAdaptRegisterAllCalled; 9 10 PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt); 11 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt); 12 PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt); 13 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt); 14 15 /*@C 16 TSAdaptRegister - adds a TSAdapt implementation 17 18 Not Collective 19 20 Input Parameters: 21 + name_scheme - name of user-defined adaptivity scheme 22 - routine_create - routine to create method context 23 24 Notes: 25 TSAdaptRegister() may be called multiple times to add several user-defined families. 26 27 Sample usage: 28 .vb 29 TSAdaptRegister("my_scheme",MySchemeCreate); 30 .ve 31 32 Then, your scheme can be chosen with the procedural interface via 33 $ TSAdaptSetType(ts,"my_scheme") 34 or at runtime via the option 35 $ -ts_adapt_type my_scheme 36 37 Level: advanced 38 39 .keywords: TSAdapt, register 40 41 .seealso: TSAdaptRegisterAll() 42 @*/ 43 PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt)) 44 { 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = PetscFunctionListAdd(&TSAdaptList,sname,function);CHKERRQ(ierr); 49 PetscFunctionReturn(0); 50 } 51 52 /*@C 53 TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt 54 55 Not Collective 56 57 Level: advanced 58 59 .keywords: TSAdapt, register, all 60 61 .seealso: TSAdaptRegisterDestroy() 62 @*/ 63 PetscErrorCode TSAdaptRegisterAll(void) 64 { 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 if (TSAdaptRegisterAllCalled) PetscFunctionReturn(0); 69 TSAdaptRegisterAllCalled = PETSC_TRUE; 70 ierr = TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);CHKERRQ(ierr); 71 ierr = TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);CHKERRQ(ierr); 72 ierr = TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);CHKERRQ(ierr); 73 ierr = TSAdaptRegister(TSADAPTGLEE ,TSAdaptCreate_GLEE);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 /*@C 78 TSAdaptFinalizePackage - This function destroys everything in the TS package. It is 79 called from PetscFinalize(). 80 81 Level: developer 82 83 .keywords: Petsc, destroy, package 84 .seealso: PetscFinalize() 85 @*/ 86 PetscErrorCode TSAdaptFinalizePackage(void) 87 { 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = PetscFunctionListDestroy(&TSAdaptList);CHKERRQ(ierr); 92 TSAdaptPackageInitialized = PETSC_FALSE; 93 TSAdaptRegisterAllCalled = PETSC_FALSE; 94 PetscFunctionReturn(0); 95 } 96 97 /*@C 98 TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is 99 called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to 100 TSAdaptCreate() when using static libraries. 101 102 Level: developer 103 104 .keywords: TSAdapt, initialize, package 105 .seealso: PetscInitialize() 106 @*/ 107 PetscErrorCode TSAdaptInitializePackage(void) 108 { 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 if (TSAdaptPackageInitialized) PetscFunctionReturn(0); 113 TSAdaptPackageInitialized = PETSC_TRUE; 114 ierr = PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);CHKERRQ(ierr); 115 ierr = TSAdaptRegisterAll();CHKERRQ(ierr); 116 ierr = PetscRegisterFinalize(TSAdaptFinalizePackage);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 /*@C 121 TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE 122 123 Logicially Collective on TSAdapt 124 125 Input Parameter: 126 + adapt - the TS adapter, most likely obtained with TSGetAdapt() 127 - type - either TSADAPTBASIC or TSADAPTNONE 128 129 Options Database: 130 . -ts_adapt_type basic or none - to setting the adapter type 131 132 Level: intermediate 133 134 .keywords: TSAdapt, create 135 136 .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType() 137 @*/ 138 PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type) 139 { 140 PetscBool match; 141 PetscErrorCode ierr,(*r)(TSAdapt); 142 143 PetscFunctionBegin; 144 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 145 ierr = PetscObjectTypeCompare((PetscObject)adapt,type,&match);CHKERRQ(ierr); 146 if (match) PetscFunctionReturn(0); 147 ierr = PetscFunctionListFind(TSAdaptList,type,&r);CHKERRQ(ierr); 148 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type); 149 if (adapt->ops->destroy) {ierr = (*adapt->ops->destroy)(adapt);CHKERRQ(ierr);} 150 ierr = PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));CHKERRQ(ierr); 151 ierr = PetscObjectChangeTypeName((PetscObject)adapt,type);CHKERRQ(ierr); 152 ierr = (*r)(adapt);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 /*@C 157 TSAdaptGetType - gets the TS adapter method type (as a string). 158 159 Not Collective 160 161 Input Parameter: 162 . adapt - The TS adapter, most likely obtained with TSGetAdapt() 163 164 Output Parameter: 165 . type - The name of TS adapter method 166 167 Level: intermediate 168 169 .keywords: TSAdapt, get, type 170 .seealso TSAdaptSetType() 171 @*/ 172 PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type) 173 { 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 176 PetscValidPointer(type,2); 177 *type = ((PetscObject)adapt)->type_name; 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[]) 182 { 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 187 ierr = PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 /*@C 192 TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView(). 193 194 Collective on PetscViewer 195 196 Input Parameters: 197 + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or 198 some related function before a call to TSAdaptLoad(). 199 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or 200 HDF5 file viewer, obtained from PetscViewerHDF5Open() 201 202 Level: intermediate 203 204 Notes: 205 The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored. 206 207 Notes for advanced users: 208 Most users should not need to know the details of the binary storage 209 format, since TSAdaptLoad() and TSAdaptView() completely hide these details. 210 But for anyone who's interested, the standard binary matrix storage 211 format is 212 .vb 213 has not yet been determined 214 .ve 215 216 .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad() 217 @*/ 218 PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer) 219 { 220 PetscErrorCode ierr; 221 PetscBool isbinary; 222 char type[256]; 223 224 PetscFunctionBegin; 225 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 226 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 227 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 228 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 229 230 ierr = PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);CHKERRQ(ierr); 231 ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr); 232 if (adapt->ops->load) { 233 ierr = (*adapt->ops->load)(adapt,viewer);CHKERRQ(ierr); 234 } 235 PetscFunctionReturn(0); 236 } 237 238 PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer) 239 { 240 PetscErrorCode ierr; 241 PetscBool iascii,isbinary,isnone; 242 243 PetscFunctionBegin; 244 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 245 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);CHKERRQ(ierr);} 246 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 247 PetscCheckSameComm(adapt,1,viewer,2); 248 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 249 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 250 if (iascii) { 251 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);CHKERRQ(ierr); 252 ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);CHKERRQ(ierr); 253 if (!isnone) { 254 if (adapt->always_accept) {ierr = PetscViewerASCIIPrintf(viewer," always accepting steps\n");CHKERRQ(ierr);} 255 ierr = PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);CHKERRQ(ierr); 260 ierr = PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);CHKERRQ(ierr); 261 } 262 if (adapt->ops->view) { 263 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 264 ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr); 265 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 266 } 267 } else if (isbinary) { 268 char type[256]; 269 270 /* need to save FILE_CLASS_ID for adapt class */ 271 ierr = PetscStrncpy(type,((PetscObject)adapt)->type_name,256);CHKERRQ(ierr); 272 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 273 } else if (adapt->ops->view) { 274 ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr); 275 } 276 PetscFunctionReturn(0); 277 } 278 279 /*@ 280 TSAdaptReset - Resets a TSAdapt context. 281 282 Collective on TS 283 284 Input Parameter: 285 . adapt - the TSAdapt context obtained from TSAdaptCreate() 286 287 Level: developer 288 289 .seealso: TSAdaptCreate(), TSAdaptDestroy() 290 @*/ 291 PetscErrorCode TSAdaptReset(TSAdapt adapt) 292 { 293 PetscErrorCode ierr; 294 295 PetscFunctionBegin; 296 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 297 if (adapt->ops->reset) {ierr = (*adapt->ops->reset)(adapt);CHKERRQ(ierr);} 298 PetscFunctionReturn(0); 299 } 300 301 PetscErrorCode TSAdaptDestroy(TSAdapt *adapt) 302 { 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 if (!*adapt) PetscFunctionReturn(0); 307 PetscValidHeaderSpecific(*adapt,TSADAPT_CLASSID,1); 308 if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; PetscFunctionReturn(0);} 309 310 ierr = TSAdaptReset(*adapt);CHKERRQ(ierr); 311 312 if ((*adapt)->ops->destroy) {ierr = (*(*adapt)->ops->destroy)(*adapt);CHKERRQ(ierr);} 313 ierr = PetscViewerDestroy(&(*adapt)->monitor);CHKERRQ(ierr); 314 ierr = PetscHeaderDestroy(adapt);CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 /*@ 319 TSAdaptSetMonitor - Monitor the choices made by the adaptive controller 320 321 Collective on TSAdapt 322 323 Input Arguments: 324 + adapt - adaptive controller context 325 - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable 326 327 Options Database Keys: 328 . -ts_adapt_monitor 329 330 Level: intermediate 331 332 .seealso: TSAdaptChoose() 333 @*/ 334 PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg) 335 { 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 340 PetscValidLogicalCollectiveBool(adapt,flg,2); 341 if (flg) { 342 if (!adapt->monitor) {ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);CHKERRQ(ierr);} 343 } else { 344 ierr = PetscViewerDestroy(&adapt->monitor);CHKERRQ(ierr); 345 } 346 PetscFunctionReturn(0); 347 } 348 349 /*@C 350 TSAdaptSetCheckStage - Set a callback to check convergence for a stage 351 352 Logically collective on TSAdapt 353 354 Input Arguments: 355 + adapt - adaptive controller context 356 - func - stage check function 357 358 Arguments of func: 359 $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept) 360 361 + adapt - adaptive controller context 362 . ts - time stepping context 363 - accept - pending choice of whether to accept, can be modified by this routine 364 365 Level: advanced 366 367 .seealso: TSAdaptChoose() 368 @*/ 369 PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*)) 370 { 371 372 PetscFunctionBegin; 373 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 374 adapt->checkstage = func; 375 PetscFunctionReturn(0); 376 } 377 378 /*@ 379 TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of 380 any error or stability condition not meeting the prescribed goal. 381 382 Logically collective on TSAdapt 383 384 Input Arguments: 385 + adapt - time step adaptivity context, usually gotten with TSGetAdapt() 386 - flag - whether to always accept steps 387 388 Options Database Keys: 389 . -ts_adapt_always_accept 390 391 Level: intermediate 392 393 .seealso: TSAdapt, TSAdaptChoose() 394 @*/ 395 PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag) 396 { 397 PetscFunctionBegin; 398 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 399 PetscValidLogicalCollectiveBool(adapt,flag,2); 400 adapt->always_accept = flag; 401 PetscFunctionReturn(0); 402 } 403 404 /*@ 405 TSAdaptSetSafety - Set safety factors 406 407 Logically collective on TSAdapt 408 409 Input Arguments: 410 + adapt - adaptive controller context 411 . safety - safety factor relative to target error/stability goal 412 + reject_safety - extra safety factor to apply if the last step was rejected 413 414 Options Database Keys: 415 + -ts_adapt_safety 416 - -ts_adapt_reject_safety 417 418 Level: intermediate 419 420 .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose() 421 @*/ 422 PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety) 423 { 424 PetscFunctionBegin; 425 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 426 PetscValidLogicalCollectiveReal(adapt,safety,2); 427 PetscValidLogicalCollectiveReal(adapt,reject_safety,3); 428 if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety); 429 if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety); 430 if (reject_safety != PETSC_DEFAULT && reject_safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be non negative",(double)reject_safety); 431 if (reject_safety != PETSC_DEFAULT && reject_safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be less than one",(double)reject_safety); 432 if (safety != PETSC_DEFAULT) adapt->safety = safety; 433 if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety; 434 PetscFunctionReturn(0); 435 } 436 437 /*@ 438 TSAdaptGetSafety - Get safety factors 439 440 Not Collective 441 442 Input Arguments: 443 . adapt - adaptive controller context 444 445 Ouput Arguments: 446 . safety - safety factor relative to target error/stability goal 447 + reject_safety - extra safety factor to apply if the last step was rejected 448 449 Level: intermediate 450 451 .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose() 452 @*/ 453 PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety) 454 { 455 PetscFunctionBegin; 456 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 457 if (safety) PetscValidRealPointer(safety,2); 458 if (reject_safety) PetscValidRealPointer(reject_safety,3); 459 if (safety) *safety = adapt->safety; 460 if (reject_safety) *reject_safety = adapt->reject_safety; 461 PetscFunctionReturn(0); 462 } 463 464 /*@ 465 TSAdaptSetClip - Sets the admissible decrease/increase factor in step size 466 467 Logically collective on TSAdapt 468 469 Input Arguments: 470 + adapt - adaptive controller context 471 . low - admissible decrease factor 472 + high - admissible increase factor 473 474 Options Database Keys: 475 . -ts_adapt_clip 476 477 Level: intermediate 478 479 .seealso: TSAdaptChoose(), TSAdaptGetClip() 480 @*/ 481 PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high) 482 { 483 PetscFunctionBegin; 484 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 485 PetscValidLogicalCollectiveReal(adapt,low,2); 486 PetscValidLogicalCollectiveReal(adapt,high,3); 487 if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low); 488 if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low); 489 if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high); 490 if (low != PETSC_DEFAULT) adapt->clip[0] = low; 491 if (high != PETSC_DEFAULT) adapt->clip[1] = high; 492 PetscFunctionReturn(0); 493 } 494 495 /*@ 496 TSAdaptGetClip - Gets the admissible decrease/increase factor in step size 497 498 Not Collective 499 500 Input Arguments: 501 . adapt - adaptive controller context 502 503 Ouput Arguments: 504 + low - optional, admissible decrease factor 505 - high - optional, admissible increase factor 506 507 Level: intermediate 508 509 .seealso: TSAdaptChoose(), TSAdaptSetClip() 510 @*/ 511 PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high) 512 { 513 PetscFunctionBegin; 514 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 515 if (low) PetscValidRealPointer(low,2); 516 if (high) PetscValidRealPointer(high,3); 517 if (low) *low = adapt->clip[0]; 518 if (high) *high = adapt->clip[1]; 519 PetscFunctionReturn(0); 520 } 521 522 /*@ 523 TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller 524 525 Logically collective on TSAdapt 526 527 Input Arguments: 528 + adapt - time step adaptivity context, usually gotten with TSGetAdapt() 529 . hmin - minimum time step 530 - hmax - maximum time step 531 532 Options Database Keys: 533 + -ts_adapt_dt_min - minimum time step 534 - -ts_adapt_dt_max - maximum time step 535 536 Level: intermediate 537 538 .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose() 539 @*/ 540 PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax) 541 { 542 543 PetscFunctionBegin; 544 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 545 PetscValidLogicalCollectiveReal(adapt,hmin,2); 546 PetscValidLogicalCollectiveReal(adapt,hmax,3); 547 if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin); 548 if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax); 549 if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin; 550 if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax; 551 hmin = adapt->dt_min; 552 hmax = adapt->dt_max; 553 if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must geather than minimum time step %g",(double)hmax,(double)hmin); 554 PetscFunctionReturn(0); 555 } 556 557 /*@ 558 TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller 559 560 Not Collective 561 562 Input Arguments: 563 . adapt - time step adaptivity context, usually gotten with TSGetAdapt() 564 565 Output Arguments: 566 + hmin - minimum time step 567 - hmax - maximum time step 568 569 Level: intermediate 570 571 .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose() 572 @*/ 573 PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax) 574 { 575 576 PetscFunctionBegin; 577 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 578 if (hmin) PetscValidRealPointer(hmin,2); 579 if (hmax) PetscValidRealPointer(hmax,3); 580 if (hmin) *hmin = adapt->dt_min; 581 if (hmax) *hmax = adapt->dt_max; 582 PetscFunctionReturn(0); 583 } 584 585 /* 586 TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options. 587 588 Collective on TSAdapt 589 590 Input Parameter: 591 . adapt - the TSAdapt context 592 593 Options Database Keys: 594 + -ts_adapt_type <type> - algorithm to use for adaptivity 595 . -ts_adapt_always_accept - always accept steps regardless of error/stability goals 596 . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal 597 . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected 598 . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors 599 . -ts_adapt_dt_min <min> - minimum timestep to use 600 . -ts_adapt_dt_max <max> - maximum timestep to use 601 . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails 602 - -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates 603 604 Level: advanced 605 606 Notes: 607 This function is automatically called by TSSetFromOptions() 608 609 .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits() 610 611 .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(), 612 TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor() 613 */ 614 PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt) 615 { 616 PetscErrorCode ierr; 617 char type[256] = TSADAPTBASIC; 618 PetscReal safety,reject_safety,clip[2],hmin,hmax; 619 PetscBool set,flg; 620 PetscInt two; 621 622 PetscFunctionBegin; 623 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 624 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this 625 * function can only be called from inside TSSetFromOptions() */ 626 ierr = PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");CHKERRQ(ierr); 627 ierr = PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);CHKERRQ(ierr); 628 if (flg || !((PetscObject)adapt)->type_name) { 629 ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr); 630 } 631 632 ierr = PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);CHKERRQ(ierr); 633 if (set) {ierr = TSAdaptSetAlwaysAccept(adapt,flg);CHKERRQ(ierr);} 634 635 safety = adapt->safety; reject_safety = adapt->reject_safety; 636 ierr = PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);CHKERRQ(ierr); 637 ierr = PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);CHKERRQ(ierr); 638 if (set || flg) {ierr = TSAdaptSetSafety(adapt,safety,reject_safety);CHKERRQ(ierr);} 639 640 two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1]; 641 ierr = PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);CHKERRQ(ierr); 642 if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip"); 643 if (set) {ierr = TSAdaptSetClip(adapt,clip[0],clip[1]);CHKERRQ(ierr);} 644 645 hmin = adapt->dt_min; hmax = adapt->dt_max; 646 ierr = PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);CHKERRQ(ierr); 647 ierr = PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);CHKERRQ(ierr); 648 if (set || flg) {ierr = TSAdaptSetStepLimits(adapt,hmin,hmax);CHKERRQ(ierr);} 649 650 ierr = PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);CHKERRQ(ierr); 651 652 ierr = PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);CHKERRQ(ierr); 653 if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported"); 654 655 ierr = PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 656 if (set) {ierr = TSAdaptSetMonitor(adapt,flg);CHKERRQ(ierr);} 657 658 if (adapt->ops->setfromoptions) {ierr = (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);CHKERRQ(ierr);} 659 ierr = PetscOptionsTail();CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 /*@ 664 TSAdaptCandidatesClear - clear any previously set candidate schemes 665 666 Logically collective on TSAdapt 667 668 Input Argument: 669 . adapt - adaptive controller 670 671 Level: developer 672 673 .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose() 674 @*/ 675 PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt) 676 { 677 PetscErrorCode ierr; 678 679 PetscFunctionBegin; 680 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 681 ierr = PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));CHKERRQ(ierr); 682 PetscFunctionReturn(0); 683 } 684 685 /*@C 686 TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from 687 688 Logically collective on TSAdapt 689 690 Input Arguments: 691 + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate() 692 . name - name of the candidate scheme to add 693 . order - order of the candidate scheme 694 . stageorder - stage order of the candidate scheme 695 . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints 696 . cost - relative measure of the amount of work required for the candidate scheme 697 - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme 698 699 Note: 700 This routine is not available in Fortran. 701 702 Level: developer 703 704 .seealso: TSAdaptCandidatesClear(), TSAdaptChoose() 705 @*/ 706 PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse) 707 { 708 PetscInt c; 709 710 PetscFunctionBegin; 711 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 712 if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order); 713 if (inuse) { 714 if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()"); 715 adapt->candidates.inuse_set = PETSC_TRUE; 716 } 717 /* first slot if this is the current scheme, otherwise the next available slot */ 718 c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n; 719 720 adapt->candidates.name[c] = name; 721 adapt->candidates.order[c] = order; 722 adapt->candidates.stageorder[c] = stageorder; 723 adapt->candidates.ccfl[c] = ccfl; 724 adapt->candidates.cost[c] = cost; 725 adapt->candidates.n++; 726 PetscFunctionReturn(0); 727 } 728 729 /*@C 730 TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost 731 732 Not Collective 733 734 Input Arguments: 735 . adapt - time step adaptivity context 736 737 Output Arguments: 738 + n - number of candidate schemes, always at least 1 739 . order - the order of each candidate scheme 740 . stageorder - the stage order of each candidate scheme 741 . ccfl - the CFL coefficient of each scheme 742 - cost - the relative cost of each scheme 743 744 Level: developer 745 746 Note: 747 The current scheme is always returned in the first slot 748 749 .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose() 750 @*/ 751 PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost) 752 { 753 PetscFunctionBegin; 754 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 755 if (n) *n = adapt->candidates.n; 756 if (order) *order = adapt->candidates.order; 757 if (stageorder) *stageorder = adapt->candidates.stageorder; 758 if (ccfl) *ccfl = adapt->candidates.ccfl; 759 if (cost) *cost = adapt->candidates.cost; 760 PetscFunctionReturn(0); 761 } 762 763 /*@C 764 TSAdaptChoose - choose which method and step size to use for the next step 765 766 Collective on TSAdapt 767 768 Input Arguments: 769 + adapt - adaptive contoller 770 - h - current step size 771 772 Output Arguments: 773 + next_sc - optional, scheme to use for the next step 774 . next_h - step size to use for the next step 775 - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size 776 777 Note: 778 The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is 779 being retried after an initial rejection. 780 781 Level: developer 782 783 .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd() 784 @*/ 785 PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept) 786 { 787 PetscErrorCode ierr; 788 PetscInt ncandidates = adapt->candidates.n; 789 PetscInt scheme = 0; 790 PetscReal wlte = -1.0; 791 PetscReal wltea = -1.0; 792 PetscReal wlter = -1.0; 793 794 PetscFunctionBegin; 795 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 796 PetscValidHeaderSpecific(ts,TS_CLASSID,2); 797 if (next_sc) PetscValidIntPointer(next_sc,4); 798 PetscValidPointer(next_h,5); 799 PetscValidIntPointer(accept,6); 800 if (next_sc) *next_sc = 0; 801 802 /* Do not mess with adaptivity while handling events*/ 803 if (ts->event && ts->event->status != TSEVENT_NONE) { 804 *next_h = h; 805 *accept = PETSC_TRUE; 806 PetscFunctionReturn(0); 807 } 808 809 ierr = (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);CHKERRQ(ierr); 810 if (scheme < 0 || (ncandidates > 0 && scheme >= ncandidates)) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",scheme,ncandidates-1); 811 if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h); 812 if (next_sc) *next_sc = scheme; 813 814 if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 815 /* Reduce time step if it overshoots max time */ 816 if (ts->ptime + ts->time_step + *next_h >= ts->max_time) { 817 PetscReal next_dt = ts->max_time - (ts->ptime + ts->time_step); 818 if (next_dt > PETSC_SMALL) *next_h = next_dt; 819 else ts->reason = TS_CONVERGED_TIME; 820 } 821 } 822 823 if (adapt->monitor) { 824 const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : ""; 825 ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 826 if (wlte < 0) { 827 ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h);CHKERRQ(ierr); 828 } else { 829 ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e wlte=%5.3g wltea=%5.3g wlter=%5.3g\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h,(double)wlte,(double)wltea,(double)wlter);CHKERRQ(ierr); 830 } 831 ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 832 } 833 PetscFunctionReturn(0); 834 } 835 836 /*@ 837 TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails) 838 839 Collective on TSAdapt 840 841 Input Arguments: 842 + adapt - adaptive controller context 843 . ts - time stepper 844 . t - Current simulation time 845 - Y - Current solution vector 846 847 Output Arguments: 848 . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject 849 850 Level: developer 851 852 .seealso: 853 @*/ 854 PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept) 855 { 856 PetscErrorCode ierr; 857 SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING; 858 859 PetscFunctionBegin; 860 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 861 PetscValidHeaderSpecific(ts,TS_CLASSID,2); 862 PetscValidIntPointer(accept,3); 863 864 if (ts->snes) {ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);} 865 if (snesreason < 0) { 866 *accept = PETSC_FALSE; 867 if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) { 868 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 869 ierr = PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr); 870 if (adapt->monitor) { 871 ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 872 ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected t=%-11g+%10.3e, nonlinear solve failures %D greater than current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)ts->time_step,ts->num_snes_failures);CHKERRQ(ierr); 873 ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 874 } 875 } 876 } else { 877 *accept = PETSC_TRUE; 878 ierr = TSFunctionDomainError(ts,t,Y,accept);CHKERRQ(ierr); 879 if(*accept && adapt->checkstage) { 880 ierr = (*adapt->checkstage)(adapt,ts,t,Y,accept);CHKERRQ(ierr); 881 } 882 } 883 884 if(!(*accept) && !ts->reason) { 885 PetscReal dt,new_dt; 886 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 887 new_dt = dt * adapt->scale_solve_failed; 888 ierr = TSSetTimeStep(ts,new_dt);CHKERRQ(ierr); 889 if (adapt->monitor) { 890 ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 891 ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)dt,(double)new_dt);CHKERRQ(ierr); 892 ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); 893 } 894 } 895 PetscFunctionReturn(0); 896 } 897 898 /*@ 899 TSAdaptCreate - create an adaptive controller context for time stepping 900 901 Collective on MPI_Comm 902 903 Input Parameter: 904 . comm - The communicator 905 906 Output Parameter: 907 . adapt - new TSAdapt object 908 909 Level: developer 910 911 Notes: 912 TSAdapt creation is handled by TS, so users should not need to call this function. 913 914 .keywords: TSAdapt, create 915 .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy() 916 @*/ 917 PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt) 918 { 919 PetscErrorCode ierr; 920 TSAdapt adapt; 921 922 PetscFunctionBegin; 923 PetscValidPointer(inadapt,1); 924 *inadapt = NULL; 925 ierr = TSAdaptInitializePackage();CHKERRQ(ierr); 926 927 ierr = PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);CHKERRQ(ierr); 928 929 adapt->always_accept = PETSC_FALSE; 930 adapt->safety = 0.9; 931 adapt->reject_safety = 0.5; 932 adapt->clip[0] = 0.1; 933 adapt->clip[1] = 10.; 934 adapt->dt_min = 1e-20; 935 adapt->dt_max = 1e+20; 936 adapt->scale_solve_failed = 0.25; 937 adapt->wnormtype = NORM_2; 938 939 ierr = TSAdaptSetType(adapt,TSADAPTBASIC);CHKERRQ(ierr); 940 941 *inadapt = adapt; 942 PetscFunctionReturn(0); 943 } 944