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: }
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: }