Actual source code: posindep.c

petsc-3.12.5 2020-03-29
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: .seealso: TSPseudoTimeStepDefault(), TSPseudoSetTimeStep()
 49: @*/
 50: PetscErrorCode  TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
 51: {
 52:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

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


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

 67:    Collective on TS

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

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

 78:    Level: advanced

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

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


 94: /*@
 95:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

 97:     Collective on TS

 99:     Input Parameters:
100: +   ts - timestep context
101: -   update - latest solution vector

103:     Output Parameters:
104: +   dt - newly computed timestep (if it had to shrink)
105: -   flag - indicates if current timestep was ok

107:     Level: advanced

109:     Notes:
110:     The routine to be called here to compute the timestep should be
111:     set by calling TSPseudoSetVerifyTimeStep().

113: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
114: @*/
115: PetscErrorCode  TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
116: {
117:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

121:   *flag = PETSC_TRUE;
122:   if(pseudo->verify) {
123:     (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
124:   }
125:   return(0);
126: }

128: /* --------------------------------------------------------------------------------*/

130: static PetscErrorCode TSStep_Pseudo(TS ts)
131: {
132:   TS_Pseudo           *pseudo = (TS_Pseudo*)ts->data;
133:   PetscInt            nits,lits,reject;
134:   PetscBool           stepok;
135:   PetscReal           next_time_step = ts->time_step;
136:   PetscErrorCode      ierr;

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

162:   VecCopy(pseudo->update,ts->vec_sol);
163:   ts->ptime += ts->time_step;
164:   ts->time_step = next_time_step;

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

184: /*------------------------------------------------------------*/
185: static PetscErrorCode TSReset_Pseudo(TS ts)
186: {
187:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

191:   VecDestroy(&pseudo->update);
192:   VecDestroy(&pseudo->func);
193:   VecDestroy(&pseudo->xdot);
194:   return(0);
195: }

197: static PetscErrorCode TSDestroy_Pseudo(TS ts)
198: {

202:   TSReset_Pseudo(ts);
203:   PetscFree(ts->data);
204:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
205:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
206:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
207:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
208:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
209:   return(0);
210: }

212: /*------------------------------------------------------------*/

214: /*
215:     Compute Xdot = (X^{n+1}-X^n)/dt) = 0
216: */
217: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
218: {
219:   TS_Pseudo         *pseudo = (TS_Pseudo*)ts->data;
220:   const PetscScalar mdt     = 1.0/ts->time_step,*xnp1,*xn;
221:   PetscScalar       *xdot;
222:   PetscErrorCode    ierr;
223:   PetscInt          i,n;

226:   *Xdot = NULL;
227:   VecGetArrayRead(ts->vec_sol,&xn);
228:   VecGetArrayRead(X,&xnp1);
229:   VecGetArray(pseudo->xdot,&xdot);
230:   VecGetLocalSize(X,&n);
231:   for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
232:   VecRestoreArrayRead(ts->vec_sol,&xn);
233:   VecRestoreArrayRead(X,&xnp1);
234:   VecRestoreArray(pseudo->xdot,&xdot);
235:   *Xdot = pseudo->xdot;
236:   return(0);
237: }

239: /*
240:     The transient residual is

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

244:     or for ODE,

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

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

262:   TSPseudoGetXdot(ts,X,&Xdot);
263:   TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
264:   return(0);
265: }

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

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

272:     and for ODE:

274:        J = I/dt - J_{Frhs}   where J_{Frhs} is the given Jacobian of Frhs.
275: */
276: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
277: {
278:   Vec            Xdot;

282:   TSPseudoGetXdot(ts,X,&Xdot);
283:   TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
284:   return(0);
285: }


288: static PetscErrorCode TSSetUp_Pseudo(TS ts)
289: {
290:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;

294:   VecDuplicate(ts->vec_sol,&pseudo->update);
295:   VecDuplicate(ts->vec_sol,&pseudo->func);
296:   VecDuplicate(ts->vec_sol,&pseudo->xdot);
297:   return(0);
298: }
299: /*------------------------------------------------------------*/

301: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
302: {
303:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
305:   PetscViewer    viewer = (PetscViewer) dummy;

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

319: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
320: {
321:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
323:   PetscBool      flg = PETSC_FALSE;
324:   PetscViewer    viewer;

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

344: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
345: {
347:   PetscBool      isascii;

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

362: /* ----------------------------------------------------------------------------- */
363: /*@C
364:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
365:    last timestep.

367:    Logically Collective on TS

369:    Input Parameters:
370: +  ts - timestep context
371: .  dt - user-defined function to verify timestep
372: -  ctx - [optional] user-defined context for private data
373:          for the timestep verification routine (may be NULL)

375:    Level: advanced

377:    Calling sequence of func:
378: $  func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool  *flag);

380: +  update - latest solution vector
381: .  ctx - [optional] timestep context
382: .  newdt - the timestep to use for the next step
383: -  flag - flag indicating whether the last time step was acceptable

385:    Notes:
386:    The routine set here will be called by TSPseudoVerifyTimeStep()
387:    during the timestepping process.

389: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
390: @*/
391: PetscErrorCode  TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
392: {

397:   PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
398:   return(0);
399: }

401: /*@
402:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
403:     dt when using the TSPseudoTimeStepDefault() routine.

405:    Logically Collective on TS

407:     Input Parameters:
408: +   ts - the timestep context
409: -   inc - the scaling factor >= 1.0

411:     Options Database Key:
412: .    -ts_pseudo_increment <increment>

414:     Level: advanced

416: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
417: @*/
418: PetscErrorCode  TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
419: {

425:   PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
426:   return(0);
427: }

429: /*@
430:     TSPseudoSetMaxTimeStep - Sets the maximum time step
431:     when using the TSPseudoTimeStepDefault() routine.

433:    Logically Collective on TS

435:     Input Parameters:
436: +   ts - the timestep context
437: -   maxdt - the maximum time step, use a non-positive value to deactivate

439:     Options Database Key:
440: .    -ts_pseudo_max_dt <increment>

442:     Level: advanced

444: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
445: @*/
446: PetscErrorCode  TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
447: {

453:   PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
454:   return(0);
455: }

457: /*@
458:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
459:     is computed via the formula
460: $         dt = initial_dt*initial_fnorm/current_fnorm
461:       rather than the default update,
462: $         dt = current_dt*previous_fnorm/current_fnorm.

464:    Logically Collective on TS

466:     Input Parameter:
467: .   ts - the timestep context

469:     Options Database Key:
470: .    -ts_pseudo_increment_dt_from_initial_dt

472:     Level: advanced

474: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
475: @*/
476: PetscErrorCode  TSPseudoIncrementDtFromInitialDt(TS ts)
477: {

482:   PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
483:   return(0);
484: }


487: /*@C
488:    TSPseudoSetTimeStep - Sets the user-defined routine to be
489:    called at each pseudo-timestep to update the timestep.

491:    Logically Collective on TS

493:    Input Parameters:
494: +  ts - timestep context
495: .  dt - function to compute timestep
496: -  ctx - [optional] user-defined context for private data
497:          required by the function (may be NULL)

499:    Level: intermediate

501:    Calling sequence of func:
502: $  func (TS ts,PetscReal *newdt,void *ctx);

504: +  newdt - the newly computed timestep
505: -  ctx - [optional] timestep context

507:    Notes:
508:    The routine set here will be called by TSPseudoComputeTimeStep()
509:    during the timestepping process.
510:    If not set then TSPseudoTimeStepDefault() is automatically used

512: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
513: @*/
514: PetscErrorCode  TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
515: {

520:   PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
521:   return(0);
522: }

524: /* ----------------------------------------------------------------------------- */

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

532:   pseudo->verify    = dt;
533:   pseudo->verifyctx = ctx;
534:   return(0);
535: }

537: static PetscErrorCode  TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
538: {
539:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

542:   pseudo->dt_increment = inc;
543:   return(0);
544: }

546: static PetscErrorCode  TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
547: {
548:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

551:   pseudo->dt_max = maxdt;
552:   return(0);
553: }

555: static PetscErrorCode  TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
556: {
557:   TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;

560:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
561:   return(0);
562: }

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

570:   pseudo->dt    = dt;
571:   pseudo->dtctx = ctx;
572:   return(0);
573: }

575: /* ----------------------------------------------------------------------------- */
576: /*MC
577:       TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping

579:   This method solves equations of the form

581: $    F(X,Xdot) = 0

583:   for steady state using the iteration

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

588:   where

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

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

595:   Options database keys:
596: +  -ts_pseudo_increment <real> - ratio of increase dt
597: .  -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
598: .  -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
599: -  -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol

601:   Level: beginner

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

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

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

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

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

620: M*/
621: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
622: {
623:   TS_Pseudo      *pseudo;
625:   SNES           snes;
626:   SNESType       stype;

629:   ts->ops->reset          = TSReset_Pseudo;
630:   ts->ops->destroy        = TSDestroy_Pseudo;
631:   ts->ops->view           = TSView_Pseudo;
632:   ts->ops->setup          = TSSetUp_Pseudo;
633:   ts->ops->step           = TSStep_Pseudo;
634:   ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
635:   ts->ops->snesfunction   = SNESTSFormFunction_Pseudo;
636:   ts->ops->snesjacobian   = SNESTSFormJacobian_Pseudo;
637:   ts->default_adapt_type  = TSADAPTNONE;
638:   ts->usessnes            = PETSC_TRUE;

640:   TSGetSNES(ts,&snes);
641:   SNESGetType(snes,&stype);
642:   if (!stype) {SNESSetType(snes,SNESKSPONLY);}

644:   PetscNewLog(ts,&pseudo);
645:   ts->data = (void*)pseudo;

647:   pseudo->dt                           = TSPseudoTimeStepDefault;
648:   pseudo->dtctx                        = NULL;
649:   pseudo->dt_increment                 = 1.1;
650:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
651:   pseudo->fnorm                        = -1;
652:   pseudo->fnorm_initial                = -1;
653:   pseudo->fnorm_previous               = -1;
654:  #if defined(PETSC_USE_REAL_SINGLE)
655:   pseudo->fatol                        = 1.e-25;
656:   pseudo->frtol                        = 1.e-5;
657: #else
658:   pseudo->fatol                        = 1.e-50;
659:   pseudo->frtol                        = 1.e-12;
660: #endif
661:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
662:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
663:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
664:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
665:   PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
666:   return(0);
667: }

669: /*@C
670:    TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
671:    Use with TSPseudoSetTimeStep().

673:    Collective on TS

675:    Input Parameters:
676: +  ts - the timestep context
677: -  dtctx - unused timestep context

679:    Output Parameter:
680: .  newdt - the timestep to use for the next step

682:    Level: advanced

684: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
685: @*/
686: PetscErrorCode  TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
687: {
688:   TS_Pseudo      *pseudo = (TS_Pseudo*)ts->data;
689:   PetscReal      inc = pseudo->dt_increment;

693:   VecZeroEntries(pseudo->xdot);
694:   TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
695:   VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
696:   if (pseudo->fnorm_initial < 0) {
697:     /* first time through so compute initial function norm */
698:     pseudo->fnorm_initial  = pseudo->fnorm;
699:     pseudo->fnorm_previous = pseudo->fnorm;
700:   }
701:   if (pseudo->fnorm == 0.0)                      *newdt = 1.e12*inc*ts->time_step;
702:   else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
703:   else                                           *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
704:   if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
705:   pseudo->fnorm_previous = pseudo->fnorm;
706:   return(0);
707: }