Actual source code: posindep.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: /*
  2:        Code for Timestepping with implicit backwards Euler.
  3: */
  4: #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/

  6: typedef struct {
  7:   Vec update;       /* work vector where new solution is formed */
  8:   Vec func;         /* work vector where F(t[i],u[i]) is stored */
  9:   Vec xdot;         /* work vector for time derivative of state */

 11:   /* information used for Pseudo-timestepping */

 13:   PetscErrorCode (*dt)(TS,PetscReal*,void*);              /* compute next timestep, and related context */
 14:   void *dtctx;
 15:   PetscErrorCode (*verify)(TS,Vec,void*,PetscReal*,PetscBool*);  /* verify previous timestep and related context */
 16:   void *verifyctx;

 18:   PetscReal fnorm_initial,fnorm;                   /* original and current norm of F(u) */
 19:   PetscReal fnorm_previous;

 21:   PetscReal dt_initial;                     /* initial time-step */
 22:   PetscReal dt_increment;                   /* scaling that dt is incremented each time-step */
 23:   PetscReal dt_max;                         /* maximum time step */
 24:   PetscBool increment_dt_from_initial_dt;
 25:   PetscReal fatol,frtol;
 26: } TS_Pseudo;

 28: /* ------------------------------------------------------------------------------*/

 32: /*@C
 33:     TSPseudoComputeTimeStep - Computes the next timestep for a currently running
 34:     pseudo-timestepping process.

 36:     Collective on TS

 38:     Input Parameter:
 39: .   ts - timestep context

 41:     Output Parameter:
 42: .   dt - newly computed timestep

 44:     Level: developer

 46:     Notes:
 47:     The routine to be called here to compute the timestep should be
 48:     set by calling TSPseudoSetTimeStep().

 50: .keywords: timestep, pseudo, compute

 52: .seealso: TSPseudoTimeStepDefault(), TSPseudoSetTimeStep()
 53: @*/
 54: PetscErrorCode  TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
 55: {
 56:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

 60:   PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);
 61:   (*pseudo->dt)(ts,dt,pseudo->dtctx);
 62:   PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);
 63:   return(0);
 64: }


 67: /* ------------------------------------------------------------------------------*/
 70: /*@C
 71:    TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.

 73:    Collective on TS

 75:    Input Parameters:
 76: +  ts - the timestep context
 77: .  dtctx - unused timestep context
 78: -  update - latest solution vector

 80:    Output Parameters:
 81: +  newdt - the timestep to use for the next step
 82: -  flag - flag indicating whether the last time step was acceptable

 84:    Level: advanced

 86:    Note:
 87:    This routine always returns a flag of 1, indicating an acceptable
 88:    timestep.

 90: .keywords: timestep, pseudo, default, verify

 92: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
 93: @*/
 94: PetscErrorCode  TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool  *flag)
 95: {
 97:   *flag = PETSC_TRUE;
 98:   return(0);
 99: }


104: /*@
105:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

107:     Collective on TS

109:     Input Parameters:
110: +   ts - timestep context
111: -   update - latest solution vector

113:     Output Parameters:
114: +   dt - newly computed timestep (if it had to shrink)
115: -   flag - indicates if current timestep was ok

117:     Level: advanced

119:     Notes:
120:     The routine to be called here to compute the timestep should be
121:     set by calling TSPseudoSetVerifyTimeStep().

123: .keywords: timestep, pseudo, verify

125: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
126: @*/
127: PetscErrorCode  TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
128: {
129:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

133:   *flag = PETSC_TRUE;
134:   if(pseudo->verify) {
135:     (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
136:   }
137:   return(0);
138: }

140: /* --------------------------------------------------------------------------------*/

144: static PetscErrorCode TSStep_Pseudo(TS ts)
145: {
146:   TS_Pseudo           *pseudo = (TS_Pseudo*)ts->data;
147:   PetscInt            nits,lits,reject;
148:   PetscBool           stepok;
149:   PetscReal           next_time_step = ts->time_step;
150:   PetscErrorCode      ierr;

153:   if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
154:   VecCopy(ts->vec_sol,pseudo->update);
155:   TSPseudoComputeTimeStep(ts,&next_time_step);
156:   for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
157:     ts->time_step = next_time_step;
158:     TSPreStage(ts,ts->ptime+ts->time_step);
159:     SNESSolve(ts->snes,NULL,pseudo->update);
160:     SNESGetIterationNumber(ts->snes,&nits);
161:     SNESGetLinearSolveIterations(ts->snes,&lits);
162:     ts->snes_its += nits; ts->ksp_its += lits;
163:     TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));
164:     TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);
165:     if (!stepok) {next_time_step = ts->time_step; continue;}
166:     pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
167:     TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
168:     if (stepok) break;
169:   }
170:   if (reject >= ts->max_reject) {
171:     ts->reason = TS_DIVERGED_STEP_REJECTED;
172:     PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
173:     return(0);
174:   }

176:   VecCopy(pseudo->update,ts->vec_sol);
177:   ts->ptime += ts->time_step;
178:   ts->time_step = next_time_step;

180:   if (pseudo->fnorm < 0) {
181:     VecZeroEntries(pseudo->xdot);
182:     TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
183:     VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
184:   }
185:   if (pseudo->fnorm < pseudo->fatol) {
186:     ts->reason = TS_CONVERGED_PSEUDO_FATOL;
187:     PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
188:     return(0);
189:   }
190:   if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
191:     ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
192:     PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
193:     return(0);
194:   }
195:   return(0);
196: }

198: /*------------------------------------------------------------*/
201: static PetscErrorCode TSReset_Pseudo(TS ts)
202: {
203:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

207:   VecDestroy(&pseudo->update);
208:   VecDestroy(&pseudo->func);
209:   VecDestroy(&pseudo->xdot);
210:   return(0);
211: }

215: static PetscErrorCode TSDestroy_Pseudo(TS ts)
216: {

220:   TSReset_Pseudo(ts);
221:   PetscFree(ts->data);
222:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
223:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
224:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
225:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
226:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
227:   return(0);
228: }

230: /*------------------------------------------------------------*/

234: /*
235:     Compute Xdot = (X^{n+1}-X^n)/dt) = 0
236: */
237: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
238: {
239:   TS_Pseudo         *pseudo = (TS_Pseudo*)ts->data;
240:   const PetscScalar mdt     = 1.0/ts->time_step,*xnp1,*xn;
241:   PetscScalar       *xdot;
242:   PetscErrorCode    ierr;
243:   PetscInt          i,n;

246:   VecGetArrayRead(ts->vec_sol,&xn);
247:   VecGetArrayRead(X,&xnp1);
248:   VecGetArray(pseudo->xdot,&xdot);
249:   VecGetLocalSize(X,&n);
250:   for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
251:   VecRestoreArrayRead(ts->vec_sol,&xn);
252:   VecRestoreArrayRead(X,&xnp1);
253:   VecRestoreArray(pseudo->xdot,&xdot);
254:   *Xdot = pseudo->xdot;
255:   return(0);
256: }

260: /*
261:     The transient residual is

263:         F(U^{n+1},(U^{n+1}-U^n)/dt) = 0

265:     or for ODE,

267:         (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0

269:     This is the function that must be evaluated for transient simulation and for
270:     finite difference Jacobians.  On the first Newton step, this algorithm uses
271:     a guess of U^{n+1} = U^n in which case the transient term vanishes and the
272:     residual is actually the steady state residual.  Pseudotransient
273:     continuation as described in the literature is a linearly implicit
274:     algorithm, it only takes this one Newton step with the steady state
275:     residual, and then advances to the next time step.
276: */
277: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
278: {
279:   Vec            Xdot;

283:   TSPseudoGetXdot(ts,X,&Xdot);
284:   TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
285:   return(0);
286: }

290: /*
291:    This constructs the Jacobian needed for SNES.  For DAE, this is

293:        dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot

295:     and for ODE:

297:        J = I/dt - J_{Frhs}   where J_{Frhs} is the given Jacobian of Frhs.
298: */
299: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
300: {
301:   Vec            Xdot;

305:   TSPseudoGetXdot(ts,X,&Xdot);
306:   TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
307:   return(0);
308: }


313: static PetscErrorCode TSSetUp_Pseudo(TS ts)
314: {
315:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

319:   VecDuplicate(ts->vec_sol,&pseudo->update);
320:   VecDuplicate(ts->vec_sol,&pseudo->func);
321:   VecDuplicate(ts->vec_sol,&pseudo->xdot);
322:   return(0);
323: }
324: /*------------------------------------------------------------*/

328: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
329: {
330:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
332:   PetscViewer    viewer = (PetscViewer) dummy;

335:   if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
336:     VecZeroEntries(pseudo->xdot);
337:     TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
338:     VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
339:   }
340:   PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
341:   PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
342:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
343:   return(0);
344: }

348: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
349: {
350:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
352:   PetscBool      flg = PETSC_FALSE;
353:   PetscViewer    viewer;

356:   PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
357:   PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
358:   if (flg) {
359:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
360:     TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
361:   }
362:   flg  = pseudo->increment_dt_from_initial_dt;
363:   PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
364:   pseudo->increment_dt_from_initial_dt = flg;
365:   PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
366:   PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
367:   PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
368:   PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
369:   PetscOptionsTail();
370:   return(0);
371: }

375: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
376: {
378:   PetscBool      isascii;

381:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
382:   if (isascii) {
383:     TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
384:     PetscViewerASCIIPrintf(viewer,"Parameters for pseudo timestepping\n");
385:     PetscViewerASCIIPrintf(viewer,"  frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
386:     PetscViewerASCIIPrintf(viewer,"  fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
387:     PetscViewerASCIIPrintf(viewer,"  dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
388:     PetscViewerASCIIPrintf(viewer,"  dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
389:     PetscViewerASCIIPrintf(viewer,"  dt_max - maximum time %g\n",(double)pseudo->dt_max);
390:   }
391:   if (ts->snes) {SNESView(ts->snes,viewer);}
392:   return(0);
393: }

395: /* ----------------------------------------------------------------------------- */
398: /*@C
399:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
400:    last timestep.

402:    Logically Collective on TS

404:    Input Parameters:
405: +  ts - timestep context
406: .  dt - user-defined function to verify timestep
407: -  ctx - [optional] user-defined context for private data
408:          for the timestep verification routine (may be NULL)

410:    Level: advanced

412:    Calling sequence of func:
413: .  func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool  *flag);

415: .  update - latest solution vector
416: .  ctx - [optional] timestep context
417: .  newdt - the timestep to use for the next step
418: .  flag - flag indicating whether the last time step was acceptable

420:    Notes:
421:    The routine set here will be called by TSPseudoVerifyTimeStep()
422:    during the timestepping process.

424: .keywords: timestep, pseudo, set, verify

426: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
427: @*/
428: PetscErrorCode  TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
429: {

434:   PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
435:   return(0);
436: }

440: /*@
441:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
442:     dt when using the TSPseudoTimeStepDefault() routine.

444:    Logically Collective on TS

446:     Input Parameters:
447: +   ts - the timestep context
448: -   inc - the scaling factor >= 1.0

450:     Options Database Key:
451: .    -ts_pseudo_increment <increment>

453:     Level: advanced

455: .keywords: timestep, pseudo, set, increment

457: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
458: @*/
459: PetscErrorCode  TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
460: {

466:   PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
467:   return(0);
468: }

472: /*@
473:     TSPseudoSetMaxTimeStep - Sets the maximum time step
474:     when using the TSPseudoTimeStepDefault() routine.

476:    Logically Collective on TS

478:     Input Parameters:
479: +   ts - the timestep context
480: -   maxdt - the maximum time step, use a non-positive value to deactivate

482:     Options Database Key:
483: .    -ts_pseudo_max_dt <increment>

485:     Level: advanced

487: .keywords: timestep, pseudo, set

489: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
490: @*/
491: PetscErrorCode  TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
492: {

498:   PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
499:   return(0);
500: }

504: /*@
505:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
506:     is computed via the formula
507: $         dt = initial_dt*initial_fnorm/current_fnorm
508:       rather than the default update,
509: $         dt = current_dt*previous_fnorm/current_fnorm.

511:    Logically Collective on TS

513:     Input Parameter:
514: .   ts - the timestep context

516:     Options Database Key:
517: .    -ts_pseudo_increment_dt_from_initial_dt

519:     Level: advanced

521: .keywords: timestep, pseudo, set, increment

523: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
524: @*/
525: PetscErrorCode  TSPseudoIncrementDtFromInitialDt(TS ts)
526: {

531:   PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
532:   return(0);
533: }


538: /*@C
539:    TSPseudoSetTimeStep - Sets the user-defined routine to be
540:    called at each pseudo-timestep to update the timestep.

542:    Logically Collective on TS

544:    Input Parameters:
545: +  ts - timestep context
546: .  dt - function to compute timestep
547: -  ctx - [optional] user-defined context for private data
548:          required by the function (may be NULL)

550:    Level: intermediate

552:    Calling sequence of func:
553: .  func (TS ts,PetscReal *newdt,void *ctx);

555: .  newdt - the newly computed timestep
556: .  ctx - [optional] timestep context

558:    Notes:
559:    The routine set here will be called by TSPseudoComputeTimeStep()
560:    during the timestepping process.
561:    If not set then TSPseudoTimeStepDefault() is automatically used

563: .keywords: timestep, pseudo, set

565: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
566: @*/
567: PetscErrorCode  TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
568: {

573:   PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
574:   return(0);
575: }

577: /* ----------------------------------------------------------------------------- */

579: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*);  /* force argument to next function to not be extern C*/
582: static PetscErrorCode  TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
583: {
584:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

587:   pseudo->verify    = dt;
588:   pseudo->verifyctx = ctx;
589:   return(0);
590: }

594: static PetscErrorCode  TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
595: {
596:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

599:   pseudo->dt_increment = inc;
600:   return(0);
601: }

605: static PetscErrorCode  TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
606: {
607:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

610:   pseudo->dt_max = maxdt;
611:   return(0);
612: }

616: static PetscErrorCode  TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
617: {
618:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

621:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
622:   return(0);
623: }

625: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
628: static PetscErrorCode  TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
629: {
630:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

633:   pseudo->dt    = dt;
634:   pseudo->dtctx = ctx;
635:   return(0);
636: }

638: /* ----------------------------------------------------------------------------- */
639: /*MC
640:       TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping

642:   This method solves equations of the form

644: $    F(X,Xdot) = 0

646:   for steady state using the iteration

648: $    [G'] S = -F(X,0)
649: $    X += S

651:   where

653: $    G(Y) = F(Y,(Y-X)/dt)

655:   This is linearly-implicit Euler with the residual always evaluated "at steady
656:   state".  See note below.

658:   Options database keys:
659: +  -ts_pseudo_increment <real> - ratio of increase dt
660: .  -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
661: .  -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
662: -  -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol

664:   Level: beginner

666:   References:
667: +  1. - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
668: -  2. - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.

670:   Notes:
671:   The residual computed by this method includes the transient term (Xdot is computed instead of
672:   always being zero), but since the prediction from the last step is always the solution from the
673:   last step, on the first Newton iteration we have

675: $  Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0

677:   Therefore, the linear system solved by the first Newton iteration is equivalent to the one
678:   described above and in the papers.  If the user chooses to perform multiple Newton iterations, the
679:   algorithm is no longer the one described in the referenced papers.

681: .seealso:  TSCreate(), TS, TSSetType()

683: M*/
686: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
687: {
688:   TS_Pseudo      *pseudo;
690:   SNES           snes;
691:   SNESType       stype;

694:   ts->ops->reset          = TSReset_Pseudo;
695:   ts->ops->destroy        = TSDestroy_Pseudo;
696:   ts->ops->view           = TSView_Pseudo;
697:   ts->ops->setup          = TSSetUp_Pseudo;
698:   ts->ops->step           = TSStep_Pseudo;
699:   ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
700:   ts->ops->snesfunction   = SNESTSFormFunction_Pseudo;
701:   ts->ops->snesjacobian   = SNESTSFormJacobian_Pseudo;

703:   TSGetSNES(ts,&snes);
704:   SNESGetType(snes,&stype);
705:   if (!stype) {SNESSetType(snes,SNESKSPONLY);}

707:   PetscNewLog(ts,&pseudo);
708:   ts->data = (void*)pseudo;

710:   pseudo->dt                           = TSPseudoTimeStepDefault;
711:   pseudo->dtctx                        = NULL;
712:   pseudo->dt_increment                 = 1.1;
713:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
714:   pseudo->fnorm                        = -1;
715:   pseudo->fnorm_initial                = -1;
716:   pseudo->fnorm_previous               = -1;
717:  #if defined(PETSC_USE_REAL_SINGLE)
718:   pseudo->fatol                        = 1.e-25;
719:   pseudo->frtol                        = 1.e-5;
720: #else
721:   pseudo->fatol                        = 1.e-50;
722:   pseudo->frtol                        = 1.e-12;
723: #endif
724:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
725:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
726:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
727:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
728:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
729:   return(0);
730: }

734: /*@C
735:    TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
736:    Use with TSPseudoSetTimeStep().

738:    Collective on TS

740:    Input Parameters:
741: .  ts - the timestep context
742: .  dtctx - unused timestep context

744:    Output Parameter:
745: .  newdt - the timestep to use for the next step

747:    Level: advanced

749: .keywords: timestep, pseudo, default

751: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
752: @*/
753: PetscErrorCode  TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
754: {
755:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
756:   PetscReal      inc = pseudo->dt_increment;

760:   VecZeroEntries(pseudo->xdot);
761:   TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
762:   VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
763:   if (pseudo->fnorm_initial < 0) {
764:     /* first time through so compute initial function norm */
765:     pseudo->fnorm_initial  = pseudo->fnorm;
766:     pseudo->fnorm_previous = pseudo->fnorm;
767:   }
768:   if (pseudo->fnorm == 0.0)                      *newdt = 1.e12*inc*ts->time_step;
769:   else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
770:   else                                           *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
771:   if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
772:   pseudo->fnorm_previous = pseudo->fnorm;
773:   return(0);
774: }