Actual source code: posindep.c
petsc-3.3-p7 2013-05-11
1: /*
2: Code for Timestepping with implicit backwards Euler.
3: */
4: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
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: } TS_Pseudo;
27: /* ------------------------------------------------------------------------------*/
31: /*@
32: TSPseudoComputeTimeStep - Computes the next timestep for a currently running
33: pseudo-timestepping process.
35: Collective on TS
37: Input Parameter:
38: . ts - timestep context
40: Output Parameter:
41: . dt - newly computed timestep
43: Level: advanced
45: Notes:
46: The routine to be called here to compute the timestep should be
47: set by calling TSPseudoSetTimeStep().
49: .keywords: timestep, pseudo, compute
51: .seealso: TSPseudoDefaultTimeStep(), TSPseudoSetTimeStep()
52: @*/
53: PetscErrorCode TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
54: {
55: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
59: PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);
60: (*pseudo->dt)(ts,dt,pseudo->dtctx);
61: PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);
62: return(0);
63: }
66: /* ------------------------------------------------------------------------------*/
69: /*@C
70: TSPseudoDefaultVerifyTimeStep - Default code to verify the quality of the last timestep.
72: Collective on TS
74: Input Parameters:
75: + ts - the timestep context
76: . dtctx - unused timestep context
77: - update - latest solution vector
79: Output Parameters:
80: + newdt - the timestep to use for the next step
81: - flag - flag indicating whether the last time step was acceptable
83: Level: advanced
85: Note:
86: This routine always returns a flag of 1, indicating an acceptable
87: timestep.
89: .keywords: timestep, pseudo, default, verify
91: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
92: @*/
93: PetscErrorCode TSPseudoDefaultVerifyTimeStep(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool *flag)
94: {
96: *flag = PETSC_TRUE;
97: return(0);
98: }
103: /*@
104: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
106: Collective on TS
108: Input Parameters:
109: + ts - timestep context
110: - update - latest solution vector
112: Output Parameters:
113: + dt - newly computed timestep (if it had to shrink)
114: - flag - indicates if current timestep was ok
116: Level: advanced
118: Notes:
119: The routine to be called here to compute the timestep should be
120: set by calling TSPseudoSetVerifyTimeStep().
122: .keywords: timestep, pseudo, verify
124: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoDefaultVerifyTimeStep()
125: @*/
126: PetscErrorCode TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
127: {
128: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
132: if (!pseudo->verify) {*flag = PETSC_TRUE; return(0);}
134: (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
136: return(0);
137: }
139: /* --------------------------------------------------------------------------------*/
143: static PetscErrorCode TSStep_Pseudo(TS ts)
144: {
145: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
146: PetscInt its,lits,reject;
147: PetscBool stepok;
148: PetscReal next_time_step;
149: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
153: if (ts->steps == 0) {
154: pseudo->dt_initial = ts->time_step;
155: }
156: VecCopy(ts->vec_sol,pseudo->update);
157: next_time_step = ts->time_step;
158: TSPseudoComputeTimeStep(ts,&next_time_step);
159: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
160: ts->time_step = next_time_step;
161: TSPreStep(ts);
162: TSPreStage(ts,ts->ptime+ts->time_step);
163: SNESSolve(ts->snes,PETSC_NULL,pseudo->update);
164: SNESGetConvergedReason(ts->snes,&snesreason);
165: SNESGetLinearSolveIterations(ts->snes,&lits);
166: SNESGetIterationNumber(ts->snes,&its);
167: ts->snes_its += its; ts->ksp_its += lits;
168: PetscInfo3(ts,"step=%D, nonlinear solve iterations=%D, linear solve iterations=%D\n",ts->steps,its,lits);
169: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
170: TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
171: if (stepok) break;
172: }
173: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
174: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
175: PetscInfo2(ts,"step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
176: return(0);
177: }
178: if (reject >= ts->max_reject) {
179: ts->reason = TS_DIVERGED_STEP_REJECTED;
180: PetscInfo2(ts,"step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
181: return(0);
182: }
183: VecCopy(pseudo->update,ts->vec_sol);
184: ts->ptime += ts->time_step;
185: ts->time_step = next_time_step;
186: ts->steps++;
187: return(0);
188: }
190: /*------------------------------------------------------------*/
193: static PetscErrorCode TSReset_Pseudo(TS ts)
194: {
195: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
199: VecDestroy(&pseudo->update);
200: VecDestroy(&pseudo->func);
201: VecDestroy(&pseudo->xdot);
202: return(0);
203: }
207: static PetscErrorCode TSDestroy_Pseudo(TS ts)
208: {
212: TSReset_Pseudo(ts);
213: PetscFree(ts->data);
214: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C","",PETSC_NULL);
215: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C","",PETSC_NULL);
216: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetMaxTimeStep_C","",PETSC_NULL);
217: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C","",PETSC_NULL);
218: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStep_C","",PETSC_NULL);
219: return(0);
220: }
222: /*------------------------------------------------------------*/
226: /*
227: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
228: */
229: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
230: {
231: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
232: const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
233: PetscScalar *xdot;
235: PetscInt i,n;
238: VecGetArrayRead(ts->vec_sol,&xn);
239: VecGetArrayRead(X,&xnp1);
240: VecGetArray(pseudo->xdot,&xdot);
241: VecGetLocalSize(X,&n);
242: for (i=0; i<n; i++) {
243: xdot[i] = mdt*(xnp1[i] - xn[i]);
244: }
245: VecRestoreArrayRead(ts->vec_sol,&xn);
246: VecRestoreArrayRead(X,&xnp1);
247: VecRestoreArray(pseudo->xdot,&xdot);
248: *Xdot = pseudo->xdot;
249: return(0);
250: }
254: /*
255: The transient residual is
257: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
259: or for ODE,
261: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
263: This is the function that must be evaluated for transient simulation and for
264: finite difference Jacobians. On the first Newton step, this algorithm uses
265: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
266: residual is actually the steady state residual. Pseudotransient
267: continuation as described in the literature is a linearly implicit
268: algorithm, it only takes this one Newton step with the steady state
269: residual, and then advances to the next time step.
270: */
271: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
272: {
273: Vec Xdot;
277: TSPseudoGetXdot(ts,X,&Xdot);
278: TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
279: return(0);
280: }
284: /*
285: This constructs the Jacobian needed for SNES. For DAE, this is
287: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
289: and for ODE:
291: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
292: */
293: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat *AA,Mat *BB,MatStructure *str,TS ts)
294: {
295: Vec Xdot;
299: TSPseudoGetXdot(ts,X,&Xdot);
300: TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,str,PETSC_FALSE);
301: return(0);
302: }
307: static PetscErrorCode TSSetUp_Pseudo(TS ts)
308: {
309: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
313: VecDuplicate(ts->vec_sol,&pseudo->update);
314: VecDuplicate(ts->vec_sol,&pseudo->func);
315: VecDuplicate(ts->vec_sol,&pseudo->xdot);
316: return(0);
317: }
318: /*------------------------------------------------------------*/
322: PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
323: {
324: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
325: PetscErrorCode ierr;
326: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ts)->comm);
329: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
330: VecZeroEntries(pseudo->xdot);
331: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
332: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
333: }
334: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
335: PetscViewerASCIIPrintf(viewer,"TS %D dt %G time %G fnorm %G\n",step,ts->time_step,ptime,pseudo->fnorm);
336: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
337: return(0);
338: }
342: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts)
343: {
344: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
345: PetscErrorCode ierr;
346: PetscBool flg = PETSC_FALSE;
347: PetscViewer viewer;
350: PetscOptionsHead("Pseudo-timestepping options");
351: PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","TSPseudoMonitorDefault",flg,&flg,PETSC_NULL);
352: if (flg) {
353: PetscViewerASCIIOpen(((PetscObject)ts)->comm,"stdout",&viewer);
354: TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
355: }
356: flg = PETSC_FALSE;
357: PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,PETSC_NULL);
358: if (flg) {
359: TSPseudoIncrementDtFromInitialDt(ts);
360: }
361: PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,0);
362: PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,0);
364: SNESSetFromOptions(ts->snes);
365: PetscOptionsTail();
366: return(0);
367: }
371: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
372: {
376: SNESView(ts->snes,viewer);
377: return(0);
378: }
380: /* ----------------------------------------------------------------------------- */
383: /*@C
384: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
385: last timestep.
387: Logically Collective on TS
389: Input Parameters:
390: + ts - timestep context
391: . dt - user-defined function to verify timestep
392: - ctx - [optional] user-defined context for private data
393: for the timestep verification routine (may be PETSC_NULL)
395: Level: advanced
397: Calling sequence of func:
398: . func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
400: . update - latest solution vector
401: . ctx - [optional] timestep context
402: . newdt - the timestep to use for the next step
403: . flag - flag indicating whether the last time step was acceptable
405: Notes:
406: The routine set here will be called by TSPseudoVerifyTimeStep()
407: during the timestepping process.
409: .keywords: timestep, pseudo, set, verify
411: .seealso: TSPseudoDefaultVerifyTimeStep(), TSPseudoVerifyTimeStep()
412: @*/
413: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool *),void* ctx)
414: {
419: PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal *,PetscBool *),void *),(ts,dt,ctx));
420: return(0);
421: }
425: /*@
426: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
427: dt when using the TSPseudoDefaultTimeStep() routine.
429: Logically Collective on TS
431: Input Parameters:
432: + ts - the timestep context
433: - inc - the scaling factor >= 1.0
435: Options Database Key:
436: $ -ts_pseudo_increment <increment>
438: Level: advanced
440: .keywords: timestep, pseudo, set, increment
442: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
443: @*/
444: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
445: {
451: PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
452: return(0);
453: }
457: /*@
458: TSPseudoSetMaxTimeStep - Sets the maximum time step
459: when using the TSPseudoDefaultTimeStep() routine.
461: Logically Collective on TS
463: Input Parameters:
464: + ts - the timestep context
465: - maxdt - the maximum time step, use a non-positive value to deactivate
467: Options Database Key:
468: $ -ts_pseudo_max_dt <increment>
470: Level: advanced
472: .keywords: timestep, pseudo, set
474: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
475: @*/
476: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
477: {
483: PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
484: return(0);
485: }
489: /*@
490: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
491: is computed via the formula
492: $ dt = initial_dt*initial_fnorm/current_fnorm
493: rather than the default update,
494: $ dt = current_dt*previous_fnorm/current_fnorm.
496: Logically Collective on TS
498: Input Parameter:
499: . ts - the timestep context
501: Options Database Key:
502: $ -ts_pseudo_increment_dt_from_initial_dt
504: Level: advanced
506: .keywords: timestep, pseudo, set, increment
508: .seealso: TSPseudoSetTimeStep(), TSPseudoDefaultTimeStep()
509: @*/
510: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
511: {
516: PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
517: return(0);
518: }
523: /*@C
524: TSPseudoSetTimeStep - Sets the user-defined routine to be
525: called at each pseudo-timestep to update the timestep.
527: Logically Collective on TS
529: Input Parameters:
530: + ts - timestep context
531: . dt - function to compute timestep
532: - ctx - [optional] user-defined context for private data
533: required by the function (may be PETSC_NULL)
535: Level: intermediate
537: Calling sequence of func:
538: . func (TS ts,PetscReal *newdt,void *ctx);
540: . newdt - the newly computed timestep
541: . ctx - [optional] timestep context
543: Notes:
544: The routine set here will be called by TSPseudoComputeTimeStep()
545: during the timestepping process.
547: .keywords: timestep, pseudo, set
549: .seealso: TSPseudoDefaultTimeStep(), TSPseudoComputeTimeStep()
550: @*/
551: PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void* ctx)
552: {
557: PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal *,void *),void *),(ts,dt,ctx));
558: return(0);
559: }
561: /* ----------------------------------------------------------------------------- */
563: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool *); /* force argument to next function to not be extern C*/
564: EXTERN_C_BEGIN
567: PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void* ctx)
568: {
569: TS_Pseudo *pseudo;
572: pseudo = (TS_Pseudo*)ts->data;
573: pseudo->verify = dt;
574: pseudo->verifyctx = ctx;
575: return(0);
576: }
577: EXTERN_C_END
579: EXTERN_C_BEGIN
582: PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
583: {
584: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
587: pseudo->dt_increment = inc;
588: return(0);
589: }
590: EXTERN_C_END
592: EXTERN_C_BEGIN
595: PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
596: {
597: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
600: pseudo->dt_max = maxdt;
601: return(0);
602: }
603: EXTERN_C_END
605: EXTERN_C_BEGIN
608: PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
609: {
610: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
613: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
614: return(0);
615: }
616: EXTERN_C_END
618: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
619: EXTERN_C_BEGIN
622: PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void* ctx)
623: {
624: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
627: pseudo->dt = dt;
628: pseudo->dtctx = ctx;
629: return(0);
630: }
631: EXTERN_C_END
633: /* ----------------------------------------------------------------------------- */
634: /*MC
635: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
637: This method solves equations of the form
639: $ F(X,Xdot) = 0
641: for steady state using the iteration
643: $ [G'] S = -F(X,0)
644: $ X += S
646: where
648: $ G(Y) = F(Y,(Y-X)/dt)
650: This is linearly-implicit Euler with the residual always evaluated "at steady
651: state". See note below.
653: Options database keys:
654: + -ts_pseudo_increment <real> - ratio of increase dt
655: - -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
657: Level: beginner
659: References:
660: Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential-Algebraic Equations, 2003.
661: C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
663: Notes:
664: The residual computed by this method includes the transient term (Xdot is computed instead of
665: always being zero), but since the prediction from the last step is always the solution from the
666: last step, on the first Newton iteration we have
668: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
670: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
671: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
672: algorithm is no longer the one described in the referenced papers.
674: .seealso: TSCreate(), TS, TSSetType()
676: M*/
677: EXTERN_C_BEGIN
680: PetscErrorCode TSCreate_Pseudo(TS ts)
681: {
682: TS_Pseudo *pseudo;
684: SNES snes;
685: const SNESType stype;
688: ts->ops->reset = TSReset_Pseudo;
689: ts->ops->destroy = TSDestroy_Pseudo;
690: ts->ops->view = TSView_Pseudo;
692: ts->ops->setup = TSSetUp_Pseudo;
693: ts->ops->step = TSStep_Pseudo;
694: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
695: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
696: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
698: TSGetSNES(ts,&snes);
699: SNESGetType(snes,&stype);
700: if (!stype) {SNESSetType(snes,SNESKSPONLY);}
702: PetscNewLog(ts,TS_Pseudo,&pseudo);
703: ts->data = (void*)pseudo;
705: pseudo->dt_increment = 1.1;
706: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
707: pseudo->dt = TSPseudoDefaultTimeStep;
708: pseudo->fnorm = -1;
710: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",
711: "TSPseudoSetVerifyTimeStep_Pseudo",
712: TSPseudoSetVerifyTimeStep_Pseudo);
713: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",
714: "TSPseudoSetTimeStepIncrement_Pseudo",
715: TSPseudoSetTimeStepIncrement_Pseudo);
716: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",
717: "TSPseudoSetMaxTimeStep_Pseudo",
718: TSPseudoSetMaxTimeStep_Pseudo);
719: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",
720: "TSPseudoIncrementDtFromInitialDt_Pseudo",
721: TSPseudoIncrementDtFromInitialDt_Pseudo);
722: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSPseudoSetTimeStep_C",
723: "TSPseudoSetTimeStep_Pseudo",
724: TSPseudoSetTimeStep_Pseudo);
725: return(0);
726: }
727: EXTERN_C_END
731: /*@C
732: TSPseudoDefaultTimeStep - Default code to compute pseudo-timestepping.
733: Use with TSPseudoSetTimeStep().
735: Collective on TS
737: Input Parameters:
738: . ts - the timestep context
739: . dtctx - unused timestep context
741: Output Parameter:
742: . newdt - the timestep to use for the next step
744: Level: advanced
746: .keywords: timestep, pseudo, default
748: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
749: @*/
750: PetscErrorCode TSPseudoDefaultTimeStep(TS ts,PetscReal* newdt,void* dtctx)
751: {
752: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
753: PetscReal inc = pseudo->dt_increment,fnorm_previous = pseudo->fnorm_previous;
757: VecZeroEntries(pseudo->xdot);
758: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
759: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
760: if (pseudo->fnorm_initial == 0.0) {
761: /* first time through so compute initial function norm */
762: pseudo->fnorm_initial = pseudo->fnorm;
763: fnorm_previous = pseudo->fnorm;
764: }
765: if (pseudo->fnorm == 0.0) {
766: *newdt = 1.e12*inc*ts->time_step;
767: } else if (pseudo->increment_dt_from_initial_dt) {
768: *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
769: } else {
770: *newdt = inc*ts->time_step*fnorm_previous/pseudo->fnorm;
771: }
772: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
773: pseudo->fnorm_previous = pseudo->fnorm;
774: return(0);
775: }