Actual source code: posindep.c
petsc-3.7.5 2017-01-01
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: PetscReal fatol,frtol;
26: } TS_Pseudo;
28: /* ------------------------------------------------------------------------------*/
32: /*@C
33: TSPseudoComputeTimeStep - Computes the next timestep for a currently running
34: pseudo-timestepping process.
36: Collective on TS
38: Input Parameter:
39: . ts - timestep context
41: Output Parameter:
42: . dt - newly computed timestep
44: Level: developer
46: Notes:
47: The routine to be called here to compute the timestep should be
48: set by calling TSPseudoSetTimeStep().
50: .keywords: timestep, pseudo, compute
52: .seealso: TSPseudoTimeStepDefault(), TSPseudoSetTimeStep()
53: @*/
54: PetscErrorCode TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
55: {
56: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
60: PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);
61: (*pseudo->dt)(ts,dt,pseudo->dtctx);
62: PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);
63: return(0);
64: }
67: /* ------------------------------------------------------------------------------*/
70: /*@C
71: TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
73: Collective on TS
75: Input Parameters:
76: + ts - the timestep context
77: . dtctx - unused timestep context
78: - update - latest solution vector
80: Output Parameters:
81: + newdt - the timestep to use for the next step
82: - flag - flag indicating whether the last time step was acceptable
84: Level: advanced
86: Note:
87: This routine always returns a flag of 1, indicating an acceptable
88: timestep.
90: .keywords: timestep, pseudo, default, verify
92: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
93: @*/
94: PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool *flag)
95: {
97: *flag = PETSC_TRUE;
98: return(0);
99: }
104: /*@
105: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
107: Collective on TS
109: Input Parameters:
110: + ts - timestep context
111: - update - latest solution vector
113: Output Parameters:
114: + dt - newly computed timestep (if it had to shrink)
115: - flag - indicates if current timestep was ok
117: Level: advanced
119: Notes:
120: The routine to be called here to compute the timestep should be
121: set by calling TSPseudoSetVerifyTimeStep().
123: .keywords: timestep, pseudo, verify
125: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
126: @*/
127: PetscErrorCode TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
128: {
129: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
133: *flag = PETSC_TRUE;
134: if(pseudo->verify) {
135: (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
136: }
137: return(0);
138: }
140: /* --------------------------------------------------------------------------------*/
144: static PetscErrorCode TSStep_Pseudo(TS ts)
145: {
146: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
147: PetscInt nits,lits,reject;
148: PetscBool stepok;
149: PetscReal next_time_step = ts->time_step;
150: PetscErrorCode ierr;
153: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
154: VecCopy(ts->vec_sol,pseudo->update);
155: TSPseudoComputeTimeStep(ts,&next_time_step);
156: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
157: ts->time_step = next_time_step;
158: TSPreStage(ts,ts->ptime+ts->time_step);
159: SNESSolve(ts->snes,NULL,pseudo->update);
160: SNESGetIterationNumber(ts->snes,&nits);
161: SNESGetLinearSolveIterations(ts->snes,&lits);
162: ts->snes_its += nits; ts->ksp_its += lits;
163: TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));
164: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);
165: if (!stepok) {next_time_step = ts->time_step; continue;}
166: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
167: TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
168: if (stepok) break;
169: }
170: if (reject >= ts->max_reject) {
171: ts->reason = TS_DIVERGED_STEP_REJECTED;
172: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
173: return(0);
174: }
176: VecCopy(pseudo->update,ts->vec_sol);
177: ts->ptime += ts->time_step;
178: ts->time_step = next_time_step;
180: if (pseudo->fnorm < 0) {
181: VecZeroEntries(pseudo->xdot);
182: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
183: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
184: }
185: if (pseudo->fnorm < pseudo->fatol) {
186: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
187: PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
188: return(0);
189: }
190: if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
191: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
192: PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
193: return(0);
194: }
195: return(0);
196: }
198: /*------------------------------------------------------------*/
201: static PetscErrorCode TSReset_Pseudo(TS ts)
202: {
203: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
207: VecDestroy(&pseudo->update);
208: VecDestroy(&pseudo->func);
209: VecDestroy(&pseudo->xdot);
210: return(0);
211: }
215: static PetscErrorCode TSDestroy_Pseudo(TS ts)
216: {
220: TSReset_Pseudo(ts);
221: PetscFree(ts->data);
222: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
223: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
224: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
225: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
226: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
227: return(0);
228: }
230: /*------------------------------------------------------------*/
234: /*
235: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
236: */
237: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
238: {
239: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
240: const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
241: PetscScalar *xdot;
242: PetscErrorCode ierr;
243: PetscInt i,n;
246: VecGetArrayRead(ts->vec_sol,&xn);
247: VecGetArrayRead(X,&xnp1);
248: VecGetArray(pseudo->xdot,&xdot);
249: VecGetLocalSize(X,&n);
250: for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
251: VecRestoreArrayRead(ts->vec_sol,&xn);
252: VecRestoreArrayRead(X,&xnp1);
253: VecRestoreArray(pseudo->xdot,&xdot);
254: *Xdot = pseudo->xdot;
255: return(0);
256: }
260: /*
261: The transient residual is
263: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
265: or for ODE,
267: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
269: This is the function that must be evaluated for transient simulation and for
270: finite difference Jacobians. On the first Newton step, this algorithm uses
271: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
272: residual is actually the steady state residual. Pseudotransient
273: continuation as described in the literature is a linearly implicit
274: algorithm, it only takes this one Newton step with the steady state
275: residual, and then advances to the next time step.
276: */
277: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
278: {
279: Vec Xdot;
283: TSPseudoGetXdot(ts,X,&Xdot);
284: TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
285: return(0);
286: }
290: /*
291: This constructs the Jacobian needed for SNES. For DAE, this is
293: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
295: and for ODE:
297: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
298: */
299: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
300: {
301: Vec Xdot;
305: TSPseudoGetXdot(ts,X,&Xdot);
306: TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
307: return(0);
308: }
313: static PetscErrorCode TSSetUp_Pseudo(TS ts)
314: {
315: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
319: VecDuplicate(ts->vec_sol,&pseudo->update);
320: VecDuplicate(ts->vec_sol,&pseudo->func);
321: VecDuplicate(ts->vec_sol,&pseudo->xdot);
322: return(0);
323: }
324: /*------------------------------------------------------------*/
328: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
329: {
330: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
332: PetscViewer viewer = (PetscViewer) dummy;
335: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
336: VecZeroEntries(pseudo->xdot);
337: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
338: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
339: }
340: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
341: PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
342: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
343: return(0);
344: }
348: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
349: {
350: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
352: PetscBool flg = PETSC_FALSE;
353: PetscViewer viewer;
356: PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
357: PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
358: if (flg) {
359: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
360: TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
361: }
362: flg = pseudo->increment_dt_from_initial_dt;
363: PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
364: pseudo->increment_dt_from_initial_dt = flg;
365: PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
366: PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
367: PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
368: PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
369: PetscOptionsTail();
370: return(0);
371: }
375: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
376: {
378: PetscBool isascii;
381: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
382: if (isascii) {
383: TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
384: PetscViewerASCIIPrintf(viewer,"Parameters for pseudo timestepping\n");
385: PetscViewerASCIIPrintf(viewer," frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
386: PetscViewerASCIIPrintf(viewer," fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
387: PetscViewerASCIIPrintf(viewer," dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
388: PetscViewerASCIIPrintf(viewer," dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
389: PetscViewerASCIIPrintf(viewer," dt_max - maximum time %g\n",(double)pseudo->dt_max);
390: }
391: if (ts->snes) {SNESView(ts->snes,viewer);}
392: return(0);
393: }
395: /* ----------------------------------------------------------------------------- */
398: /*@C
399: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
400: last timestep.
402: Logically Collective on TS
404: Input Parameters:
405: + ts - timestep context
406: . dt - user-defined function to verify timestep
407: - ctx - [optional] user-defined context for private data
408: for the timestep verification routine (may be NULL)
410: Level: advanced
412: Calling sequence of func:
413: . func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
415: . update - latest solution vector
416: . ctx - [optional] timestep context
417: . newdt - the timestep to use for the next step
418: . flag - flag indicating whether the last time step was acceptable
420: Notes:
421: The routine set here will be called by TSPseudoVerifyTimeStep()
422: during the timestepping process.
424: .keywords: timestep, pseudo, set, verify
426: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
427: @*/
428: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
429: {
434: PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
435: return(0);
436: }
440: /*@
441: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
442: dt when using the TSPseudoTimeStepDefault() routine.
444: Logically Collective on TS
446: Input Parameters:
447: + ts - the timestep context
448: - inc - the scaling factor >= 1.0
450: Options Database Key:
451: . -ts_pseudo_increment <increment>
453: Level: advanced
455: .keywords: timestep, pseudo, set, increment
457: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
458: @*/
459: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
460: {
466: PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
467: return(0);
468: }
472: /*@
473: TSPseudoSetMaxTimeStep - Sets the maximum time step
474: when using the TSPseudoTimeStepDefault() routine.
476: Logically Collective on TS
478: Input Parameters:
479: + ts - the timestep context
480: - maxdt - the maximum time step, use a non-positive value to deactivate
482: Options Database Key:
483: . -ts_pseudo_max_dt <increment>
485: Level: advanced
487: .keywords: timestep, pseudo, set
489: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
490: @*/
491: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
492: {
498: PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
499: return(0);
500: }
504: /*@
505: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
506: is computed via the formula
507: $ dt = initial_dt*initial_fnorm/current_fnorm
508: rather than the default update,
509: $ dt = current_dt*previous_fnorm/current_fnorm.
511: Logically Collective on TS
513: Input Parameter:
514: . ts - the timestep context
516: Options Database Key:
517: . -ts_pseudo_increment_dt_from_initial_dt
519: Level: advanced
521: .keywords: timestep, pseudo, set, increment
523: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
524: @*/
525: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
526: {
531: PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
532: return(0);
533: }
538: /*@C
539: TSPseudoSetTimeStep - Sets the user-defined routine to be
540: called at each pseudo-timestep to update the timestep.
542: Logically Collective on TS
544: Input Parameters:
545: + ts - timestep context
546: . dt - function to compute timestep
547: - ctx - [optional] user-defined context for private data
548: required by the function (may be NULL)
550: Level: intermediate
552: Calling sequence of func:
553: . func (TS ts,PetscReal *newdt,void *ctx);
555: . newdt - the newly computed timestep
556: . ctx - [optional] timestep context
558: Notes:
559: The routine set here will be called by TSPseudoComputeTimeStep()
560: during the timestepping process.
561: If not set then TSPseudoTimeStepDefault() is automatically used
563: .keywords: timestep, pseudo, set
565: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
566: @*/
567: PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
568: {
573: PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
574: return(0);
575: }
577: /* ----------------------------------------------------------------------------- */
579: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/
582: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
583: {
584: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
587: pseudo->verify = dt;
588: pseudo->verifyctx = ctx;
589: return(0);
590: }
594: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
595: {
596: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
599: pseudo->dt_increment = inc;
600: return(0);
601: }
605: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
606: {
607: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
610: pseudo->dt_max = maxdt;
611: return(0);
612: }
616: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
617: {
618: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
621: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
622: return(0);
623: }
625: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
628: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
629: {
630: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
633: pseudo->dt = dt;
634: pseudo->dtctx = ctx;
635: return(0);
636: }
638: /* ----------------------------------------------------------------------------- */
639: /*MC
640: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
642: This method solves equations of the form
644: $ F(X,Xdot) = 0
646: for steady state using the iteration
648: $ [G'] S = -F(X,0)
649: $ X += S
651: where
653: $ G(Y) = F(Y,(Y-X)/dt)
655: This is linearly-implicit Euler with the residual always evaluated "at steady
656: state". See note below.
658: Options database keys:
659: + -ts_pseudo_increment <real> - ratio of increase dt
660: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
661: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
662: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
664: Level: beginner
666: References:
667: + 1. - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
668: - 2. - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
670: Notes:
671: The residual computed by this method includes the transient term (Xdot is computed instead of
672: always being zero), but since the prediction from the last step is always the solution from the
673: last step, on the first Newton iteration we have
675: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
677: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
678: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
679: algorithm is no longer the one described in the referenced papers.
681: .seealso: TSCreate(), TS, TSSetType()
683: M*/
686: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
687: {
688: TS_Pseudo *pseudo;
690: SNES snes;
691: SNESType stype;
694: ts->ops->reset = TSReset_Pseudo;
695: ts->ops->destroy = TSDestroy_Pseudo;
696: ts->ops->view = TSView_Pseudo;
697: ts->ops->setup = TSSetUp_Pseudo;
698: ts->ops->step = TSStep_Pseudo;
699: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
700: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
701: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
703: TSGetSNES(ts,&snes);
704: SNESGetType(snes,&stype);
705: if (!stype) {SNESSetType(snes,SNESKSPONLY);}
707: PetscNewLog(ts,&pseudo);
708: ts->data = (void*)pseudo;
710: pseudo->dt = TSPseudoTimeStepDefault;
711: pseudo->dtctx = NULL;
712: pseudo->dt_increment = 1.1;
713: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
714: pseudo->fnorm = -1;
715: pseudo->fnorm_initial = -1;
716: pseudo->fnorm_previous = -1;
717: #if defined(PETSC_USE_REAL_SINGLE)
718: pseudo->fatol = 1.e-25;
719: pseudo->frtol = 1.e-5;
720: #else
721: pseudo->fatol = 1.e-50;
722: pseudo->frtol = 1.e-12;
723: #endif
724: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
725: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
726: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
727: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
728: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
729: return(0);
730: }
734: /*@C
735: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
736: Use with TSPseudoSetTimeStep().
738: Collective on TS
740: Input Parameters:
741: . ts - the timestep context
742: . dtctx - unused timestep context
744: Output Parameter:
745: . newdt - the timestep to use for the next step
747: Level: advanced
749: .keywords: timestep, pseudo, default
751: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
752: @*/
753: PetscErrorCode TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
754: {
755: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
756: PetscReal inc = pseudo->dt_increment;
760: VecZeroEntries(pseudo->xdot);
761: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
762: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
763: if (pseudo->fnorm_initial < 0) {
764: /* first time through so compute initial function norm */
765: pseudo->fnorm_initial = pseudo->fnorm;
766: pseudo->fnorm_previous = pseudo->fnorm;
767: }
768: if (pseudo->fnorm == 0.0) *newdt = 1.e12*inc*ts->time_step;
769: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
770: else *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
771: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
772: pseudo->fnorm_previous = pseudo->fnorm;
773: return(0);
774: }