1 #define PETSCTS_DLL 2 3 /* 4 Code for Timestepping with implicit backwards Euler. 5 */ 6 #include "src/ts/tsimpl.h" /*I "petscts.h" I*/ 7 8 typedef struct { 9 Vec update; /* work vector where new solution is formed */ 10 Vec func; /* work vector where F(t[i],u[i]) is stored */ 11 Vec rhs; /* work vector for RHS; vec_sol/dt */ 12 13 /* information used for Pseudo-timestepping */ 14 15 PetscErrorCode (*dt)(TS,PetscReal*,void*); /* compute next timestep, and related context */ 16 void *dtctx; 17 PetscErrorCode (*verify)(TS,Vec,void*,PetscReal*,PetscTruth*); /* verify previous timestep and related context */ 18 void *verifyctx; 19 20 PetscReal initial_fnorm,fnorm; /* original and current norm of F(u) */ 21 PetscReal fnorm_previous; 22 23 PetscReal dt_increment; /* scaling that dt is incremented each time-step */ 24 PetscTruth increment_dt_from_initial_dt; 25 } TS_Pseudo; 26 27 /* ------------------------------------------------------------------------------*/ 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "TSPseudoComputeTimeStep" 31 /*@ 32 TSPseudoComputeTimeStep - Computes the next timestep for a currently running 33 pseudo-timestepping process. 34 35 Collective on TS 36 37 Input Parameter: 38 . ts - timestep context 39 40 Output Parameter: 41 . dt - newly computed timestep 42 43 Level: advanced 44 45 Notes: 46 The routine to be called here to compute the timestep should be 47 set by calling TSPseudoSetTimeStep(). 48 49 .keywords: timestep, pseudo, compute 50 51 .seealso: TSPseudoDefaultTimeStep(), TSPseudoSetTimeStep() 52 @*/ 53 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoComputeTimeStep(TS ts,PetscReal *dt) 54 { 55 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);CHKERRQ(ierr); 60 ierr = (*pseudo->dt)(ts,dt,pseudo->dtctx);CHKERRQ(ierr); 61 ierr = PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 66 /* ------------------------------------------------------------------------------*/ 67 #undef __FUNCT__ 68 #define __FUNCT__ "TSPseudoDefaultVerifyTimeStep" 69 /*@C 70 TSPseudoDefaultVerifyTimeStep - Default code to verify the quality of the last timestep. 71 72 Collective on TS 73 74 Input Parameters: 75 + ts - the timestep context 76 . dtctx - unused timestep context 77 - update - latest solution vector 78 79 Output Parameters: 80 + newdt - the timestep to use for the next step 81 - flag - flag indicating whether the last time step was acceptable 82 83 Level: advanced 84 85 Note: 86 This routine always returns a flag of 1, indicating an acceptable 87 timestep. 88 89 .keywords: timestep, pseudo, default, verify 90 91 .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep() 92 @*/ 93 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoDefaultVerifyTimeStep(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscTruth *flag) 94 { 95 PetscFunctionBegin; 96 *flag = PETSC_TRUE; 97 PetscFunctionReturn(0); 98 } 99 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "TSPseudoVerifyTimeStep" 103 /*@ 104 TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable. 105 106 Collective on TS 107 108 Input Parameters: 109 + ts - timestep context 110 - update - latest solution vector 111 112 Output Parameters: 113 + dt - newly computed timestep (if it had to shrink) 114 - flag - indicates if current timestep was ok 115 116 Level: advanced 117 118 Notes: 119 The routine to be called here to compute the timestep should be 120 set by calling TSPseudoSetVerifyTimeStep(). 121 122 .keywords: timestep, pseudo, verify 123 124 .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoDefaultVerifyTimeStep() 125 @*/ 126 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscTruth *flag) 127 { 128 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 if (!pseudo->verify) {*flag = PETSC_TRUE; PetscFunctionReturn(0);} 133 134 ierr = (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);CHKERRQ(ierr); 135 136 PetscFunctionReturn(0); 137 } 138 139 /* --------------------------------------------------------------------------------*/ 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "TSStep_Pseudo" 143 static PetscErrorCode TSStep_Pseudo(TS ts,PetscInt *steps,PetscReal *ptime) 144 { 145 Vec sol = ts->vec_sol; 146 PetscErrorCode ierr; 147 PetscInt i,max_steps = ts->max_steps,its,lits; 148 PetscTruth ok; 149 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 150 PetscReal current_time_step; 151 152 PetscFunctionBegin; 153 *steps = -ts->steps; 154 155 ierr = VecCopy(sol,pseudo->update);CHKERRQ(ierr); 156 for (i=0; i<max_steps && ts->ptime < ts->max_time; i++) { 157 ierr = TSPseudoComputeTimeStep(ts,&ts->time_step);CHKERRQ(ierr); 158 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 159 current_time_step = ts->time_step; 160 while (PETSC_TRUE) { 161 ts->ptime += current_time_step; 162 ierr = SNESSolve(ts->snes,PETSC_NULL,pseudo->update);CHKERRQ(ierr); 163 ierr = SNESGetNumberLinearIterations(ts->snes,&lits);CHKERRQ(ierr); 164 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 165 ts->nonlinear_its += its; ts->linear_its += lits; 166 ierr = TSPseudoVerifyTimeStep(ts,pseudo->update,&ts->time_step,&ok);CHKERRQ(ierr); 167 if (ok) break; 168 ts->ptime -= current_time_step; 169 current_time_step = ts->time_step; 170 } 171 ierr = VecCopy(pseudo->update,sol);CHKERRQ(ierr); 172 ts->steps++; 173 } 174 ierr = TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,pseudo->func);CHKERRQ(ierr); 175 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);CHKERRQ(ierr); 176 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 177 178 *steps += ts->steps; 179 *ptime = ts->ptime; 180 PetscFunctionReturn(0); 181 } 182 183 /*------------------------------------------------------------*/ 184 #undef __FUNCT__ 185 #define __FUNCT__ "TSDestroy_Pseudo" 186 static PetscErrorCode TSDestroy_Pseudo(TS ts) 187 { 188 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 if (pseudo->update) {ierr = VecDestroy(pseudo->update);CHKERRQ(ierr);} 193 if (pseudo->func) {ierr = VecDestroy(pseudo->func);CHKERRQ(ierr);} 194 if (pseudo->rhs) {ierr = VecDestroy(pseudo->rhs);CHKERRQ(ierr);} 195 ierr = PetscFree(pseudo);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 200 /*------------------------------------------------------------*/ 201 202 /* 203 This defines the nonlinear equation that is to be solved with SNES 204 205 (U^{n+1} - U^{n})/dt - F(U^{n+1}) 206 */ 207 #undef __FUNCT__ 208 #define __FUNCT__ "TSPseudoFunction" 209 PetscErrorCode TSPseudoFunction(SNES snes,Vec x,Vec y,void *ctx) 210 { 211 TS ts = (TS) ctx; 212 PetscScalar mdt = 1.0/ts->time_step,*unp1,*un,*Funp1; 213 PetscErrorCode ierr; 214 PetscInt i,n; 215 216 PetscFunctionBegin; 217 /* apply user provided function */ 218 ierr = TSComputeRHSFunction(ts,ts->ptime,x,y);CHKERRQ(ierr); 219 /* compute (u^{n+1) - u^{n})/dt - F(u^{n+1}) */ 220 ierr = VecGetArray(ts->vec_sol,&un);CHKERRQ(ierr); 221 ierr = VecGetArray(x,&unp1);CHKERRQ(ierr); 222 ierr = VecGetArray(y,&Funp1);CHKERRQ(ierr); 223 ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr); 224 for (i=0; i<n; i++) { 225 Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i]; 226 } 227 ierr = VecRestoreArray(ts->vec_sol,&un);CHKERRQ(ierr); 228 ierr = VecRestoreArray(x,&unp1);CHKERRQ(ierr); 229 ierr = VecRestoreArray(y,&Funp1);CHKERRQ(ierr); 230 231 PetscFunctionReturn(0); 232 } 233 234 /* 235 This constructs the Jacobian needed for SNES 236 237 J = I/dt - J_{F} where J_{F} is the given Jacobian of F. 238 */ 239 #undef __FUNCT__ 240 #define __FUNCT__ "TSPseudoJacobian" 241 PetscErrorCode TSPseudoJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx) 242 { 243 TS ts = (TS) ctx; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 /* construct users Jacobian */ 248 ierr = TSComputeRHSJacobian(ts,ts->ptime,x,AA,BB,str);CHKERRQ(ierr); 249 250 /* shift and scale Jacobian */ 251 ierr = TSScaleShiftMatrices(ts,*AA,*BB,*str);CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "TSSetUp_Pseudo" 258 static PetscErrorCode TSSetUp_Pseudo(TS ts) 259 { 260 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 /* ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); */ 265 ierr = VecDuplicate(ts->vec_sol,&pseudo->update);CHKERRQ(ierr); 266 ierr = VecDuplicate(ts->vec_sol,&pseudo->func);CHKERRQ(ierr); 267 ierr = SNESSetFunction(ts->snes,pseudo->func,TSPseudoFunction,ts);CHKERRQ(ierr); 268 ierr = SNESSetJacobian(ts->snes,ts->Arhs,ts->B,TSPseudoJacobian,ts);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 /*------------------------------------------------------------*/ 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "TSPseudoMonitorDefault" 275 PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx) 276 { 277 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 ierr = (*PetscHelpPrintf)(ts->comm,"TS %D dt %G time %G fnorm %G\n",step,ts->time_step,ptime,pseudo->fnorm);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "TSSetFromOptions_Pseudo" 287 static PetscErrorCode TSSetFromOptions_Pseudo(TS ts) 288 { 289 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 290 PetscErrorCode ierr; 291 PetscTruth flg; 292 293 PetscFunctionBegin; 294 295 ierr = PetscOptionsHead("Pseudo-timestepping options");CHKERRQ(ierr); 296 ierr = PetscOptionsName("-ts_monitor","Monitor convergence","TSPseudoMonitorDefault",&flg);CHKERRQ(ierr); 297 if (flg) { 298 ierr = TSMonitorSet(ts,TSPseudoMonitorDefault,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 299 } 300 ierr = PetscOptionsName("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",&flg);CHKERRQ(ierr); 301 if (flg) { 302 ierr = TSPseudoIncrementDtFromInitialDt(ts);CHKERRQ(ierr); 303 } 304 ierr = PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,0);CHKERRQ(ierr); 305 ierr = PetscOptionsTail();CHKERRQ(ierr); 306 307 PetscFunctionReturn(0); 308 } 309 310 #undef __FUNCT__ 311 #define __FUNCT__ "TSView_Pseudo" 312 static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer) 313 { 314 PetscFunctionBegin; 315 PetscFunctionReturn(0); 316 } 317 318 /* ----------------------------------------------------------------------------- */ 319 #undef __FUNCT__ 320 #define __FUNCT__ "TSPseudoSetVerifyTimeStep" 321 /*@C 322 TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the 323 last timestep. 324 325 Collective on TS 326 327 Input Parameters: 328 + ts - timestep context 329 . dt - user-defined function to verify timestep 330 - ctx - [optional] user-defined context for private data 331 for the timestep verification routine (may be PETSC_NULL) 332 333 Level: advanced 334 335 Calling sequence of func: 336 . func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscTruth *flag); 337 338 . update - latest solution vector 339 . ctx - [optional] timestep context 340 . newdt - the timestep to use for the next step 341 . flag - flag indicating whether the last time step was acceptable 342 343 Notes: 344 The routine set here will be called by TSPseudoVerifyTimeStep() 345 during the timestepping process. 346 347 .keywords: timestep, pseudo, set, verify 348 349 .seealso: TSPseudoDefaultVerifyTimeStep(), TSPseudoVerifyTimeStep() 350 @*/ 351 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscTruth*),void* ctx) 352 { 353 PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal *,PetscTruth *),void *); 354 355 PetscFunctionBegin; 356 PetscValidHeaderSpecific(ts,TS_COOKIE,1); 357 358 ierr = PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",(void (**)(void))&f);CHKERRQ(ierr); 359 if (f) { 360 ierr = (*f)(ts,dt,ctx);CHKERRQ(ierr); 361 } 362 PetscFunctionReturn(0); 363 } 364 365 #undef __FUNCT__ 366 #define __FUNCT__ "TSPseudoSetTimeStepIncrement" 367 /*@ 368 TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to 369 dt when using the TSPseudoDefaultTimeStep() routine. 370 371 Collective on TS 372 373 Input Parameters: 374 + ts - the timestep context 375 - inc - the scaling factor >= 1.0 376 377 Options Database Key: 378 $ -ts_pseudo_increment <increment> 379 380 Level: advanced 381 382 .keywords: timestep, pseudo, set, increment 383 384 .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep() 385 @*/ 386 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc) 387 { 388 PetscErrorCode ierr,(*f)(TS,PetscReal); 389 390 PetscFunctionBegin; 391 PetscValidHeaderSpecific(ts,TS_COOKIE,1); 392 393 ierr = PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",(void (**)(void))&f);CHKERRQ(ierr); 394 if (f) { 395 ierr = (*f)(ts,inc);CHKERRQ(ierr); 396 } 397 PetscFunctionReturn(0); 398 } 399 400 #undef __FUNCT__ 401 #define __FUNCT__ "TSPseudoIncrementDtFromInitialDt" 402 /*@ 403 TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep 404 is computed via the formula 405 $ dt = initial_dt*initial_fnorm/current_fnorm 406 rather than the default update, 407 $ dt = current_dt*previous_fnorm/current_fnorm. 408 409 Collective on TS 410 411 Input Parameter: 412 . ts - the timestep context 413 414 Options Database Key: 415 $ -ts_pseudo_increment_dt_from_initial_dt 416 417 Level: advanced 418 419 .keywords: timestep, pseudo, set, increment 420 421 .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep() 422 @*/ 423 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoIncrementDtFromInitialDt(TS ts) 424 { 425 PetscErrorCode ierr,(*f)(TS); 426 427 PetscFunctionBegin; 428 PetscValidHeaderSpecific(ts,TS_COOKIE,1); 429 430 ierr = PetscObjectQueryFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",(void (**)(void))&f);CHKERRQ(ierr); 431 if (f) { 432 ierr = (*f)(ts);CHKERRQ(ierr); 433 } 434 PetscFunctionReturn(0); 435 } 436 437 438 #undef __FUNCT__ 439 #define __FUNCT__ "TSPseudoSetTimeStep" 440 /*@C 441 TSPseudoSetTimeStep - Sets the user-defined routine to be 442 called at each pseudo-timestep to update the timestep. 443 444 Collective on TS 445 446 Input Parameters: 447 + ts - timestep context 448 . dt - function to compute timestep 449 - ctx - [optional] user-defined context for private data 450 required by the function (may be PETSC_NULL) 451 452 Level: intermediate 453 454 Calling sequence of func: 455 . func (TS ts,PetscReal *newdt,void *ctx); 456 457 . newdt - the newly computed timestep 458 . ctx - [optional] timestep context 459 460 Notes: 461 The routine set here will be called by TSPseudoComputeTimeStep() 462 during the timestepping process. 463 464 .keywords: timestep, pseudo, set 465 466 .seealso: TSPseudoDefaultTimeStep(), TSPseudoComputeTimeStep() 467 @*/ 468 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void* ctx) 469 { 470 PetscErrorCode ierr,(*f)(TS,PetscErrorCode (*)(TS,PetscReal *,void *),void *); 471 472 PetscFunctionBegin; 473 PetscValidHeaderSpecific(ts,TS_COOKIE,1); 474 475 ierr = PetscObjectQueryFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",(void (**)(void))&f);CHKERRQ(ierr); 476 if (f) { 477 ierr = (*f)(ts,dt,ctx);CHKERRQ(ierr); 478 } 479 PetscFunctionReturn(0); 480 } 481 482 /* ----------------------------------------------------------------------------- */ 483 484 typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscTruth*); /* force argument to next function to not be extern C*/ 485 EXTERN_C_BEGIN 486 #undef __FUNCT__ 487 #define __FUNCT__ "TSPseudoSetVerifyTimeStep_Pseudo" 488 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void* ctx) 489 { 490 TS_Pseudo *pseudo; 491 492 PetscFunctionBegin; 493 pseudo = (TS_Pseudo*)ts->data; 494 pseudo->verify = dt; 495 pseudo->verifyctx = ctx; 496 PetscFunctionReturn(0); 497 } 498 EXTERN_C_END 499 500 EXTERN_C_BEGIN 501 #undef __FUNCT__ 502 #define __FUNCT__ "TSPseudoSetTimeStepIncrement_Pseudo" 503 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc) 504 { 505 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 506 507 PetscFunctionBegin; 508 pseudo->dt_increment = inc; 509 PetscFunctionReturn(0); 510 } 511 EXTERN_C_END 512 513 EXTERN_C_BEGIN 514 #undef __FUNCT__ 515 #define __FUNCT__ "TSPseudoIncrementDtFromInitialDt_Pseudo" 516 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts) 517 { 518 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 519 520 PetscFunctionBegin; 521 pseudo->increment_dt_from_initial_dt = PETSC_TRUE; 522 PetscFunctionReturn(0); 523 } 524 EXTERN_C_END 525 526 typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/ 527 EXTERN_C_BEGIN 528 #undef __FUNCT__ 529 #define __FUNCT__ "TSPseudoSetTimeStep_Pseudo" 530 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void* ctx) 531 { 532 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 533 534 PetscFunctionBegin; 535 pseudo->dt = dt; 536 pseudo->dtctx = ctx; 537 PetscFunctionReturn(0); 538 } 539 EXTERN_C_END 540 541 /* ----------------------------------------------------------------------------- */ 542 543 EXTERN_C_BEGIN 544 #undef __FUNCT__ 545 #define __FUNCT__ "TSCreate_Pseudo" 546 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Pseudo(TS ts) 547 { 548 TS_Pseudo *pseudo; 549 PetscErrorCode ierr; 550 551 PetscFunctionBegin; 552 ts->ops->destroy = TSDestroy_Pseudo; 553 ts->ops->view = TSView_Pseudo; 554 555 if (ts->problem_type == TS_LINEAR) { 556 SETERRQ(PETSC_ERR_ARG_WRONG,"Only for nonlinear problems"); 557 } 558 if (!ts->Arhs) { 559 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set Jacobian"); 560 } 561 562 ts->ops->setup = TSSetUp_Pseudo; 563 ts->ops->step = TSStep_Pseudo; 564 ts->ops->setfromoptions = TSSetFromOptions_Pseudo; 565 566 /* create the required nonlinear solver context */ 567 ierr = SNESCreate(ts->comm,&ts->snes);CHKERRQ(ierr); 568 569 ierr = PetscNew(TS_Pseudo,&pseudo);CHKERRQ(ierr); 570 ierr = PetscLogObjectMemory(ts,sizeof(TS_Pseudo));CHKERRQ(ierr); 571 572 ierr = PetscMemzero(pseudo,sizeof(TS_Pseudo));CHKERRQ(ierr); 573 ts->data = (void*)pseudo; 574 575 pseudo->dt_increment = 1.1; 576 pseudo->increment_dt_from_initial_dt = PETSC_FALSE; 577 pseudo->dt = TSPseudoDefaultTimeStep; 578 579 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C", 580 "TSPseudoSetVerifyTimeStep_Pseudo", 581 TSPseudoSetVerifyTimeStep_Pseudo);CHKERRQ(ierr); 582 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C", 583 "TSPseudoSetTimeStepIncrement_Pseudo", 584 TSPseudoSetTimeStepIncrement_Pseudo);CHKERRQ(ierr); 585 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C", 586 "TSPseudoIncrementDtFromInitialDt_Pseudo", 587 TSPseudoIncrementDtFromInitialDt_Pseudo);CHKERRQ(ierr); 588 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStep_C", 589 "TSPseudoSetTimeStep_Pseudo", 590 TSPseudoSetTimeStep_Pseudo);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 EXTERN_C_END 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "TSPseudoDefaultTimeStep" 597 /*@C 598 TSPseudoDefaultTimeStep - Default code to compute pseudo-timestepping. 599 Use with TSPseudoSetTimeStep(). 600 601 Collective on TS 602 603 Input Parameters: 604 . ts - the timestep context 605 . dtctx - unused timestep context 606 607 Output Parameter: 608 . newdt - the timestep to use for the next step 609 610 Level: advanced 611 612 .keywords: timestep, pseudo, default 613 614 .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep() 615 @*/ 616 PetscErrorCode PETSCTS_DLLEXPORT TSPseudoDefaultTimeStep(TS ts,PetscReal* newdt,void* dtctx) 617 { 618 TS_Pseudo *pseudo = (TS_Pseudo*)ts->data; 619 PetscReal inc = pseudo->dt_increment,fnorm_previous = pseudo->fnorm_previous; 620 PetscErrorCode ierr; 621 622 PetscFunctionBegin; 623 ierr = TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,pseudo->func);CHKERRQ(ierr); 624 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);CHKERRQ(ierr); 625 if (pseudo->initial_fnorm == 0.0) { 626 /* first time through so compute initial function norm */ 627 pseudo->initial_fnorm = pseudo->fnorm; 628 fnorm_previous = pseudo->fnorm; 629 } 630 if (pseudo->fnorm == 0.0) { 631 *newdt = 1.e12*inc*ts->time_step; 632 } else if (pseudo->increment_dt_from_initial_dt) { 633 *newdt = inc*ts->initial_time_step*pseudo->initial_fnorm/pseudo->fnorm; 634 } else { 635 *newdt = inc*ts->time_step*fnorm_previous/pseudo->fnorm; 636 } 637 pseudo->fnorm_previous = pseudo->fnorm; 638 PetscFunctionReturn(0); 639 } 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657