Actual source code: posindep.c

  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;

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

 60: /* ------------------------------------------------------------------------------*/
 61: /*@C
 62:    TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.

 64:    Collective on TS

 66:    Input Parameters:
 67: +  ts - the timestep context
 68: .  dtctx - unused timestep context
 69: -  update - latest solution vector

 71:    Output Parameters:
 72: +  newdt - the timestep to use for the next step
 73: -  flag - flag indicating whether the last time step was acceptable

 75:    Level: advanced

 77:    Note:
 78:    This routine always returns a flag of 1, indicating an acceptable
 79:    timestep.

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

 89: /*@
 90:     TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.

 92:     Collective on TS

 94:     Input Parameters:
 95: +   ts - timestep context
 96: -   update - latest solution vector

 98:     Output Parameters:
 99: +   dt - newly computed timestep (if it had to shrink)
100: -   flag - indicates if current timestep was ok

102:     Level: advanced

104:     Notes:
105:     The routine to be called here to compute the timestep should be
106:     set by calling TSPseudoSetVerifyTimeStep().

108: .seealso: `TSPseudoSetVerifyTimeStep()`, `TSPseudoVerifyTimeStepDefault()`
109: @*/
110: PetscErrorCode TSPseudoVerifyTimeStep(TS ts, Vec update, PetscReal *dt, PetscBool *flag)
111: {
112:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

114:   *flag = PETSC_TRUE;
115:   if (pseudo->verify) (*pseudo->verify)(ts, update, pseudo->verifyctx, dt, flag);
116:   return 0;
117: }

119: /* --------------------------------------------------------------------------------*/

121: static PetscErrorCode TSStep_Pseudo(TS ts)
122: {
123:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
124:   PetscInt   nits, lits, reject;
125:   PetscBool  stepok;
126:   PetscReal  next_time_step = ts->time_step;

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

155:   VecCopy(pseudo->update, ts->vec_sol);
156:   ts->ptime += ts->time_step;
157:   ts->time_step = next_time_step;

159:   if (pseudo->fnorm < 0) {
160:     VecZeroEntries(pseudo->xdot);
161:     TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE);
162:     VecNorm(pseudo->func, NORM_2, &pseudo->fnorm);
163:   }
164:   if (pseudo->fnorm < pseudo->fatol) {
165:     ts->reason = TS_CONVERGED_PSEUDO_FATOL;
166:     PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->frtol);
167:     return 0;
168:   }
169:   if (pseudo->fnorm / pseudo->fnorm_initial < pseudo->frtol) {
170:     ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
171:     PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g / fnorm_initial %g < frtol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->fnorm_initial, (double)pseudo->fatol);
172:     return 0;
173:   }
174:   return 0;
175: }

177: /*------------------------------------------------------------*/
178: static PetscErrorCode TSReset_Pseudo(TS ts)
179: {
180:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

182:   VecDestroy(&pseudo->update);
183:   VecDestroy(&pseudo->func);
184:   VecDestroy(&pseudo->xdot);
185:   return 0;
186: }

188: static PetscErrorCode TSDestroy_Pseudo(TS ts)
189: {
190:   TSReset_Pseudo(ts);
191:   PetscFree(ts->data);
192:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL);
193:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL);
194:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL);
195:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL);
196:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL);
197:   return 0;
198: }

200: /*------------------------------------------------------------*/

202: /*
203:     Compute Xdot = (X^{n+1}-X^n)/dt) = 0
204: */
205: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
206: {
207:   TS_Pseudo        *pseudo = (TS_Pseudo *)ts->data;
208:   const PetscScalar mdt    = 1.0 / ts->time_step, *xnp1, *xn;
209:   PetscScalar      *xdot;
210:   PetscInt          i, n;

212:   *Xdot = NULL;
213:   VecGetArrayRead(ts->vec_sol, &xn);
214:   VecGetArrayRead(X, &xnp1);
215:   VecGetArray(pseudo->xdot, &xdot);
216:   VecGetLocalSize(X, &n);
217:   for (i = 0; i < n; i++) xdot[i] = mdt * (xnp1[i] - xn[i]);
218:   VecRestoreArrayRead(ts->vec_sol, &xn);
219:   VecRestoreArrayRead(X, &xnp1);
220:   VecRestoreArray(pseudo->xdot, &xdot);
221:   *Xdot = pseudo->xdot;
222:   return 0;
223: }

225: /*
226:     The transient residual is

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

230:     or for ODE,

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

234:     This is the function that must be evaluated for transient simulation and for
235:     finite difference Jacobians.  On the first Newton step, this algorithm uses
236:     a guess of U^{n+1} = U^n in which case the transient term vanishes and the
237:     residual is actually the steady state residual.  Pseudotransient
238:     continuation as described in the literature is a linearly implicit
239:     algorithm, it only takes this one Newton step with the steady state
240:     residual, and then advances to the next time step.
241: */
242: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
243: {
244:   Vec Xdot;

246:   TSPseudoGetXdot(ts, X, &Xdot);
247:   TSComputeIFunction(ts, ts->ptime + ts->time_step, X, Xdot, Y, PETSC_FALSE);
248:   return 0;
249: }

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

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

256:     and for ODE:

258:        J = I/dt - J_{Frhs}   where J_{Frhs} is the given Jacobian of Frhs.
259: */
260: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
261: {
262:   Vec Xdot;

264:   TSPseudoGetXdot(ts, X, &Xdot);
265:   TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE);
266:   return 0;
267: }

269: static PetscErrorCode TSSetUp_Pseudo(TS ts)
270: {
271:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

273:   VecDuplicate(ts->vec_sol, &pseudo->update);
274:   VecDuplicate(ts->vec_sol, &pseudo->func);
275:   VecDuplicate(ts->vec_sol, &pseudo->xdot);
276:   return 0;
277: }
278: /*------------------------------------------------------------*/

280: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
281: {
282:   TS_Pseudo  *pseudo = (TS_Pseudo *)ts->data;
283:   PetscViewer viewer = (PetscViewer)dummy;

285:   if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
286:     VecZeroEntries(pseudo->xdot);
287:     TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE);
288:     VecNorm(pseudo->func, NORM_2, &pseudo->fnorm);
289:   }
290:   PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel);
291:   PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)pseudo->fnorm);
292:   PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel);
293:   return 0;
294: }

296: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems *PetscOptionsObject)
297: {
298:   TS_Pseudo  *pseudo = (TS_Pseudo *)ts->data;
299:   PetscBool   flg    = PETSC_FALSE;
300:   PetscViewer viewer;

302:   PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
303:   PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL);
304:   if (flg) {
305:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer);
306:     TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscErrorCode(*)(void **))PetscViewerDestroy);
307:   }
308:   flg = pseudo->increment_dt_from_initial_dt;
309:   PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL);
310:   pseudo->increment_dt_from_initial_dt = flg;
311:   PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL);
312:   PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL);
313:   PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL);
314:   PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL);
315:   PetscOptionsHeadEnd();
316:   return 0;
317: }

319: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
320: {
321:   PetscBool isascii;

323:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
324:   if (isascii) {
325:     TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
326:     PetscViewerASCIIPrintf(viewer, "  frtol - relative tolerance in function value %g\n", (double)pseudo->frtol);
327:     PetscViewerASCIIPrintf(viewer, "  fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol);
328:     PetscViewerASCIIPrintf(viewer, "  dt_initial - initial timestep %g\n", (double)pseudo->dt_initial);
329:     PetscViewerASCIIPrintf(viewer, "  dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment);
330:     PetscViewerASCIIPrintf(viewer, "  dt_max - maximum time %g\n", (double)pseudo->dt_max);
331:   }
332:   return 0;
333: }

335: /* ----------------------------------------------------------------------------- */
336: /*@C
337:    TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
338:    last timestep.

340:    Logically Collective on TS

342:    Input Parameters:
343: +  ts - timestep context
344: .  dt - user-defined function to verify timestep
345: -  ctx - [optional] user-defined context for private data
346:          for the timestep verification routine (may be NULL)

348:    Level: advanced

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

353: +  update - latest solution vector
354: .  ctx - [optional] timestep context
355: .  newdt - the timestep to use for the next step
356: -  flag - flag indicating whether the last time step was acceptable

358:    Notes:
359:    The routine set here will be called by TSPseudoVerifyTimeStep()
360:    during the timestepping process.

362: .seealso: `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
363: @*/
364: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS, Vec, void *, PetscReal *, PetscBool *), void *ctx)
365: {
367:   PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode(*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
368:   return 0;
369: }

371: /*@
372:     TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
373:     dt when using the TSPseudoTimeStepDefault() routine.

375:    Logically Collective on TS

377:     Input Parameters:
378: +   ts - the timestep context
379: -   inc - the scaling factor >= 1.0

381:     Options Database Key:
382: .    -ts_pseudo_increment <increment> - set pseudo increment

384:     Level: advanced

386: .seealso: `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
387: @*/
388: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
389: {
392:   PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
393:   return 0;
394: }

396: /*@
397:     TSPseudoSetMaxTimeStep - Sets the maximum time step
398:     when using the TSPseudoTimeStepDefault() routine.

400:    Logically Collective on TS

402:     Input Parameters:
403: +   ts - the timestep context
404: -   maxdt - the maximum time step, use a non-positive value to deactivate

406:     Options Database Key:
407: .    -ts_pseudo_max_dt <increment> - set pseudo max dt

409:     Level: advanced

411: .seealso: `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
412: @*/
413: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
414: {
417:   PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
418:   return 0;
419: }

421: /*@
422:     TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
423:     is computed via the formula
424: $         dt = initial_dt*initial_fnorm/current_fnorm
425:       rather than the default update,
426: $         dt = current_dt*previous_fnorm/current_fnorm.

428:    Logically Collective on TS

430:     Input Parameter:
431: .   ts - the timestep context

433:     Options Database Key:
434: .    -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial dt to determine increment

436:     Level: advanced

438: .seealso: `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
439: @*/
440: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
441: {
443:   PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
444:   return 0;
445: }

447: /*@C
448:    TSPseudoSetTimeStep - Sets the user-defined routine to be
449:    called at each pseudo-timestep to update the timestep.

451:    Logically Collective on TS

453:    Input Parameters:
454: +  ts - timestep context
455: .  dt - function to compute timestep
456: -  ctx - [optional] user-defined context for private data
457:          required by the function (may be NULL)

459:    Level: intermediate

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

464: +  newdt - the newly computed timestep
465: -  ctx - [optional] timestep context

467:    Notes:
468:    The routine set here will be called by TSPseudoComputeTimeStep()
469:    during the timestepping process.
470:    If not set then TSPseudoTimeStepDefault() is automatically used

472: .seealso: `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
473: @*/
474: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS, PetscReal *, void *), void *ctx)
475: {
477:   PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode(*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
478:   return 0;
479: }

481: /* ----------------------------------------------------------------------------- */

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

488:   pseudo->verify    = dt;
489:   pseudo->verifyctx = ctx;
490:   return 0;
491: }

493: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
494: {
495:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

497:   pseudo->dt_increment = inc;
498:   return 0;
499: }

501: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
502: {
503:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

505:   pseudo->dt_max = maxdt;
506:   return 0;
507: }

509: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
510: {
511:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;

513:   pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
514:   return 0;
515: }

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

522:   pseudo->dt    = dt;
523:   pseudo->dtctx = ctx;
524:   return 0;
525: }

527: /* ----------------------------------------------------------------------------- */
528: /*MC
529:       TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping

531:   This method solves equations of the form

533: $    F(X,Xdot) = 0

535:   for steady state using the iteration

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

540:   where

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

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

547:   Options database keys:
548: +  -ts_pseudo_increment <real> - ratio of increase dt
549: .  -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
550: .  -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
551: -  -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol

553:   Level: beginner

555:   References:
556: +  * - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
557: -  * - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.

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

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

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

570: .seealso: `TSCreate()`, `TS`, `TSSetType()`

572: M*/
573: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
574: {
575:   TS_Pseudo *pseudo;
576:   SNES       snes;
577:   SNESType   stype;

579:   ts->ops->reset          = TSReset_Pseudo;
580:   ts->ops->destroy        = TSDestroy_Pseudo;
581:   ts->ops->view           = TSView_Pseudo;
582:   ts->ops->setup          = TSSetUp_Pseudo;
583:   ts->ops->step           = TSStep_Pseudo;
584:   ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
585:   ts->ops->snesfunction   = SNESTSFormFunction_Pseudo;
586:   ts->ops->snesjacobian   = SNESTSFormJacobian_Pseudo;
587:   ts->default_adapt_type  = TSADAPTNONE;
588:   ts->usessnes            = PETSC_TRUE;

590:   TSGetSNES(ts, &snes);
591:   SNESGetType(snes, &stype);
592:   if (!stype) SNESSetType(snes, SNESKSPONLY);

594:   PetscNew(&pseudo);
595:   ts->data = (void *)pseudo;

597:   pseudo->dt                           = TSPseudoTimeStepDefault;
598:   pseudo->dtctx                        = NULL;
599:   pseudo->dt_increment                 = 1.1;
600:   pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
601:   pseudo->fnorm                        = -1;
602:   pseudo->fnorm_initial                = -1;
603:   pseudo->fnorm_previous               = -1;
604: #if defined(PETSC_USE_REAL_SINGLE)
605:   pseudo->fatol = 1.e-25;
606:   pseudo->frtol = 1.e-5;
607: #else
608:   pseudo->fatol = 1.e-50;
609:   pseudo->frtol = 1.e-12;
610: #endif
611:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo);
612:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo);
613:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo);
614:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo);
615:   PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo);
616:   return 0;
617: }

619: /*@C
620:    TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
621:    Use with TSPseudoSetTimeStep().

623:    Collective on TS

625:    Input Parameters:
626: +  ts - the timestep context
627: -  dtctx - unused timestep context

629:    Output Parameter:
630: .  newdt - the timestep to use for the next step

632:    Level: advanced

634: .seealso: `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`
635: @*/
636: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
637: {
638:   TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
639:   PetscReal  inc    = pseudo->dt_increment;

641:   VecZeroEntries(pseudo->xdot);
642:   TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE);
643:   VecNorm(pseudo->func, NORM_2, &pseudo->fnorm);
644:   if (pseudo->fnorm_initial < 0) {
645:     /* first time through so compute initial function norm */
646:     pseudo->fnorm_initial  = pseudo->fnorm;
647:     pseudo->fnorm_previous = pseudo->fnorm;
648:   }
649:   if (pseudo->fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
650:   else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / pseudo->fnorm;
651:   else *newdt = inc * ts->time_step * pseudo->fnorm_previous / pseudo->fnorm;
652:   if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
653:   pseudo->fnorm_previous = pseudo->fnorm;
654:   return 0;
655: }