Actual source code: posindep.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1: /*
  2:        Code for Timestepping with implicit backwards Euler.
  3: */
  4:  #include <petsc/private/tsimpl.h>

  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: /* ------------------------------------------------------------------------------*/

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

 34:     Collective on TS

 36:     Input Parameter:
 37: .   ts - timestep context

 39:     Output Parameter:
 40: .   dt - newly computed timestep

 42:     Level: developer

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

 48: .keywords: timestep, pseudo, compute

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

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


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

 69:    Collective on TS

 71:    Input Parameters:
 72: +  ts - the timestep context
 73: .  dtctx - unused timestep context
 74: -  update - latest solution vector

 76:    Output Parameters:
 77: +  newdt - the timestep to use for the next step
 78: -  flag - flag indicating whether the last time step was acceptable

 80:    Level: advanced

 82:    Note:
 83:    This routine always returns a flag of 1, indicating an acceptable
 84:    timestep.

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

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


 98: /*@
 99:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

101:     Collective on TS

103:     Input Parameters:
104: +   ts - timestep context
105: -   update - latest solution vector

107:     Output Parameters:
108: +   dt - newly computed timestep (if it had to shrink)
109: -   flag - indicates if current timestep was ok

111:     Level: advanced

113:     Notes:
114:     The routine to be called here to compute the timestep should be
115:     set by calling TSPseudoSetVerifyTimeStep().

117: .keywords: timestep, pseudo, verify

119: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
120: @*/
121: PetscErrorCode  TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
122: {
123:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

127:   *flag = PETSC_TRUE;
128:   if(pseudo->verify) {
129:     (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
130:   }
131:   return(0);
132: }

134: /* --------------------------------------------------------------------------------*/

136: static PetscErrorCode TSStep_Pseudo(TS ts)
137: {
138:   TS_Pseudo           *pseudo = (TS_Pseudo*)ts->data;
139:   PetscInt            nits,lits,reject;
140:   PetscBool           stepok;
141:   PetscReal           next_time_step = ts->time_step;
142:   PetscErrorCode      ierr;

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

168:   VecCopy(pseudo->update,ts->vec_sol);
169:   ts->ptime += ts->time_step;
170:   ts->time_step = next_time_step;

172:   if (pseudo->fnorm < 0) {
173:     VecZeroEntries(pseudo->xdot);
174:     TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
175:     VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
176:   }
177:   if (pseudo->fnorm < pseudo->fatol) {
178:     ts->reason = TS_CONVERGED_PSEUDO_FATOL;
179:     PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
180:     return(0);
181:   }
182:   if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
183:     ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
184:     PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
185:     return(0);
186:   }
187:   return(0);
188: }

190: /*------------------------------------------------------------*/
191: static PetscErrorCode TSReset_Pseudo(TS ts)
192: {
193:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

197:   VecDestroy(&pseudo->update);
198:   VecDestroy(&pseudo->func);
199:   VecDestroy(&pseudo->xdot);
200:   return(0);
201: }

203: static PetscErrorCode TSDestroy_Pseudo(TS ts)
204: {

208:   TSReset_Pseudo(ts);
209:   PetscFree(ts->data);
210:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
211:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
212:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
213:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
214:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
215:   return(0);
216: }

218: /*------------------------------------------------------------*/

220: /*
221:     Compute Xdot = (X^{n+1}-X^n)/dt) = 0
222: */
223: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
224: {
225:   TS_Pseudo         *pseudo = (TS_Pseudo*)ts->data;
226:   const PetscScalar mdt     = 1.0/ts->time_step,*xnp1,*xn;
227:   PetscScalar       *xdot;
228:   PetscErrorCode    ierr;
229:   PetscInt          i,n;

232:   *Xdot = NULL;
233:   VecGetArrayRead(ts->vec_sol,&xn);
234:   VecGetArrayRead(X,&xnp1);
235:   VecGetArray(pseudo->xdot,&xdot);
236:   VecGetLocalSize(X,&n);
237:   for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
238:   VecRestoreArrayRead(ts->vec_sol,&xn);
239:   VecRestoreArrayRead(X,&xnp1);
240:   VecRestoreArray(pseudo->xdot,&xdot);
241:   *Xdot = pseudo->xdot;
242:   return(0);
243: }

245: /*
246:     The transient residual is

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

250:     or for ODE,

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

254:     This is the function that must be evaluated for transient simulation and for
255:     finite difference Jacobians.  On the first Newton step, this algorithm uses
256:     a guess of U^{n+1} = U^n in which case the transient term vanishes and the
257:     residual is actually the steady state residual.  Pseudotransient
258:     continuation as described in the literature is a linearly implicit
259:     algorithm, it only takes this one Newton step with the steady state
260:     residual, and then advances to the next time step.
261: */
262: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
263: {
264:   Vec            Xdot;

268:   TSPseudoGetXdot(ts,X,&Xdot);
269:   TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
270:   return(0);
271: }

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

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

278:     and for ODE:

280:        J = I/dt - J_{Frhs}   where J_{Frhs} is the given Jacobian of Frhs.
281: */
282: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
283: {
284:   Vec            Xdot;

288:   TSPseudoGetXdot(ts,X,&Xdot);
289:   TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
290:   return(0);
291: }


294: static PetscErrorCode TSSetUp_Pseudo(TS ts)
295: {
296:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

300:   VecDuplicate(ts->vec_sol,&pseudo->update);
301:   VecDuplicate(ts->vec_sol,&pseudo->func);
302:   VecDuplicate(ts->vec_sol,&pseudo->xdot);
303:   return(0);
304: }
305: /*------------------------------------------------------------*/

307: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
308: {
309:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
311:   PetscViewer    viewer = (PetscViewer) dummy;

314:   if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
315:     VecZeroEntries(pseudo->xdot);
316:     TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
317:     VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
318:   }
319:   PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
320:   PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
321:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
322:   return(0);
323: }

325: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
326: {
327:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
329:   PetscBool      flg = PETSC_FALSE;
330:   PetscViewer    viewer;

333:   PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
334:   PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
335:   if (flg) {
336:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
337:     TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
338:   }
339:   flg  = pseudo->increment_dt_from_initial_dt;
340:   PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
341:   pseudo->increment_dt_from_initial_dt = flg;
342:   PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
343:   PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
344:   PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
345:   PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
346:   PetscOptionsTail();
347:   return(0);
348: }

350: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
351: {
353:   PetscBool      isascii;

356:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
357:   if (isascii) {
358:     TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
359:     PetscViewerASCIIPrintf(viewer,"  frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
360:     PetscViewerASCIIPrintf(viewer,"  fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
361:     PetscViewerASCIIPrintf(viewer,"  dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
362:     PetscViewerASCIIPrintf(viewer,"  dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
363:     PetscViewerASCIIPrintf(viewer,"  dt_max - maximum time %g\n",(double)pseudo->dt_max);
364:   }
365:   return(0);
366: }

368: /* ----------------------------------------------------------------------------- */
369: /*@C
370:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
371:    last timestep.

373:    Logically Collective on TS

375:    Input Parameters:
376: +  ts - timestep context
377: .  dt - user-defined function to verify timestep
378: -  ctx - [optional] user-defined context for private data
379:          for the timestep verification routine (may be NULL)

381:    Level: advanced

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

386: .  update - latest solution vector
387: .  ctx - [optional] timestep context
388: .  newdt - the timestep to use for the next step
389: .  flag - flag indicating whether the last time step was acceptable

391:    Notes:
392:    The routine set here will be called by TSPseudoVerifyTimeStep()
393:    during the timestepping process.

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

397: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
398: @*/
399: PetscErrorCode  TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
400: {

405:   PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
406:   return(0);
407: }

409: /*@
410:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
411:     dt when using the TSPseudoTimeStepDefault() routine.

413:    Logically Collective on TS

415:     Input Parameters:
416: +   ts - the timestep context
417: -   inc - the scaling factor >= 1.0

419:     Options Database Key:
420: .    -ts_pseudo_increment <increment>

422:     Level: advanced

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

426: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
427: @*/
428: PetscErrorCode  TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
429: {

435:   PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
436:   return(0);
437: }

439: /*@
440:     TSPseudoSetMaxTimeStep - Sets the maximum time step
441:     when using the TSPseudoTimeStepDefault() routine.

443:    Logically Collective on TS

445:     Input Parameters:
446: +   ts - the timestep context
447: -   maxdt - the maximum time step, use a non-positive value to deactivate

449:     Options Database Key:
450: .    -ts_pseudo_max_dt <increment>

452:     Level: advanced

454: .keywords: timestep, pseudo, set

456: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
457: @*/
458: PetscErrorCode  TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
459: {

465:   PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
466:   return(0);
467: }

469: /*@
470:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
471:     is computed via the formula
472: $         dt = initial_dt*initial_fnorm/current_fnorm
473:       rather than the default update,
474: $         dt = current_dt*previous_fnorm/current_fnorm.

476:    Logically Collective on TS

478:     Input Parameter:
479: .   ts - the timestep context

481:     Options Database Key:
482: .    -ts_pseudo_increment_dt_from_initial_dt

484:     Level: advanced

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

488: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
489: @*/
490: PetscErrorCode  TSPseudoIncrementDtFromInitialDt(TS ts)
491: {

496:   PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
497:   return(0);
498: }


501: /*@C
502:    TSPseudoSetTimeStep - Sets the user-defined routine to be
503:    called at each pseudo-timestep to update the timestep.

505:    Logically Collective on TS

507:    Input Parameters:
508: +  ts - timestep context
509: .  dt - function to compute timestep
510: -  ctx - [optional] user-defined context for private data
511:          required by the function (may be NULL)

513:    Level: intermediate

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

518: .  newdt - the newly computed timestep
519: .  ctx - [optional] timestep context

521:    Notes:
522:    The routine set here will be called by TSPseudoComputeTimeStep()
523:    during the timestepping process.
524:    If not set then TSPseudoTimeStepDefault() is automatically used

526: .keywords: timestep, pseudo, set

528: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
529: @*/
530: PetscErrorCode  TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
531: {

536:   PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
537:   return(0);
538: }

540: /* ----------------------------------------------------------------------------- */

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

548:   pseudo->verify    = dt;
549:   pseudo->verifyctx = ctx;
550:   return(0);
551: }

553: static PetscErrorCode  TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
554: {
555:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

558:   pseudo->dt_increment = inc;
559:   return(0);
560: }

562: static PetscErrorCode  TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
563: {
564:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

567:   pseudo->dt_max = maxdt;
568:   return(0);
569: }

571: static PetscErrorCode  TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
572: {
573:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

576:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
577:   return(0);
578: }

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

586:   pseudo->dt    = dt;
587:   pseudo->dtctx = ctx;
588:   return(0);
589: }

591: /* ----------------------------------------------------------------------------- */
592: /*MC
593:       TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping

595:   This method solves equations of the form

597: $    F(X,Xdot) = 0

599:   for steady state using the iteration

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

604:   where

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

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

611:   Options database keys:
612: +  -ts_pseudo_increment <real> - ratio of increase dt
613: .  -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
614: .  -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
615: -  -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol

617:   Level: beginner

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

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

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

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

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

636: M*/
637: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
638: {
639:   TS_Pseudo      *pseudo;
641:   SNES           snes;
642:   SNESType       stype;

645:   ts->ops->reset          = TSReset_Pseudo;
646:   ts->ops->destroy        = TSDestroy_Pseudo;
647:   ts->ops->view           = TSView_Pseudo;
648:   ts->ops->setup          = TSSetUp_Pseudo;
649:   ts->ops->step           = TSStep_Pseudo;
650:   ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
651:   ts->ops->snesfunction   = SNESTSFormFunction_Pseudo;
652:   ts->ops->snesjacobian   = SNESTSFormJacobian_Pseudo;
653:   ts->default_adapt_type  = TSADAPTNONE;

655:   TSGetSNES(ts,&snes);
656:   SNESGetType(snes,&stype);
657:   if (!stype) {SNESSetType(snes,SNESKSPONLY);}

659:   PetscNewLog(ts,&pseudo);
660:   ts->data = (void*)pseudo;

662:   pseudo->dt                           = TSPseudoTimeStepDefault;
663:   pseudo->dtctx                        = NULL;
664:   pseudo->dt_increment                 = 1.1;
665:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
666:   pseudo->fnorm                        = -1;
667:   pseudo->fnorm_initial                = -1;
668:   pseudo->fnorm_previous               = -1;
669:  #if defined(PETSC_USE_REAL_SINGLE)
670:   pseudo->fatol                        = 1.e-25;
671:   pseudo->frtol                        = 1.e-5;
672: #else
673:   pseudo->fatol                        = 1.e-50;
674:   pseudo->frtol                        = 1.e-12;
675: #endif
676:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
677:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
678:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
679:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
680:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
681:   return(0);
682: }

684: /*@C
685:    TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
686:    Use with TSPseudoSetTimeStep().

688:    Collective on TS

690:    Input Parameters:
691: .  ts - the timestep context
692: .  dtctx - unused timestep context

694:    Output Parameter:
695: .  newdt - the timestep to use for the next step

697:    Level: advanced

699: .keywords: timestep, pseudo, default

701: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
702: @*/
703: PetscErrorCode  TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
704: {
705:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
706:   PetscReal      inc = pseudo->dt_increment;

710:   VecZeroEntries(pseudo->xdot);
711:   TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
712:   VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
713:   if (pseudo->fnorm_initial < 0) {
714:     /* first time through so compute initial function norm */
715:     pseudo->fnorm_initial  = pseudo->fnorm;
716:     pseudo->fnorm_previous = pseudo->fnorm;
717:   }
718:   if (pseudo->fnorm == 0.0)                      *newdt = 1.e12*inc*ts->time_step;
719:   else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
720:   else                                           *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
721:   if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
722:   pseudo->fnorm_previous = pseudo->fnorm;
723:   return(0);
724: }