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) {
116: (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
117: }
118: return 0;
119: }
121: /* --------------------------------------------------------------------------------*/
123: static PetscErrorCode TSStep_Pseudo(TS ts)
124: {
125: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
126: PetscInt nits,lits,reject;
127: PetscBool stepok;
128: PetscReal next_time_step = ts->time_step;
130: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
131: VecCopy(ts->vec_sol,pseudo->update);
132: TSPseudoComputeTimeStep(ts,&next_time_step);
133: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
134: ts->time_step = next_time_step;
135: TSPreStage(ts,ts->ptime+ts->time_step);
136: SNESSolve(ts->snes,NULL,pseudo->update);
137: SNESGetIterationNumber(ts->snes,&nits);
138: SNESGetLinearSolveIterations(ts->snes,&lits);
139: ts->snes_its += nits; ts->ksp_its += lits;
140: TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));
141: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);
142: if (!stepok) {next_time_step = ts->time_step; continue;}
143: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
144: TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
145: if (stepok) break;
146: }
147: if (reject >= ts->max_reject) {
148: ts->reason = TS_DIVERGED_STEP_REJECTED;
149: PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
150: return 0;
151: }
153: VecCopy(pseudo->update,ts->vec_sol);
154: ts->ptime += ts->time_step;
155: ts->time_step = next_time_step;
157: if (pseudo->fnorm < 0) {
158: VecZeroEntries(pseudo->xdot);
159: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
160: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
161: }
162: if (pseudo->fnorm < pseudo->fatol) {
163: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
164: PetscInfo(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
165: return 0;
166: }
167: if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
168: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
169: PetscInfo(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
170: return 0;
171: }
172: return 0;
173: }
175: /*------------------------------------------------------------*/
176: static PetscErrorCode TSReset_Pseudo(TS ts)
177: {
178: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
180: VecDestroy(&pseudo->update);
181: VecDestroy(&pseudo->func);
182: VecDestroy(&pseudo->xdot);
183: return 0;
184: }
186: static PetscErrorCode TSDestroy_Pseudo(TS ts)
187: {
188: TSReset_Pseudo(ts);
189: PetscFree(ts->data);
190: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
191: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
192: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
193: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
194: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
195: return 0;
196: }
198: /*------------------------------------------------------------*/
200: /*
201: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
202: */
203: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
204: {
205: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
206: const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
207: PetscScalar *xdot;
208: PetscInt i,n;
210: *Xdot = NULL;
211: VecGetArrayRead(ts->vec_sol,&xn);
212: VecGetArrayRead(X,&xnp1);
213: VecGetArray(pseudo->xdot,&xdot);
214: VecGetLocalSize(X,&n);
215: for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
216: VecRestoreArrayRead(ts->vec_sol,&xn);
217: VecRestoreArrayRead(X,&xnp1);
218: VecRestoreArray(pseudo->xdot,&xdot);
219: *Xdot = pseudo->xdot;
220: return 0;
221: }
223: /*
224: The transient residual is
226: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
228: or for ODE,
230: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
232: This is the function that must be evaluated for transient simulation and for
233: finite difference Jacobians. On the first Newton step, this algorithm uses
234: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
235: residual is actually the steady state residual. Pseudotransient
236: continuation as described in the literature is a linearly implicit
237: algorithm, it only takes this one Newton step with the steady state
238: residual, and then advances to the next time step.
239: */
240: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
241: {
242: Vec Xdot;
244: TSPseudoGetXdot(ts,X,&Xdot);
245: TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
246: return 0;
247: }
249: /*
250: This constructs the Jacobian needed for SNES. For DAE, this is
252: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
254: and for ODE:
256: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
257: */
258: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
259: {
260: Vec Xdot;
262: TSPseudoGetXdot(ts,X,&Xdot);
263: TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
264: return 0;
265: }
267: static PetscErrorCode TSSetUp_Pseudo(TS ts)
268: {
269: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
271: VecDuplicate(ts->vec_sol,&pseudo->update);
272: VecDuplicate(ts->vec_sol,&pseudo->func);
273: VecDuplicate(ts->vec_sol,&pseudo->xdot);
274: return 0;
275: }
276: /*------------------------------------------------------------*/
278: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
279: {
280: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
281: PetscViewer viewer = (PetscViewer) dummy;
283: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
284: VecZeroEntries(pseudo->xdot);
285: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
286: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
287: }
288: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
289: PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
290: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
291: return 0;
292: }
294: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
295: {
296: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
297: PetscBool flg = PETSC_FALSE;
298: PetscViewer viewer;
300: PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
301: PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
302: if (flg) {
303: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
304: TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
305: }
306: flg = pseudo->increment_dt_from_initial_dt;
307: PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
308: pseudo->increment_dt_from_initial_dt = flg;
309: PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
310: PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
311: PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
312: PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
313: PetscOptionsTail();
314: return 0;
315: }
317: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
318: {
319: PetscBool isascii;
321: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
322: if (isascii) {
323: TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
324: PetscViewerASCIIPrintf(viewer," frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
325: PetscViewerASCIIPrintf(viewer," fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
326: PetscViewerASCIIPrintf(viewer," dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
327: PetscViewerASCIIPrintf(viewer," dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
328: PetscViewerASCIIPrintf(viewer," dt_max - maximum time %g\n",(double)pseudo->dt_max);
329: }
330: return 0;
331: }
333: /* ----------------------------------------------------------------------------- */
334: /*@C
335: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
336: last timestep.
338: Logically Collective on TS
340: Input Parameters:
341: + ts - timestep context
342: . dt - user-defined function to verify timestep
343: - ctx - [optional] user-defined context for private data
344: for the timestep verification routine (may be NULL)
346: Level: advanced
348: Calling sequence of func:
349: $ func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
351: + update - latest solution vector
352: . ctx - [optional] timestep context
353: . newdt - the timestep to use for the next step
354: - flag - flag indicating whether the last time step was acceptable
356: Notes:
357: The routine set here will be called by TSPseudoVerifyTimeStep()
358: during the timestepping process.
360: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
361: @*/
362: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
363: {
365: PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
366: return 0;
367: }
369: /*@
370: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
371: dt when using the TSPseudoTimeStepDefault() routine.
373: Logically Collective on TS
375: Input Parameters:
376: + ts - the timestep context
377: - inc - the scaling factor >= 1.0
379: Options Database Key:
380: . -ts_pseudo_increment <increment> - set pseudo increment
382: Level: advanced
384: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
385: @*/
386: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
387: {
390: PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
391: return 0;
392: }
394: /*@
395: TSPseudoSetMaxTimeStep - Sets the maximum time step
396: when using the TSPseudoTimeStepDefault() routine.
398: Logically Collective on TS
400: Input Parameters:
401: + ts - the timestep context
402: - maxdt - the maximum time step, use a non-positive value to deactivate
404: Options Database Key:
405: . -ts_pseudo_max_dt <increment> - set pseudo max dt
407: Level: advanced
409: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
410: @*/
411: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
412: {
415: PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
416: return 0;
417: }
419: /*@
420: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
421: is computed via the formula
422: $ dt = initial_dt*initial_fnorm/current_fnorm
423: rather than the default update,
424: $ dt = current_dt*previous_fnorm/current_fnorm.
426: Logically Collective on TS
428: Input Parameter:
429: . ts - the timestep context
431: Options Database Key:
432: . -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial dt to determine increment
434: Level: advanced
436: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
437: @*/
438: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
439: {
441: PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
442: return 0;
443: }
445: /*@C
446: TSPseudoSetTimeStep - Sets the user-defined routine to be
447: called at each pseudo-timestep to update the timestep.
449: Logically Collective on TS
451: Input Parameters:
452: + ts - timestep context
453: . dt - function to compute timestep
454: - ctx - [optional] user-defined context for private data
455: required by the function (may be NULL)
457: Level: intermediate
459: Calling sequence of func:
460: $ func (TS ts,PetscReal *newdt,void *ctx);
462: + newdt - the newly computed timestep
463: - ctx - [optional] timestep context
465: Notes:
466: The routine set here will be called by TSPseudoComputeTimeStep()
467: during the timestepping process.
468: If not set then TSPseudoTimeStepDefault() is automatically used
470: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
471: @*/
472: PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
473: {
475: PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
476: return 0;
477: }
479: /* ----------------------------------------------------------------------------- */
481: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/
482: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
483: {
484: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
486: pseudo->verify = dt;
487: pseudo->verifyctx = ctx;
488: return 0;
489: }
491: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
492: {
493: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
495: pseudo->dt_increment = inc;
496: return 0;
497: }
499: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
500: {
501: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
503: pseudo->dt_max = maxdt;
504: return 0;
505: }
507: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
508: {
509: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
511: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
512: return 0;
513: }
515: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
516: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
517: {
518: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
520: pseudo->dt = dt;
521: pseudo->dtctx = ctx;
522: return 0;
523: }
525: /* ----------------------------------------------------------------------------- */
526: /*MC
527: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
529: This method solves equations of the form
531: $ F(X,Xdot) = 0
533: for steady state using the iteration
535: $ [G'] S = -F(X,0)
536: $ X += S
538: where
540: $ G(Y) = F(Y,(Y-X)/dt)
542: This is linearly-implicit Euler with the residual always evaluated "at steady
543: state". See note below.
545: Options database keys:
546: + -ts_pseudo_increment <real> - ratio of increase dt
547: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
548: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
549: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
551: Level: beginner
553: References:
554: + * - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
555: - * - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
557: Notes:
558: The residual computed by this method includes the transient term (Xdot is computed instead of
559: always being zero), but since the prediction from the last step is always the solution from the
560: last step, on the first Newton iteration we have
562: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
564: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
565: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
566: algorithm is no longer the one described in the referenced papers.
568: .seealso: TSCreate(), TS, TSSetType()
570: M*/
571: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
572: {
573: TS_Pseudo *pseudo;
574: SNES snes;
575: SNESType stype;
577: ts->ops->reset = TSReset_Pseudo;
578: ts->ops->destroy = TSDestroy_Pseudo;
579: ts->ops->view = TSView_Pseudo;
580: ts->ops->setup = TSSetUp_Pseudo;
581: ts->ops->step = TSStep_Pseudo;
582: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
583: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
584: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
585: ts->default_adapt_type = TSADAPTNONE;
586: ts->usessnes = PETSC_TRUE;
588: TSGetSNES(ts,&snes);
589: SNESGetType(snes,&stype);
590: if (!stype) SNESSetType(snes,SNESKSPONLY);
592: PetscNewLog(ts,&pseudo);
593: ts->data = (void*)pseudo;
595: pseudo->dt = TSPseudoTimeStepDefault;
596: pseudo->dtctx = NULL;
597: pseudo->dt_increment = 1.1;
598: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
599: pseudo->fnorm = -1;
600: pseudo->fnorm_initial = -1;
601: pseudo->fnorm_previous = -1;
602: #if defined(PETSC_USE_REAL_SINGLE)
603: pseudo->fatol = 1.e-25;
604: pseudo->frtol = 1.e-5;
605: #else
606: pseudo->fatol = 1.e-50;
607: pseudo->frtol = 1.e-12;
608: #endif
609: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
610: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
611: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
612: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
613: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
614: return 0;
615: }
617: /*@C
618: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
619: Use with TSPseudoSetTimeStep().
621: Collective on TS
623: Input Parameters:
624: + ts - the timestep context
625: - dtctx - unused timestep context
627: Output Parameter:
628: . newdt - the timestep to use for the next step
630: Level: advanced
632: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
633: @*/
634: PetscErrorCode TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
635: {
636: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
637: PetscReal inc = pseudo->dt_increment;
639: VecZeroEntries(pseudo->xdot);
640: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
641: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
642: if (pseudo->fnorm_initial < 0) {
643: /* first time through so compute initial function norm */
644: pseudo->fnorm_initial = pseudo->fnorm;
645: pseudo->fnorm_previous = pseudo->fnorm;
646: }
647: if (pseudo->fnorm == 0.0) *newdt = 1.e12*inc*ts->time_step;
648: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
649: else *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
650: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
651: pseudo->fnorm_previous = pseudo->fnorm;
652: return 0;
653: }