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;
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: }
62: /* ------------------------------------------------------------------------------*/
63: /*@C
64: TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
66: Collective on TS
68: Input Parameters:
69: + ts - the timestep context
70: . dtctx - unused timestep context
71: - update - latest solution vector
73: Output Parameters:
74: + newdt - the timestep to use for the next step
75: - flag - flag indicating whether the last time step was acceptable
77: Level: advanced
79: Note:
80: This routine always returns a flag of 1, indicating an acceptable
81: timestep.
83: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
84: @*/
85: PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool *flag)
86: {
88: *flag = PETSC_TRUE;
89: return(0);
90: }
92: /*@
93: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
95: Collective on TS
97: Input Parameters:
98: + ts - timestep context
99: - update - latest solution vector
101: Output Parameters:
102: + dt - newly computed timestep (if it had to shrink)
103: - flag - indicates if current timestep was ok
105: Level: advanced
107: Notes:
108: The routine to be called here to compute the timestep should be
109: set by calling TSPseudoSetVerifyTimeStep().
111: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
112: @*/
113: PetscErrorCode TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
114: {
115: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
119: *flag = PETSC_TRUE;
120: if (pseudo->verify) {
121: (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
122: }
123: return(0);
124: }
126: /* --------------------------------------------------------------------------------*/
128: static PetscErrorCode TSStep_Pseudo(TS ts)
129: {
130: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
131: PetscInt nits,lits,reject;
132: PetscBool stepok;
133: PetscReal next_time_step = ts->time_step;
134: PetscErrorCode ierr;
137: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
138: VecCopy(ts->vec_sol,pseudo->update);
139: TSPseudoComputeTimeStep(ts,&next_time_step);
140: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
141: ts->time_step = next_time_step;
142: TSPreStage(ts,ts->ptime+ts->time_step);
143: SNESSolve(ts->snes,NULL,pseudo->update);
144: SNESGetIterationNumber(ts->snes,&nits);
145: SNESGetLinearSolveIterations(ts->snes,&lits);
146: ts->snes_its += nits; ts->ksp_its += lits;
147: TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));
148: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);
149: if (!stepok) {next_time_step = ts->time_step; continue;}
150: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
151: TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
152: if (stepok) break;
153: }
154: if (reject >= ts->max_reject) {
155: ts->reason = TS_DIVERGED_STEP_REJECTED;
156: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
157: return(0);
158: }
160: VecCopy(pseudo->update,ts->vec_sol);
161: ts->ptime += ts->time_step;
162: ts->time_step = next_time_step;
164: if (pseudo->fnorm < 0) {
165: VecZeroEntries(pseudo->xdot);
166: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
167: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
168: }
169: if (pseudo->fnorm < pseudo->fatol) {
170: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
171: PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
172: return(0);
173: }
174: if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
175: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
176: PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
177: return(0);
178: }
179: return(0);
180: }
182: /*------------------------------------------------------------*/
183: static PetscErrorCode TSReset_Pseudo(TS ts)
184: {
185: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
189: VecDestroy(&pseudo->update);
190: VecDestroy(&pseudo->func);
191: VecDestroy(&pseudo->xdot);
192: return(0);
193: }
195: static PetscErrorCode TSDestroy_Pseudo(TS ts)
196: {
200: TSReset_Pseudo(ts);
201: PetscFree(ts->data);
202: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
203: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
204: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
205: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
206: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
207: return(0);
208: }
210: /*------------------------------------------------------------*/
212: /*
213: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
214: */
215: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
216: {
217: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
218: const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
219: PetscScalar *xdot;
220: PetscErrorCode ierr;
221: PetscInt i,n;
224: *Xdot = NULL;
225: VecGetArrayRead(ts->vec_sol,&xn);
226: VecGetArrayRead(X,&xnp1);
227: VecGetArray(pseudo->xdot,&xdot);
228: VecGetLocalSize(X,&n);
229: for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
230: VecRestoreArrayRead(ts->vec_sol,&xn);
231: VecRestoreArrayRead(X,&xnp1);
232: VecRestoreArray(pseudo->xdot,&xdot);
233: *Xdot = pseudo->xdot;
234: return(0);
235: }
237: /*
238: The transient residual is
240: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
242: or for ODE,
244: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
246: This is the function that must be evaluated for transient simulation and for
247: finite difference Jacobians. On the first Newton step, this algorithm uses
248: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
249: residual is actually the steady state residual. Pseudotransient
250: continuation as described in the literature is a linearly implicit
251: algorithm, it only takes this one Newton step with the steady state
252: residual, and then advances to the next time step.
253: */
254: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
255: {
256: Vec Xdot;
260: TSPseudoGetXdot(ts,X,&Xdot);
261: TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
262: return(0);
263: }
265: /*
266: This constructs the Jacobian needed for SNES. For DAE, this is
268: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
270: and for ODE:
272: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
273: */
274: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
275: {
276: Vec Xdot;
280: TSPseudoGetXdot(ts,X,&Xdot);
281: TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
282: return(0);
283: }
285: static PetscErrorCode TSSetUp_Pseudo(TS ts)
286: {
287: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
291: VecDuplicate(ts->vec_sol,&pseudo->update);
292: VecDuplicate(ts->vec_sol,&pseudo->func);
293: VecDuplicate(ts->vec_sol,&pseudo->xdot);
294: return(0);
295: }
296: /*------------------------------------------------------------*/
298: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
299: {
300: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
302: PetscViewer viewer = (PetscViewer) dummy;
305: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
306: VecZeroEntries(pseudo->xdot);
307: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
308: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
309: }
310: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
311: PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
312: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
313: return(0);
314: }
316: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
317: {
318: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
320: PetscBool flg = PETSC_FALSE;
321: PetscViewer viewer;
324: PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
325: PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
326: if (flg) {
327: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
328: TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
329: }
330: flg = pseudo->increment_dt_from_initial_dt;
331: PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
332: pseudo->increment_dt_from_initial_dt = flg;
333: PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
334: PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
335: PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
336: PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
337: PetscOptionsTail();
338: return(0);
339: }
341: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
342: {
344: PetscBool isascii;
347: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
348: if (isascii) {
349: TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
350: PetscViewerASCIIPrintf(viewer," frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
351: PetscViewerASCIIPrintf(viewer," fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
352: PetscViewerASCIIPrintf(viewer," dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
353: PetscViewerASCIIPrintf(viewer," dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
354: PetscViewerASCIIPrintf(viewer," dt_max - maximum time %g\n",(double)pseudo->dt_max);
355: }
356: return(0);
357: }
359: /* ----------------------------------------------------------------------------- */
360: /*@C
361: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
362: last timestep.
364: Logically Collective on TS
366: Input Parameters:
367: + ts - timestep context
368: . dt - user-defined function to verify timestep
369: - ctx - [optional] user-defined context for private data
370: for the timestep verification routine (may be NULL)
372: Level: advanced
374: Calling sequence of func:
375: $ func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
377: + update - latest solution vector
378: . ctx - [optional] timestep context
379: . newdt - the timestep to use for the next step
380: - flag - flag indicating whether the last time step was acceptable
382: Notes:
383: The routine set here will be called by TSPseudoVerifyTimeStep()
384: during the timestepping process.
386: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
387: @*/
388: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
389: {
394: PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
395: return(0);
396: }
398: /*@
399: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
400: dt when using the TSPseudoTimeStepDefault() routine.
402: Logically Collective on TS
404: Input Parameters:
405: + ts - the timestep context
406: - inc - the scaling factor >= 1.0
408: Options Database Key:
409: . -ts_pseudo_increment <increment>
411: Level: advanced
413: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
414: @*/
415: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
416: {
422: PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
423: return(0);
424: }
426: /*@
427: TSPseudoSetMaxTimeStep - Sets the maximum time step
428: when using the TSPseudoTimeStepDefault() routine.
430: Logically Collective on TS
432: Input Parameters:
433: + ts - the timestep context
434: - maxdt - the maximum time step, use a non-positive value to deactivate
436: Options Database Key:
437: . -ts_pseudo_max_dt <increment>
439: Level: advanced
441: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
442: @*/
443: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
444: {
450: PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
451: return(0);
452: }
454: /*@
455: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
456: is computed via the formula
457: $ dt = initial_dt*initial_fnorm/current_fnorm
458: rather than the default update,
459: $ dt = current_dt*previous_fnorm/current_fnorm.
461: Logically Collective on TS
463: Input Parameter:
464: . ts - the timestep context
466: Options Database Key:
467: . -ts_pseudo_increment_dt_from_initial_dt
469: Level: advanced
471: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
472: @*/
473: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
474: {
479: PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
480: return(0);
481: }
483: /*@C
484: TSPseudoSetTimeStep - Sets the user-defined routine to be
485: called at each pseudo-timestep to update the timestep.
487: Logically Collective on TS
489: Input Parameters:
490: + ts - timestep context
491: . dt - function to compute timestep
492: - ctx - [optional] user-defined context for private data
493: required by the function (may be NULL)
495: Level: intermediate
497: Calling sequence of func:
498: $ func (TS ts,PetscReal *newdt,void *ctx);
500: + newdt - the newly computed timestep
501: - ctx - [optional] timestep context
503: Notes:
504: The routine set here will be called by TSPseudoComputeTimeStep()
505: during the timestepping process.
506: If not set then TSPseudoTimeStepDefault() is automatically used
508: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
509: @*/
510: PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
511: {
516: PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
517: return(0);
518: }
520: /* ----------------------------------------------------------------------------- */
522: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/
523: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
524: {
525: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
528: pseudo->verify = dt;
529: pseudo->verifyctx = ctx;
530: return(0);
531: }
533: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
534: {
535: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
538: pseudo->dt_increment = inc;
539: return(0);
540: }
542: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
543: {
544: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
547: pseudo->dt_max = maxdt;
548: return(0);
549: }
551: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
552: {
553: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
556: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
557: return(0);
558: }
560: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
561: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
562: {
563: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
566: pseudo->dt = dt;
567: pseudo->dtctx = ctx;
568: return(0);
569: }
571: /* ----------------------------------------------------------------------------- */
572: /*MC
573: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
575: This method solves equations of the form
577: $ F(X,Xdot) = 0
579: for steady state using the iteration
581: $ [G'] S = -F(X,0)
582: $ X += S
584: where
586: $ G(Y) = F(Y,(Y-X)/dt)
588: This is linearly-implicit Euler with the residual always evaluated "at steady
589: state". See note below.
591: Options database keys:
592: + -ts_pseudo_increment <real> - ratio of increase dt
593: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
594: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
595: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
597: Level: beginner
599: References:
600: + 1. - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
601: - 2. - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
603: Notes:
604: The residual computed by this method includes the transient term (Xdot is computed instead of
605: always being zero), but since the prediction from the last step is always the solution from the
606: last step, on the first Newton iteration we have
608: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
610: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
611: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
612: algorithm is no longer the one described in the referenced papers.
614: .seealso: TSCreate(), TS, TSSetType()
616: M*/
617: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
618: {
619: TS_Pseudo *pseudo;
621: SNES snes;
622: SNESType stype;
625: ts->ops->reset = TSReset_Pseudo;
626: ts->ops->destroy = TSDestroy_Pseudo;
627: ts->ops->view = TSView_Pseudo;
628: ts->ops->setup = TSSetUp_Pseudo;
629: ts->ops->step = TSStep_Pseudo;
630: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
631: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
632: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
633: ts->default_adapt_type = TSADAPTNONE;
634: ts->usessnes = PETSC_TRUE;
636: TSGetSNES(ts,&snes);
637: SNESGetType(snes,&stype);
638: if (!stype) {SNESSetType(snes,SNESKSPONLY);}
640: PetscNewLog(ts,&pseudo);
641: ts->data = (void*)pseudo;
643: pseudo->dt = TSPseudoTimeStepDefault;
644: pseudo->dtctx = NULL;
645: pseudo->dt_increment = 1.1;
646: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
647: pseudo->fnorm = -1;
648: pseudo->fnorm_initial = -1;
649: pseudo->fnorm_previous = -1;
650: #if defined(PETSC_USE_REAL_SINGLE)
651: pseudo->fatol = 1.e-25;
652: pseudo->frtol = 1.e-5;
653: #else
654: pseudo->fatol = 1.e-50;
655: pseudo->frtol = 1.e-12;
656: #endif
657: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
658: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
659: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
660: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
661: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
662: return(0);
663: }
665: /*@C
666: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
667: Use with TSPseudoSetTimeStep().
669: Collective on TS
671: Input Parameters:
672: + ts - the timestep context
673: - dtctx - unused timestep context
675: Output Parameter:
676: . newdt - the timestep to use for the next step
678: Level: advanced
680: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
681: @*/
682: PetscErrorCode TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
683: {
684: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
685: PetscReal inc = pseudo->dt_increment;
689: VecZeroEntries(pseudo->xdot);
690: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
691: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
692: if (pseudo->fnorm_initial < 0) {
693: /* first time through so compute initial function norm */
694: pseudo->fnorm_initial = pseudo->fnorm;
695: pseudo->fnorm_previous = pseudo->fnorm;
696: }
697: if (pseudo->fnorm == 0.0) *newdt = 1.e12*inc*ts->time_step;
698: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
699: else *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
700: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
701: pseudo->fnorm_previous = pseudo->fnorm;
702: return(0);
703: }