1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 PetscClassId TSADAPT_CLASSID; 4 5 static PetscFunctionList TSAdaptList; 6 static PetscBool TSAdaptPackageInitialized; 7 static PetscBool TSAdaptRegisterAllCalled; 8 9 PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt); 10 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt); 11 PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt); 12 PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt); 13 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt); 14 PETSC_EXTERN PetscErrorCode TSAdaptCreate_History(TSAdapt); 15 16 /*@C 17 TSAdaptRegister - adds a TSAdapt implementation 18 19 Not Collective, No Fortran Support 20 21 Input Parameters: 22 + sname - name of user-defined adaptivity scheme 23 - function - routine to create method context 24 25 Level: advanced 26 27 Notes: 28 `TSAdaptRegister()` may be called multiple times to add several user-defined families. 29 30 Example Usage: 31 .vb 32 TSAdaptRegister("my_scheme", MySchemeCreate); 33 .ve 34 35 Then, your scheme can be chosen with the procedural interface via 36 .vb 37 TSAdaptSetType(ts, "my_scheme") 38 .ve 39 or at runtime via the option 40 .vb 41 -ts_adapt_type my_scheme 42 .ve 43 44 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdaptRegisterAll()` 45 @*/ 46 PetscErrorCode TSAdaptRegister(const char sname[], PetscErrorCode (*function)(TSAdapt)) 47 { 48 PetscFunctionBegin; 49 PetscCall(TSAdaptInitializePackage()); 50 PetscCall(PetscFunctionListAdd(&TSAdaptList, sname, function)); 51 PetscFunctionReturn(PETSC_SUCCESS); 52 } 53 54 /*@C 55 TSAdaptRegisterAll - Registers all of the adaptivity schemes in `TSAdapt` 56 57 Not Collective 58 59 Level: advanced 60 61 .seealso: [](ch_ts), `TSAdaptRegisterDestroy()` 62 @*/ 63 PetscErrorCode TSAdaptRegisterAll(void) 64 { 65 PetscFunctionBegin; 66 if (TSAdaptRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 67 TSAdaptRegisterAllCalled = PETSC_TRUE; 68 PetscCall(TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None)); 69 PetscCall(TSAdaptRegister(TSADAPTBASIC, TSAdaptCreate_Basic)); 70 PetscCall(TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP)); 71 PetscCall(TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL)); 72 PetscCall(TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE)); 73 PetscCall(TSAdaptRegister(TSADAPTHISTORY, TSAdaptCreate_History)); 74 PetscFunctionReturn(PETSC_SUCCESS); 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 .seealso: [](ch_ts), `PetscFinalize()` 84 @*/ 85 PetscErrorCode TSAdaptFinalizePackage(void) 86 { 87 PetscFunctionBegin; 88 PetscCall(PetscFunctionListDestroy(&TSAdaptList)); 89 TSAdaptPackageInitialized = PETSC_FALSE; 90 TSAdaptRegisterAllCalled = PETSC_FALSE; 91 PetscFunctionReturn(PETSC_SUCCESS); 92 } 93 94 /*@C 95 TSAdaptInitializePackage - This function initializes everything in the `TSAdapt` package. It is 96 called from `TSInitializePackage()`. 97 98 Level: developer 99 100 .seealso: [](ch_ts), `PetscInitialize()` 101 @*/ 102 PetscErrorCode TSAdaptInitializePackage(void) 103 { 104 PetscFunctionBegin; 105 if (TSAdaptPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 106 TSAdaptPackageInitialized = PETSC_TRUE; 107 PetscCall(PetscClassIdRegister("TSAdapt", &TSADAPT_CLASSID)); 108 PetscCall(TSAdaptRegisterAll()); 109 PetscCall(PetscRegisterFinalize(TSAdaptFinalizePackage)); 110 PetscFunctionReturn(PETSC_SUCCESS); 111 } 112 113 /*@ 114 TSAdaptSetType - sets the approach used for the error adapter 115 116 Logicially Collective 117 118 Input Parameters: 119 + adapt - the `TS` adapter, most likely obtained with `TSGetAdapt()` 120 - type - one of the `TSAdaptType` 121 122 Options Database Key: 123 . -ts_adapt_type <basic or dsp or none> - to set the adapter type 124 125 Level: intermediate 126 127 .seealso: [](ch_ts), [](sec_ts_error_control), `TSGetAdapt()`, `TSAdaptDestroy()`, `TSAdaptType`, `TSAdaptGetType()` 128 @*/ 129 PetscErrorCode TSAdaptSetType(TSAdapt adapt, TSAdaptType type) 130 { 131 PetscBool match; 132 PetscErrorCode (*r)(TSAdapt); 133 134 PetscFunctionBegin; 135 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 136 PetscAssertPointer(type, 2); 137 PetscCall(PetscObjectTypeCompare((PetscObject)adapt, type, &match)); 138 if (match) PetscFunctionReturn(PETSC_SUCCESS); 139 PetscCall(PetscFunctionListFind(TSAdaptList, type, &r)); 140 PetscCheck(r, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSAdapt type \"%s\" given", type); 141 PetscTryTypeMethod(adapt, destroy); 142 PetscCall(PetscMemzero(adapt->ops, sizeof(struct _TSAdaptOps))); 143 PetscCall(PetscObjectChangeTypeName((PetscObject)adapt, type)); 144 PetscCall((*r)(adapt)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 /*@ 149 TSAdaptGetType - gets the `TS` adapter method type (as a string). 150 151 Not Collective 152 153 Input Parameter: 154 . adapt - The `TS` adapter, most likely obtained with `TSGetAdapt()` 155 156 Output Parameter: 157 . type - The name of `TS` adapter method 158 159 Level: intermediate 160 161 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptType`, `TSAdaptSetType()` 162 @*/ 163 PetscErrorCode TSAdaptGetType(TSAdapt adapt, TSAdaptType *type) 164 { 165 PetscFunctionBegin; 166 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 167 PetscAssertPointer(type, 2); 168 *type = ((PetscObject)adapt)->type_name; 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt, const char prefix[]) 173 { 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 176 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 /*@ 181 TSAdaptLoad - Loads a TSAdapt that has been stored in binary with `TSAdaptView()`. 182 183 Collective 184 185 Input Parameters: 186 + adapt - the newly loaded `TSAdapt`, this needs to have been created with `TSAdaptCreate()` or 187 some related function before a call to `TSAdaptLoad()`. 188 - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or 189 HDF5 file viewer, obtained from `PetscViewerHDF5Open()` 190 191 Level: intermediate 192 193 Note: 194 The type is determined by the data in the file, any type set into the `TSAdapt` before this call is ignored. 195 196 .seealso: [](ch_ts), `PetscViewerBinaryOpen()`, `TSAdaptView()`, `MatLoad()`, `VecLoad()`, `TSAdapt` 197 @*/ 198 PetscErrorCode TSAdaptLoad(TSAdapt adapt, PetscViewer viewer) 199 { 200 PetscBool isbinary; 201 char type[256]; 202 203 PetscFunctionBegin; 204 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 205 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 206 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 207 PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 208 209 PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR)); 210 PetscCall(TSAdaptSetType(adapt, type)); 211 PetscTryTypeMethod(adapt, load, viewer); 212 PetscFunctionReturn(PETSC_SUCCESS); 213 } 214 215 /*@ 216 TSAdaptView - Prints the `TSAdapt` data structure. 217 218 Collective 219 220 Input Parameters: 221 + adapt - the `TSAdapt` context obtained from `TSGetAdapt()` 222 - viewer - visualization context 223 224 Options Database Key: 225 . -ts_view - calls `TSView()` at end of `TSStep()` 226 227 Level: advanced 228 229 Notes: 230 This is called by `TSView()` so rarely called directly. 231 232 The available visualization contexts include 233 + `PETSC_VIEWER_STDOUT_SELF` - standard output (default) 234 - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard 235 output where only the first processor opens 236 the file. All other processes send their 237 data to the first process to print. 238 239 The user can open an alternative visualization context with 240 `PetscViewerASCIIOpen()` - output to a specified file. 241 242 In the debugger you can do call `TSAdaptView`(adapt,0) to display the `TSAdapt`. (The same holds for any PETSc object viewer). 243 244 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSView()`, `PetscViewer`, `PetscViewerASCIIOpen()` 245 @*/ 246 PetscErrorCode TSAdaptView(TSAdapt adapt, PetscViewer viewer) 247 { 248 PetscBool isascii, isbinary, isnone, isglee; 249 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 252 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt), &viewer)); 253 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 254 PetscCheckSameComm(adapt, 1, viewer, 2); 255 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 256 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 257 if (isascii) { 258 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer)); 259 PetscCall(PetscObjectTypeCompare((PetscObject)adapt, TSADAPTNONE, &isnone)); 260 PetscCall(PetscObjectTypeCompare((PetscObject)adapt, TSADAPTGLEE, &isglee)); 261 if (!isnone) { 262 if (adapt->always_accept) PetscCall(PetscViewerASCIIPrintf(viewer, " always accepting steps\n")); 263 PetscCall(PetscViewerASCIIPrintf(viewer, " safety factor %g\n", (double)adapt->safety)); 264 PetscCall(PetscViewerASCIIPrintf(viewer, " extra safety factor after step rejection %g\n", (double)adapt->reject_safety)); 265 PetscCall(PetscViewerASCIIPrintf(viewer, " clip fastest increase %g\n", (double)adapt->clip[1])); 266 PetscCall(PetscViewerASCIIPrintf(viewer, " clip fastest decrease %g\n", (double)adapt->clip[0])); 267 PetscCall(PetscViewerASCIIPrintf(viewer, " maximum allowed timestep %g\n", (double)adapt->dt_max)); 268 PetscCall(PetscViewerASCIIPrintf(viewer, " minimum allowed timestep %g\n", (double)adapt->dt_min)); 269 PetscCall(PetscViewerASCIIPrintf(viewer, " maximum solution absolute value to be ignored %g\n", (double)adapt->ignore_max)); 270 } 271 if (isglee) { 272 if (adapt->glee_use_local) { 273 PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE uses local error control\n")); 274 } else { 275 PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE uses global error control\n")); 276 } 277 } 278 PetscCall(PetscViewerASCIIPushTab(viewer)); 279 PetscTryTypeMethod(adapt, view, viewer); 280 PetscCall(PetscViewerASCIIPopTab(viewer)); 281 } else if (isbinary) { 282 char type[256]; 283 284 /* need to save FILE_CLASS_ID for adapt class */ 285 PetscCall(PetscStrncpy(type, ((PetscObject)adapt)->type_name, 256)); 286 PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR)); 287 } else PetscTryTypeMethod(adapt, view, viewer); 288 PetscFunctionReturn(PETSC_SUCCESS); 289 } 290 291 /*@ 292 TSAdaptReset - Resets a `TSAdapt` context to its defaults 293 294 Collective 295 296 Input Parameter: 297 . adapt - the `TSAdapt` context obtained from `TSGetAdapt()` or `TSAdaptCreate()` 298 299 Level: developer 300 301 .seealso: [](ch_ts), [](sec_ts_error_control), `TSGetAdapt()`, `TSAdapt`, `TSAdaptCreate()`, `TSAdaptDestroy()` 302 @*/ 303 PetscErrorCode TSAdaptReset(TSAdapt adapt) 304 { 305 PetscFunctionBegin; 306 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 307 PetscTryTypeMethod(adapt, reset); 308 PetscFunctionReturn(PETSC_SUCCESS); 309 } 310 311 PetscErrorCode TSAdaptDestroy(TSAdapt *adapt) 312 { 313 PetscFunctionBegin; 314 if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS); 315 PetscValidHeaderSpecific(*adapt, TSADAPT_CLASSID, 1); 316 if (--((PetscObject)*adapt)->refct > 0) { 317 *adapt = NULL; 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320 321 PetscCall(TSAdaptReset(*adapt)); 322 323 PetscTryTypeMethod(*adapt, destroy); 324 PetscCall(PetscViewerDestroy(&(*adapt)->monitor)); 325 PetscCall(PetscHeaderDestroy(adapt)); 326 PetscFunctionReturn(PETSC_SUCCESS); 327 } 328 329 /*@ 330 TSAdaptSetMonitor - Monitor the choices made by the adaptive controller 331 332 Collective 333 334 Input Parameters: 335 + adapt - adaptive controller context 336 - flg - `PETSC_TRUE` to active a monitor, `PETSC_FALSE` to disable 337 338 Options Database Key: 339 . -ts_adapt_monitor - to turn on monitoring 340 341 Level: intermediate 342 343 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSGetAdapt()`, `TSAdaptChoose()` 344 @*/ 345 PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt, PetscBool flg) 346 { 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 349 PetscValidLogicalCollectiveBool(adapt, flg, 2); 350 if (flg) { 351 if (!adapt->monitor) PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt), "stdout", &adapt->monitor)); 352 } else { 353 PetscCall(PetscViewerDestroy(&adapt->monitor)); 354 } 355 PetscFunctionReturn(PETSC_SUCCESS); 356 } 357 358 /*@C 359 TSAdaptSetCheckStage - Set a callback to check convergence for a stage 360 361 Logically Collective 362 363 Input Parameters: 364 + adapt - adaptive controller context 365 - func - stage check function 366 367 Calling sequence: 368 + adapt - adaptive controller context 369 . ts - time stepping context 370 . t - current time 371 . Y - current solution vector 372 - accept - pending choice of whether to accept, can be modified by this routine 373 374 Level: advanced 375 376 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSGetAdapt()`, `TSAdaptChoose()` 377 @*/ 378 PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt, PetscErrorCode (*func)(TSAdapt adapt, TS ts, PetscReal t, Vec Y, PetscBool *accept)) 379 { 380 PetscFunctionBegin; 381 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 382 adapt->checkstage = func; 383 PetscFunctionReturn(PETSC_SUCCESS); 384 } 385 386 /*@ 387 TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of 388 any error or stability condition not meeting the prescribed goal. 389 390 Logically Collective 391 392 Input Parameters: 393 + adapt - time step adaptivity context, usually gotten with `TSGetAdapt()` 394 - flag - whether to always accept steps 395 396 Options Database Key: 397 . -ts_adapt_always_accept - to always accept steps 398 399 Level: intermediate 400 401 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSGetAdapt()`, `TSAdaptChoose()` 402 @*/ 403 PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt, PetscBool flag) 404 { 405 PetscFunctionBegin; 406 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 407 PetscValidLogicalCollectiveBool(adapt, flag, 2); 408 adapt->always_accept = flag; 409 PetscFunctionReturn(PETSC_SUCCESS); 410 } 411 412 /*@ 413 TSAdaptSetSafety - Set safety factors for time step adaptor 414 415 Logically Collective 416 417 Input Parameters: 418 + adapt - adaptive controller context 419 . safety - safety factor relative to target error/stability goal 420 - reject_safety - extra safety factor to apply if the last step was rejected 421 422 Options Database Keys: 423 + -ts_adapt_safety <safety> - to set safety factor 424 - -ts_adapt_reject_safety <reject_safety> - to set reject safety factor 425 426 Level: intermediate 427 428 Note: 429 Use `PETSC_CURRENT` to keep the current value for either parameter 430 431 Fortran Note: 432 Use `PETSC_CURRENT_REAL` 433 434 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptGetSafety()`, `TSAdaptChoose()` 435 @*/ 436 PetscErrorCode TSAdaptSetSafety(TSAdapt adapt, PetscReal safety, PetscReal reject_safety) 437 { 438 PetscFunctionBegin; 439 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 440 PetscValidLogicalCollectiveReal(adapt, safety, 2); 441 PetscValidLogicalCollectiveReal(adapt, reject_safety, 3); 442 PetscCheck(safety == (PetscReal)PETSC_CURRENT || safety >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Safety factor %g must be non negative", (double)safety); 443 PetscCheck(safety == (PetscReal)PETSC_CURRENT || safety <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Safety factor %g must be less than one", (double)safety); 444 PetscCheck(reject_safety == (PetscReal)PETSC_CURRENT || reject_safety >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reject safety factor %g must be non negative", (double)reject_safety); 445 PetscCheck(reject_safety == (PetscReal)PETSC_CURRENT || reject_safety <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reject safety factor %g must be less than one", (double)reject_safety); 446 if (safety != (PetscReal)PETSC_CURRENT) adapt->safety = safety; 447 if (reject_safety != (PetscReal)PETSC_CURRENT) adapt->reject_safety = reject_safety; 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 /*@ 452 TSAdaptGetSafety - Get safety factors for time step adapter 453 454 Not Collective 455 456 Input Parameter: 457 . adapt - adaptive controller context 458 459 Output Parameters: 460 + safety - safety factor relative to target error/stability goal 461 - reject_safety - extra safety factor to apply if the last step was rejected 462 463 Level: intermediate 464 465 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptSetSafety()`, `TSAdaptChoose()` 466 @*/ 467 PetscErrorCode TSAdaptGetSafety(TSAdapt adapt, PetscReal *safety, PetscReal *reject_safety) 468 { 469 PetscFunctionBegin; 470 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 471 if (safety) PetscAssertPointer(safety, 2); 472 if (reject_safety) PetscAssertPointer(reject_safety, 3); 473 if (safety) *safety = adapt->safety; 474 if (reject_safety) *reject_safety = adapt->reject_safety; 475 PetscFunctionReturn(PETSC_SUCCESS); 476 } 477 478 /*@ 479 TSAdaptSetMaxIgnore - Set error estimation threshold. Solution components below this threshold value will not be considered when computing error norms 480 for time step adaptivity (in absolute value). A negative value (default) of the threshold leads to considering all solution components. 481 482 Logically Collective 483 484 Input Parameters: 485 + adapt - adaptive controller context 486 - max_ignore - threshold for solution components that are ignored during error estimation 487 488 Options Database Key: 489 . -ts_adapt_max_ignore <max_ignore> - to set the threshold 490 491 Level: intermediate 492 493 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptGetMaxIgnore()`, `TSAdaptChoose()` 494 @*/ 495 PetscErrorCode TSAdaptSetMaxIgnore(TSAdapt adapt, PetscReal max_ignore) 496 { 497 PetscFunctionBegin; 498 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 499 PetscValidLogicalCollectiveReal(adapt, max_ignore, 2); 500 adapt->ignore_max = max_ignore; 501 PetscFunctionReturn(PETSC_SUCCESS); 502 } 503 504 /*@ 505 TSAdaptGetMaxIgnore - Get error estimation threshold. Solution components below this threshold value will not be considered when computing error norms 506 for time step adaptivity (in absolute value). 507 508 Not Collective 509 510 Input Parameter: 511 . adapt - adaptive controller context 512 513 Output Parameter: 514 . max_ignore - threshold for solution components that are ignored during error estimation 515 516 Level: intermediate 517 518 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptSetMaxIgnore()`, `TSAdaptChoose()` 519 @*/ 520 PetscErrorCode TSAdaptGetMaxIgnore(TSAdapt adapt, PetscReal *max_ignore) 521 { 522 PetscFunctionBegin; 523 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 524 PetscAssertPointer(max_ignore, 2); 525 *max_ignore = adapt->ignore_max; 526 PetscFunctionReturn(PETSC_SUCCESS); 527 } 528 529 /*@ 530 TSAdaptSetClip - Sets the admissible decrease/increase factor in step size in the time step adapter 531 532 Logically collective 533 534 Input Parameters: 535 + adapt - adaptive controller context 536 . low - admissible decrease factor 537 - high - admissible increase factor 538 539 Options Database Key: 540 . -ts_adapt_clip <low>,<high> - to set admissible time step decrease and increase factors 541 542 Level: intermediate 543 544 Note: 545 Use `PETSC_CURRENT` to keep the current value for either parameter 546 547 Fortran Note: 548 Use `PETSC_CURRENT_REAL` 549 550 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptChoose()`, `TSAdaptGetClip()`, `TSAdaptSetScaleSolveFailed()` 551 @*/ 552 PetscErrorCode TSAdaptSetClip(TSAdapt adapt, PetscReal low, PetscReal high) 553 { 554 PetscFunctionBegin; 555 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 556 PetscValidLogicalCollectiveReal(adapt, low, 2); 557 PetscValidLogicalCollectiveReal(adapt, high, 3); 558 PetscCheck(low == (PetscReal)PETSC_CURRENT || low >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Decrease factor %g must be non negative", (double)low); 559 PetscCheck(low == (PetscReal)PETSC_CURRENT || low <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Decrease factor %g must be less than one", (double)low); 560 PetscCheck(high == (PetscReal)PETSC_CURRENT || high >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Increase factor %g must be greater than one", (double)high); 561 if (low != (PetscReal)PETSC_CURRENT) adapt->clip[0] = low; 562 if (high != (PetscReal)PETSC_CURRENT) adapt->clip[1] = high; 563 PetscFunctionReturn(PETSC_SUCCESS); 564 } 565 566 /*@ 567 TSAdaptGetClip - Gets the admissible decrease/increase factor in step size in the time step adapter 568 569 Not Collective 570 571 Input Parameter: 572 . adapt - adaptive controller context 573 574 Output Parameters: 575 + low - optional, admissible decrease factor 576 - high - optional, admissible increase factor 577 578 Level: intermediate 579 580 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptChoose()`, `TSAdaptSetClip()`, `TSAdaptSetScaleSolveFailed()` 581 @*/ 582 PetscErrorCode TSAdaptGetClip(TSAdapt adapt, PetscReal *low, PetscReal *high) 583 { 584 PetscFunctionBegin; 585 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 586 if (low) PetscAssertPointer(low, 2); 587 if (high) PetscAssertPointer(high, 3); 588 if (low) *low = adapt->clip[0]; 589 if (high) *high = adapt->clip[1]; 590 PetscFunctionReturn(PETSC_SUCCESS); 591 } 592 593 /*@ 594 TSAdaptSetScaleSolveFailed - Scale step size by this factor if solve fails 595 596 Logically Collective 597 598 Input Parameters: 599 + adapt - adaptive controller context 600 - scale - scale 601 602 Options Database Key: 603 . -ts_adapt_scale_solve_failed <scale> - to set scale step by this factor if solve fails 604 605 Level: intermediate 606 607 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptChoose()`, `TSAdaptGetScaleSolveFailed()`, `TSAdaptGetClip()` 608 @*/ 609 PetscErrorCode TSAdaptSetScaleSolveFailed(TSAdapt adapt, PetscReal scale) 610 { 611 PetscFunctionBegin; 612 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 613 PetscValidLogicalCollectiveReal(adapt, scale, 2); 614 PetscCheck(scale > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scale factor %g must be positive", (double)scale); 615 PetscCheck(scale <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scale factor %g must be less than one", (double)scale); 616 adapt->scale_solve_failed = scale; 617 PetscFunctionReturn(PETSC_SUCCESS); 618 } 619 620 /*@ 621 TSAdaptGetScaleSolveFailed - Gets the admissible decrease/increase factor in step size 622 623 Not Collective 624 625 Input Parameter: 626 . adapt - adaptive controller context 627 628 Output Parameter: 629 . scale - scale factor 630 631 Level: intermediate 632 633 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptChoose()`, `TSAdaptSetScaleSolveFailed()`, `TSAdaptSetClip()` 634 @*/ 635 PetscErrorCode TSAdaptGetScaleSolveFailed(TSAdapt adapt, PetscReal *scale) 636 { 637 PetscFunctionBegin; 638 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 639 if (scale) PetscAssertPointer(scale, 2); 640 if (scale) *scale = adapt->scale_solve_failed; 641 PetscFunctionReturn(PETSC_SUCCESS); 642 } 643 644 /*@ 645 TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the time step controller 646 647 Logically Collective 648 649 Input Parameters: 650 + adapt - time step adaptivity context, usually gotten with `TSGetAdapt()` 651 . hmin - minimum time step 652 - hmax - maximum time step 653 654 Options Database Keys: 655 + -ts_adapt_dt_min <min> - to set minimum time step 656 - -ts_adapt_dt_max <max> - to set maximum time step 657 658 Level: intermediate 659 660 Note: 661 Use `PETSC_CURRENT` to keep the current value for either parameter 662 663 Fortran Note: 664 Use `PETSC_CURRENT_REAL` 665 666 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptGetStepLimits()`, `TSAdaptChoose()` 667 @*/ 668 PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt, PetscReal hmin, PetscReal hmax) 669 { 670 PetscFunctionBegin; 671 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 672 PetscValidLogicalCollectiveReal(adapt, hmin, 2); 673 PetscValidLogicalCollectiveReal(adapt, hmax, 3); 674 PetscCheck(hmin == (PetscReal)PETSC_CURRENT || hmin >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum time step %g must be non negative", (double)hmin); 675 PetscCheck(hmax == (PetscReal)PETSC_CURRENT || hmax >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum time step %g must be non negative", (double)hmax); 676 if (hmin != (PetscReal)PETSC_CURRENT) adapt->dt_min = hmin; 677 if (hmax != (PetscReal)PETSC_CURRENT) adapt->dt_max = hmax; 678 hmin = adapt->dt_min; 679 hmax = adapt->dt_max; 680 PetscCheck(hmax > hmin, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum time step %g must greater than minimum time step %g", (double)hmax, (double)hmin); 681 PetscFunctionReturn(PETSC_SUCCESS); 682 } 683 684 /*@ 685 TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the time step controller 686 687 Not Collective 688 689 Input Parameter: 690 . adapt - time step adaptivity context, usually gotten with `TSGetAdapt()` 691 692 Output Parameters: 693 + hmin - minimum time step 694 - hmax - maximum time step 695 696 Level: intermediate 697 698 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptSetStepLimits()`, `TSAdaptChoose()` 699 @*/ 700 PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt, PetscReal *hmin, PetscReal *hmax) 701 { 702 PetscFunctionBegin; 703 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 704 if (hmin) PetscAssertPointer(hmin, 2); 705 if (hmax) PetscAssertPointer(hmax, 3); 706 if (hmin) *hmin = adapt->dt_min; 707 if (hmax) *hmax = adapt->dt_max; 708 PetscFunctionReturn(PETSC_SUCCESS); 709 } 710 711 /*@C 712 TSAdaptSetFromOptions - Sets various `TSAdapt` parameters from user options. 713 714 Collective 715 716 Input Parameters: 717 + adapt - the `TSAdapt` context 718 - PetscOptionsObject - object created by `PetscOptionsBegin()` 719 720 Options Database Keys: 721 + -ts_adapt_type <type> - algorithm to use for adaptivity 722 . -ts_adapt_always_accept - always accept steps regardless of error/stability goals 723 . -ts_adapt_safety <safety> - safety factor relative to target error/stability goal 724 . -ts_adapt_reject_safety <safety> - extra safety factor to apply if the last step was rejected 725 . -ts_adapt_clip <low,high> - admissible time step decrease and increase factors 726 . -ts_adapt_dt_min <min> - minimum timestep to use 727 . -ts_adapt_dt_max <max> - maximum timestep to use 728 . -ts_adapt_scale_solve_failed <scale> - scale timestep by this factor if a solve fails 729 . -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates 730 - -ts_adapt_time_step_increase_delay - number of timesteps to delay increasing the time step after it has been decreased due to failed solver 731 732 Level: advanced 733 734 Note: 735 This function is automatically called by `TSSetFromOptions()` 736 737 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSGetAdapt()`, `TSAdaptSetType()`, `TSAdaptSetAlwaysAccept()`, `TSAdaptSetSafety()`, 738 `TSAdaptSetClip()`, `TSAdaptSetScaleSolveFailed()`, `TSAdaptSetStepLimits()`, `TSAdaptSetMonitor()` 739 @*/ 740 PetscErrorCode TSAdaptSetFromOptions(TSAdapt adapt, PetscOptionItems PetscOptionsObject) 741 { 742 char type[256] = TSADAPTBASIC; 743 PetscReal safety, reject_safety, clip[2], scale, hmin, hmax; 744 PetscBool set, flg; 745 PetscInt two; 746 747 PetscFunctionBegin; 748 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 749 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this 750 * function can only be called from inside TSSetFromOptions() */ 751 PetscOptionsHeadBegin(PetscOptionsObject, "TS Adaptivity options"); 752 PetscCall(PetscOptionsFList("-ts_adapt_type", "Algorithm to use for adaptivity", "TSAdaptSetType", TSAdaptList, ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type, type, sizeof(type), &flg)); 753 if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, type)); 754 755 PetscCall(PetscOptionsBool("-ts_adapt_always_accept", "Always accept the step", "TSAdaptSetAlwaysAccept", adapt->always_accept, &flg, &set)); 756 if (set) PetscCall(TSAdaptSetAlwaysAccept(adapt, flg)); 757 758 safety = adapt->safety; 759 reject_safety = adapt->reject_safety; 760 PetscCall(PetscOptionsReal("-ts_adapt_safety", "Safety factor relative to target error/stability goal", "TSAdaptSetSafety", safety, &safety, &set)); 761 PetscCall(PetscOptionsReal("-ts_adapt_reject_safety", "Extra safety factor to apply if the last step was rejected", "TSAdaptSetSafety", reject_safety, &reject_safety, &flg)); 762 if (set || flg) PetscCall(TSAdaptSetSafety(adapt, safety, reject_safety)); 763 764 two = 2; 765 clip[0] = adapt->clip[0]; 766 clip[1] = adapt->clip[1]; 767 PetscCall(PetscOptionsRealArray("-ts_adapt_clip", "Admissible decrease/increase factor in step size", "TSAdaptSetClip", clip, &two, &set)); 768 PetscCheck(!set || (two == 2), PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Must give exactly two values to -ts_adapt_clip"); 769 if (set) PetscCall(TSAdaptSetClip(adapt, clip[0], clip[1])); 770 771 hmin = adapt->dt_min; 772 hmax = adapt->dt_max; 773 PetscCall(PetscOptionsReal("-ts_adapt_dt_min", "Minimum time step considered", "TSAdaptSetStepLimits", hmin, &hmin, &set)); 774 PetscCall(PetscOptionsReal("-ts_adapt_dt_max", "Maximum time step considered", "TSAdaptSetStepLimits", hmax, &hmax, &flg)); 775 if (set || flg) PetscCall(TSAdaptSetStepLimits(adapt, hmin, hmax)); 776 777 PetscCall(PetscOptionsReal("-ts_adapt_max_ignore", "Adaptor ignores (absolute) solution values smaller than this value", "", adapt->ignore_max, &adapt->ignore_max, &set)); 778 PetscCall(PetscOptionsBool("-ts_adapt_glee_use_local", "GLEE adaptor uses local error estimation for step control", "", adapt->glee_use_local, &adapt->glee_use_local, &set)); 779 780 PetscCall(PetscOptionsReal("-ts_adapt_scale_solve_failed", "Scale step by this factor if solve fails", "TSAdaptSetScaleSolveFailed", adapt->scale_solve_failed, &scale, &set)); 781 if (set) PetscCall(TSAdaptSetScaleSolveFailed(adapt, scale)); 782 783 PetscCall(PetscOptionsEnum("-ts_adapt_wnormtype", "Type of norm computed for error estimation", "", NormTypes, (PetscEnum)adapt->wnormtype, (PetscEnum *)&adapt->wnormtype, NULL)); 784 PetscCheck(adapt->wnormtype == NORM_2 || adapt->wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)adapt), PETSC_ERR_SUP, "Only 2-norm and infinite norm supported"); 785 786 PetscCall(PetscOptionsInt("-ts_adapt_time_step_increase_delay", "Number of timesteps to delay increasing the time step after it has been decreased due to failed solver", "TSAdaptSetTimeStepIncreaseDelay", adapt->timestepjustdecreased_delay, &adapt->timestepjustdecreased_delay, NULL)); 787 788 PetscCall(PetscOptionsBool("-ts_adapt_monitor", "Print choices made by adaptive controller", "TSAdaptSetMonitor", adapt->monitor ? PETSC_TRUE : PETSC_FALSE, &flg, &set)); 789 if (set) PetscCall(TSAdaptSetMonitor(adapt, flg)); 790 791 PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject); 792 PetscOptionsHeadEnd(); 793 PetscFunctionReturn(PETSC_SUCCESS); 794 } 795 796 /*@ 797 TSAdaptCandidatesClear - clear any previously set candidate schemes 798 799 Logically Collective 800 801 Input Parameter: 802 . adapt - adaptive controller 803 804 Level: developer 805 806 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptCreate()`, `TSAdaptCandidateAdd()`, `TSAdaptChoose()` 807 @*/ 808 PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt) 809 { 810 PetscFunctionBegin; 811 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 812 PetscCall(PetscMemzero(&adapt->candidates, sizeof(adapt->candidates))); 813 PetscFunctionReturn(PETSC_SUCCESS); 814 } 815 816 /*@C 817 TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from 818 819 Logically Collective; No Fortran Support 820 821 Input Parameters: 822 + adapt - time step adaptivity context, obtained with `TSGetAdapt()` or `TSAdaptCreate()` 823 . name - name of the candidate scheme to add 824 . order - order of the candidate scheme 825 . stageorder - stage order of the candidate scheme 826 . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints 827 . cost - relative measure of the amount of work required for the candidate scheme 828 - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme 829 830 Level: developer 831 832 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptCandidatesClear()`, `TSAdaptChoose()` 833 @*/ 834 PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt, const char name[], PetscInt order, PetscInt stageorder, PetscReal ccfl, PetscReal cost, PetscBool inuse) 835 { 836 PetscInt c; 837 838 PetscFunctionBegin; 839 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 840 PetscCheck(order >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Classical order %" PetscInt_FMT " must be a positive integer", order); 841 if (inuse) { 842 PetscCheck(!adapt->candidates.inuse_set, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()"); 843 adapt->candidates.inuse_set = PETSC_TRUE; 844 } 845 /* first slot if this is the current scheme, otherwise the next available slot */ 846 c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n; 847 848 adapt->candidates.name[c] = name; 849 adapt->candidates.order[c] = order; 850 adapt->candidates.stageorder[c] = stageorder; 851 adapt->candidates.ccfl[c] = ccfl; 852 adapt->candidates.cost[c] = cost; 853 adapt->candidates.n++; 854 PetscFunctionReturn(PETSC_SUCCESS); 855 } 856 857 /*@C 858 TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost 859 860 Not Collective 861 862 Input Parameter: 863 . adapt - time step adaptivity context 864 865 Output Parameters: 866 + n - number of candidate schemes, always at least 1 867 . order - the order of each candidate scheme 868 . stageorder - the stage order of each candidate scheme 869 . ccfl - the CFL coefficient of each scheme 870 - cost - the relative cost of each scheme 871 872 Level: developer 873 874 Note: 875 The current scheme is always returned in the first slot 876 877 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptCandidatesClear()`, `TSAdaptCandidateAdd()`, `TSAdaptChoose()` 878 @*/ 879 PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt, PetscInt *n, const PetscInt **order, const PetscInt **stageorder, const PetscReal **ccfl, const PetscReal **cost) 880 { 881 PetscFunctionBegin; 882 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 883 if (n) *n = adapt->candidates.n; 884 if (order) *order = adapt->candidates.order; 885 if (stageorder) *stageorder = adapt->candidates.stageorder; 886 if (ccfl) *ccfl = adapt->candidates.ccfl; 887 if (cost) *cost = adapt->candidates.cost; 888 PetscFunctionReturn(PETSC_SUCCESS); 889 } 890 891 /*@C 892 TSAdaptChoose - choose which method and step size to use for the next step 893 894 Collective 895 896 Input Parameters: 897 + adapt - adaptive controller 898 . ts - time stepper 899 - h - current step size 900 901 Output Parameters: 902 + next_sc - optional, scheme to use for the next step 903 . next_h - step size to use for the next step 904 - accept - `PETSC_TRUE` to accept the current step, `PETSC_FALSE` to repeat the current step with the new step size 905 906 Level: developer 907 908 Note: 909 The input value of parameter accept is retained from the last time step, so it will be `PETSC_FALSE` if the step is 910 being retried after an initial rejection. 911 912 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSAdaptCandidatesClear()`, `TSAdaptCandidateAdd()` 913 @*/ 914 PetscErrorCode TSAdaptChoose(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept) 915 { 916 PetscInt ncandidates = adapt->candidates.n; 917 PetscInt scheme = 0; 918 PetscReal wlte = -1.0; 919 PetscReal wltea = -1.0; 920 PetscReal wlter = -1.0; 921 922 PetscFunctionBegin; 923 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 924 PetscValidHeaderSpecific(ts, TS_CLASSID, 2); 925 if (next_sc) PetscAssertPointer(next_sc, 4); 926 PetscAssertPointer(next_h, 5); 927 PetscAssertPointer(accept, 6); 928 if (next_sc) *next_sc = 0; 929 930 /* Do not mess with adaptivity while handling events */ 931 if (ts->event && ts->event->processing) { 932 *next_h = h; 933 *accept = PETSC_TRUE; 934 if (adapt->monitor) { 935 PetscCall(PetscViewerASCIIAddTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 936 937 if (ts->event->iterctr == 0) { 938 /* 939 An event has been found, now finalising the event processing: performing the 1st and 2nd post-event steps. 940 Entering this if-branch means both these steps (set to either PETSC_DECIDE or numerical value) are managed 941 by the event handler. In this case the 1st post-event step is always accepted, without interference of TSAdapt. 942 Note: if the 2nd post-event step is not managed by the event handler (e.g. given 1st = numerical, 2nd = PETSC_DECIDE), 943 this if-branch is not entered, and TSAdapt may reject/adjust the proposed 1st post-event step. 944 */ 945 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt does not interfere, step %3" PetscInt_FMT " accepted. Processing post-event steps: 1-st accepted just now, 2-nd yet to come\n", ts->steps)); 946 } else PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt does not interfere, step %3" PetscInt_FMT " accepted. Event handling in progress\n", ts->steps)); 947 948 PetscCall(PetscViewerASCIISubtractTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 949 } 950 PetscFunctionReturn(PETSC_SUCCESS); 951 } 952 953 PetscUseTypeMethod(adapt, choose, ts, h, &scheme, next_h, accept, &wlte, &wltea, &wlter); 954 PetscCheck(scheme >= 0 && (ncandidates <= 0 || scheme < ncandidates), PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Chosen scheme %" PetscInt_FMT " not in valid range 0..%" PetscInt_FMT, scheme, ncandidates - 1); 955 PetscCheck(*next_h >= 0, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Computed step size %g must be positive", (double)*next_h); 956 if (next_sc) *next_sc = scheme; 957 958 if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 959 /* Increase/reduce step size if end time of next step is close to or overshoots max time */ 960 PetscReal t = ts->ptime + ts->time_step, tend, tmax, h1, hmax; 961 PetscReal a = (PetscReal)(1.0 + adapt->matchstepfac[0]); 962 PetscReal b = adapt->matchstepfac[1]; 963 964 /* 965 Logic in using 'dt_span_cached': 966 1. It always overrides *next_h, except (any of): 967 a) the current step was rejected, 968 b) the adaptor proposed to decrease the next step, 969 c) the adaptor proposed *next_h > dt_span_cached. 970 2. If *next_h was adjusted by eval_times points (or the final point): 971 -- when dt_span_cached is filled (>0), it keeps its value, 972 -- when dt_span_cached is clear (==0), it gets the unadjusted version of *next_h. 973 3. If *next_h was not adjusted as in (2), dt_span_cached is cleared. 974 Note, if a combination (1.b || 1.c) && (3) takes place, this means that 975 dt_span_cached remains unused at the moment of clearing. 976 If (1.a) takes place, dt_span_cached keeps its value. 977 Also, dt_span_cached can be updated by the event handler, see tsevent.c. 978 */ 979 if (h <= *next_h && *next_h <= adapt->dt_eval_times_cached) *next_h = adapt->dt_eval_times_cached; /* try employing the cache */ 980 h1 = *next_h; 981 tend = t + h1; 982 983 if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points) { 984 PetscCheck(ts->eval_times->worktol == 0, PetscObjectComm((PetscObject)adapt), PETSC_ERR_PLIB, "Unexpected state (tspan->worktol != 0) in TSAdaptChoose()"); 985 ts->eval_times->worktol = ts->eval_times->reltol * h1 + ts->eval_times->abstol; 986 if (PetscIsCloseAtTol(t, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) /* hit a span time point */ 987 if (ts->eval_times->time_point_idx + 1 < ts->eval_times->num_time_points) tmax = ts->eval_times->time_points[ts->eval_times->time_point_idx + 1]; 988 else tmax = ts->max_time; /* hit the last span time point */ 989 else tmax = ts->eval_times->time_points[ts->eval_times->time_point_idx]; 990 } else tmax = ts->max_time; 991 tmax = PetscMin(tmax, ts->max_time); 992 hmax = tmax - t; 993 994 if (t < tmax && tend > tmax) *next_h = hmax; 995 if (t < tmax && tend < tmax && h1 * b > hmax) *next_h = hmax / 2; 996 if (t < tmax && tend < tmax && h1 * a > hmax) *next_h = hmax; 997 if (ts->eval_times && h1 != *next_h && !adapt->dt_eval_times_cached) adapt->dt_eval_times_cached = h1; /* cache the step size if it is to be changed */ 998 if (ts->eval_times && h1 == *next_h && adapt->dt_eval_times_cached) adapt->dt_eval_times_cached = 0; /* clear the cache if the step size is unchanged */ 999 } 1000 if (adapt->monitor) { 1001 const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : ""; 1002 PetscCall(PetscViewerASCIIAddTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1003 if (wlte < 0) { 1004 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt %s %s %" PetscInt_FMT ":%s step %3" PetscInt_FMT " %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", 1005 (double)ts->ptime, (double)h, (double)*next_h)); 1006 } else { 1007 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt %s %s %" PetscInt_FMT ":%s step %3" PetscInt_FMT " %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", 1008 (double)ts->ptime, (double)h, (double)*next_h, (double)wlte, (double)wltea, (double)wlter)); 1009 } 1010 PetscCall(PetscViewerASCIISubtractTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1011 } 1012 PetscFunctionReturn(PETSC_SUCCESS); 1013 } 1014 1015 /*@ 1016 TSAdaptSetTimeStepIncreaseDelay - The number of timesteps to wait after a decrease in the timestep due to failed solver 1017 before increasing the time step. 1018 1019 Logicially Collective 1020 1021 Input Parameters: 1022 + adapt - adaptive controller context 1023 - cnt - the number of timesteps 1024 1025 Options Database Key: 1026 . -ts_adapt_time_step_increase_delay cnt - number of steps to delay the increase 1027 1028 Level: advanced 1029 1030 Notes: 1031 This is to prevent an adaptor from bouncing back and forth between two nearby timesteps. The default is 0. 1032 1033 The successful use of this option is problem dependent 1034 1035 Developer Notes: 1036 There is no theory to support this option 1037 1038 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt` 1039 @*/ 1040 PetscErrorCode TSAdaptSetTimeStepIncreaseDelay(TSAdapt adapt, PetscInt cnt) 1041 { 1042 PetscFunctionBegin; 1043 adapt->timestepjustdecreased_delay = cnt; 1044 PetscFunctionReturn(PETSC_SUCCESS); 1045 } 1046 1047 /*@ 1048 TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails or solution vector is infeasible) 1049 1050 Collective 1051 1052 Input Parameters: 1053 + adapt - adaptive controller context 1054 . ts - time stepper 1055 . t - Current simulation time 1056 - Y - Current solution vector 1057 1058 Output Parameter: 1059 . accept - `PETSC_TRUE` to accept the stage, `PETSC_FALSE` to reject 1060 1061 Level: developer 1062 1063 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt` 1064 @*/ 1065 PetscErrorCode TSAdaptCheckStage(TSAdapt adapt, TS ts, PetscReal t, Vec Y, PetscBool *accept) 1066 { 1067 SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING; 1068 PetscBool func_accept, snes_div_func; 1069 1070 PetscFunctionBegin; 1071 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 1072 PetscValidHeaderSpecific(ts, TS_CLASSID, 2); 1073 PetscAssertPointer(accept, 5); 1074 1075 PetscCall(TSFunctionDomainError(ts, t, Y, &func_accept)); 1076 if (ts->snes) PetscCall(SNESGetConvergedReason(ts->snes, &snesreason)); 1077 snes_div_func = (PetscBool)(snesreason == SNES_DIVERGED_FUNCTION_DOMAIN); 1078 if (func_accept && snesreason < 0 && !snes_div_func) { 1079 *accept = PETSC_FALSE; 1080 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failure: %s\n", ts->steps, SNESConvergedReasons[snesreason])); 1081 if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures != PETSC_UNLIMITED) { 1082 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 1083 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures)); 1084 if (adapt->monitor) { 1085 PetscCall(PetscViewerASCIIAddTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1086 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt %s step %3" PetscInt_FMT " stage rejected t=%-11g+%10.3e, nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed\n", ((PetscObject)adapt)->type_name, ts->steps, 1087 (double)ts->ptime, (double)ts->time_step, ts->num_snes_failures)); 1088 PetscCall(PetscViewerASCIISubtractTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1089 } 1090 } 1091 } else { 1092 *accept = (PetscBool)(func_accept && !snes_div_func); 1093 if (*accept && adapt->checkstage) PetscCall((*adapt->checkstage)(adapt, ts, t, Y, accept)); 1094 if (!*accept) { 1095 const char *user_func = !func_accept ? "TSSetFunctionDomainError()" : "TSAdaptSetCheckStage"; 1096 const char *snes_err = "SNES invalid function domain"; 1097 const char *err_msg = snes_div_func && func_accept ? snes_err : user_func; 1098 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", solution rejected by %s\n", ts->steps, err_msg)); 1099 if (adapt->monitor) { 1100 PetscCall(PetscViewerASCIIAddTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1101 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt %s step %3" PetscInt_FMT " stage rejected by %s\n", ((PetscObject)adapt)->type_name, ts->steps, err_msg)); 1102 PetscCall(PetscViewerASCIISubtractTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1103 } 1104 } 1105 } 1106 1107 if (!*accept && !ts->reason) { 1108 PetscReal dt, new_dt; 1109 PetscCall(TSGetTimeStep(ts, &dt)); 1110 new_dt = dt * adapt->scale_solve_failed; 1111 PetscCall(TSSetTimeStep(ts, new_dt)); 1112 adapt->timestepjustdecreased += adapt->timestepjustdecreased_delay; 1113 if (adapt->monitor) { 1114 PetscCall(PetscViewerASCIIAddTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1115 PetscCall(PetscViewerASCIIPrintf(adapt->monitor, " TSAdapt %s step %3" PetscInt_FMT " stage rejected (SNES reason %s) t=%-11g+%10.3e retrying with dt=%-10.3e\n", ((PetscObject)adapt)->type_name, ts->steps, SNESConvergedReasons[snesreason], 1116 (double)ts->ptime, (double)dt, (double)new_dt)); 1117 PetscCall(PetscViewerASCIISubtractTab(adapt->monitor, ((PetscObject)adapt)->tablevel)); 1118 } 1119 } 1120 PetscFunctionReturn(PETSC_SUCCESS); 1121 } 1122 1123 /*@ 1124 TSAdaptCreate - create an adaptive controller context for time stepping 1125 1126 Collective 1127 1128 Input Parameter: 1129 . comm - The communicator 1130 1131 Output Parameter: 1132 . inadapt - new `TSAdapt` object 1133 1134 Level: developer 1135 1136 Note: 1137 `TSAdapt` creation is handled by `TS`, so users should not need to call this function. 1138 1139 .seealso: [](ch_ts), [](sec_ts_error_control), `TSAdapt`, `TSGetAdapt()`, `TSAdaptSetType()`, `TSAdaptDestroy()` 1140 @*/ 1141 PetscErrorCode TSAdaptCreate(MPI_Comm comm, TSAdapt *inadapt) 1142 { 1143 TSAdapt adapt; 1144 1145 PetscFunctionBegin; 1146 PetscAssertPointer(inadapt, 2); 1147 PetscCall(TSAdaptInitializePackage()); 1148 1149 PetscCall(PetscHeaderCreate(adapt, TSADAPT_CLASSID, "TSAdapt", "Time stepping adaptivity", "TS", comm, TSAdaptDestroy, TSAdaptView)); 1150 adapt->always_accept = PETSC_FALSE; 1151 adapt->safety = 0.9; 1152 adapt->reject_safety = 0.5; 1153 adapt->clip[0] = 0.1; 1154 adapt->clip[1] = 10.; 1155 adapt->dt_min = 1e-20; 1156 adapt->dt_max = 1e+20; 1157 adapt->ignore_max = -1.0; 1158 adapt->glee_use_local = PETSC_TRUE; 1159 adapt->scale_solve_failed = 0.25; 1160 /* these two safety factors are not public, and they are used only in the TS_EXACTFINALTIME_MATCHSTEP case 1161 to prevent from situations were unreasonably small time steps are taken in order to match the final time */ 1162 adapt->matchstepfac[0] = 0.01; /* allow 1% step size increase in the last step */ 1163 adapt->matchstepfac[1] = 2.0; /* halve last step if it is greater than what remains divided this factor */ 1164 adapt->wnormtype = NORM_2; 1165 adapt->timestepjustdecreased_delay = 0; 1166 *inadapt = adapt; 1167 PetscFunctionReturn(PETSC_SUCCESS); 1168 } 1169