Actual source code: posindep.c
petsc-3.9.4 2018-09-11
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: .keywords: timestep, pseudo, compute
50: .seealso: TSPseudoTimeStepDefault(), TSPseudoSetTimeStep()
51: @*/
52: PetscErrorCode TSPseudoComputeTimeStep(TS ts,PetscReal *dt)
53: {
54: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
58: PetscLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0);
59: (*pseudo->dt)(ts,dt,pseudo->dtctx);
60: PetscLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0);
61: return(0);
62: }
65: /* ------------------------------------------------------------------------------*/
66: /*@C
67: TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
69: Collective on TS
71: Input Parameters:
72: + ts - the timestep context
73: . dtctx - unused timestep context
74: - update - latest solution vector
76: Output Parameters:
77: + newdt - the timestep to use for the next step
78: - flag - flag indicating whether the last time step was acceptable
80: Level: advanced
82: Note:
83: This routine always returns a flag of 1, indicating an acceptable
84: timestep.
86: .keywords: timestep, pseudo, default, verify
88: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStep()
89: @*/
90: PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts,Vec update,void *dtctx,PetscReal *newdt,PetscBool *flag)
91: {
93: *flag = PETSC_TRUE;
94: return(0);
95: }
98: /*@
99: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
101: Collective on TS
103: Input Parameters:
104: + ts - timestep context
105: - update - latest solution vector
107: Output Parameters:
108: + dt - newly computed timestep (if it had to shrink)
109: - flag - indicates if current timestep was ok
111: Level: advanced
113: Notes:
114: The routine to be called here to compute the timestep should be
115: set by calling TSPseudoSetVerifyTimeStep().
117: .keywords: timestep, pseudo, verify
119: .seealso: TSPseudoSetVerifyTimeStep(), TSPseudoVerifyTimeStepDefault()
120: @*/
121: PetscErrorCode TSPseudoVerifyTimeStep(TS ts,Vec update,PetscReal *dt,PetscBool *flag)
122: {
123: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
127: *flag = PETSC_TRUE;
128: if(pseudo->verify) {
129: (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag);
130: }
131: return(0);
132: }
134: /* --------------------------------------------------------------------------------*/
136: static PetscErrorCode TSStep_Pseudo(TS ts)
137: {
138: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
139: PetscInt nits,lits,reject;
140: PetscBool stepok;
141: PetscReal next_time_step = ts->time_step;
142: PetscErrorCode ierr;
145: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
146: VecCopy(ts->vec_sol,pseudo->update);
147: TSPseudoComputeTimeStep(ts,&next_time_step);
148: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
149: ts->time_step = next_time_step;
150: TSPreStage(ts,ts->ptime+ts->time_step);
151: SNESSolve(ts->snes,NULL,pseudo->update);
152: SNESGetIterationNumber(ts->snes,&nits);
153: SNESGetLinearSolveIterations(ts->snes,&lits);
154: ts->snes_its += nits; ts->ksp_its += lits;
155: TSPostStage(ts,ts->ptime+ts->time_step,0,&(pseudo->update));
156: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,pseudo->update,&stepok);
157: if (!stepok) {next_time_step = ts->time_step; continue;}
158: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
159: TSPseudoVerifyTimeStep(ts,pseudo->update,&next_time_step,&stepok);
160: if (stepok) break;
161: }
162: if (reject >= ts->max_reject) {
163: ts->reason = TS_DIVERGED_STEP_REJECTED;
164: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
165: return(0);
166: }
168: VecCopy(pseudo->update,ts->vec_sol);
169: ts->ptime += ts->time_step;
170: ts->time_step = next_time_step;
172: if (pseudo->fnorm < 0) {
173: VecZeroEntries(pseudo->xdot);
174: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
175: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
176: }
177: if (pseudo->fnorm < pseudo->fatol) {
178: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
179: PetscInfo3(ts,"Step=%D, converged since fnorm %g < fatol %g\n",ts->steps,pseudo->fnorm,pseudo->frtol);
180: return(0);
181: }
182: if (pseudo->fnorm/pseudo->fnorm_initial < pseudo->frtol) {
183: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
184: PetscInfo4(ts,"Step=%D, converged since fnorm %g / fnorm_initial %g < frtol %g\n",ts->steps,pseudo->fnorm,pseudo->fnorm_initial,pseudo->fatol);
185: return(0);
186: }
187: return(0);
188: }
190: /*------------------------------------------------------------*/
191: static PetscErrorCode TSReset_Pseudo(TS ts)
192: {
193: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
197: VecDestroy(&pseudo->update);
198: VecDestroy(&pseudo->func);
199: VecDestroy(&pseudo->xdot);
200: return(0);
201: }
203: static PetscErrorCode TSDestroy_Pseudo(TS ts)
204: {
208: TSReset_Pseudo(ts);
209: PetscFree(ts->data);
210: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",NULL);
211: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",NULL);
212: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",NULL);
213: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",NULL);
214: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",NULL);
215: return(0);
216: }
218: /*------------------------------------------------------------*/
220: /*
221: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
222: */
223: static PetscErrorCode TSPseudoGetXdot(TS ts,Vec X,Vec *Xdot)
224: {
225: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
226: const PetscScalar mdt = 1.0/ts->time_step,*xnp1,*xn;
227: PetscScalar *xdot;
228: PetscErrorCode ierr;
229: PetscInt i,n;
232: *Xdot = NULL;
233: VecGetArrayRead(ts->vec_sol,&xn);
234: VecGetArrayRead(X,&xnp1);
235: VecGetArray(pseudo->xdot,&xdot);
236: VecGetLocalSize(X,&n);
237: for (i=0; i<n; i++) xdot[i] = mdt*(xnp1[i] - xn[i]);
238: VecRestoreArrayRead(ts->vec_sol,&xn);
239: VecRestoreArrayRead(X,&xnp1);
240: VecRestoreArray(pseudo->xdot,&xdot);
241: *Xdot = pseudo->xdot;
242: return(0);
243: }
245: /*
246: The transient residual is
248: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
250: or for ODE,
252: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
254: This is the function that must be evaluated for transient simulation and for
255: finite difference Jacobians. On the first Newton step, this algorithm uses
256: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
257: residual is actually the steady state residual. Pseudotransient
258: continuation as described in the literature is a linearly implicit
259: algorithm, it only takes this one Newton step with the steady state
260: residual, and then advances to the next time step.
261: */
262: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes,Vec X,Vec Y,TS ts)
263: {
264: Vec Xdot;
268: TSPseudoGetXdot(ts,X,&Xdot);
269: TSComputeIFunction(ts,ts->ptime+ts->time_step,X,Xdot,Y,PETSC_FALSE);
270: return(0);
271: }
273: /*
274: This constructs the Jacobian needed for SNES. For DAE, this is
276: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
278: and for ODE:
280: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
281: */
282: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes,Vec X,Mat AA,Mat BB,TS ts)
283: {
284: Vec Xdot;
288: TSPseudoGetXdot(ts,X,&Xdot);
289: TSComputeIJacobian(ts,ts->ptime+ts->time_step,X,Xdot,1./ts->time_step,AA,BB,PETSC_FALSE);
290: return(0);
291: }
294: static PetscErrorCode TSSetUp_Pseudo(TS ts)
295: {
296: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
300: VecDuplicate(ts->vec_sol,&pseudo->update);
301: VecDuplicate(ts->vec_sol,&pseudo->func);
302: VecDuplicate(ts->vec_sol,&pseudo->xdot);
303: return(0);
304: }
305: /*------------------------------------------------------------*/
307: static PetscErrorCode TSPseudoMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy)
308: {
309: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
311: PetscViewer viewer = (PetscViewer) dummy;
314: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
315: VecZeroEntries(pseudo->xdot);
316: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
317: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
318: }
319: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
320: PetscViewerASCIIPrintf(viewer,"TS %D dt %g time %g fnorm %g\n",step,(double)ts->time_step,(double)ptime,(double)pseudo->fnorm);
321: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
322: return(0);
323: }
325: static PetscErrorCode TSSetFromOptions_Pseudo(PetscOptionItems *PetscOptionsObject,TS ts)
326: {
327: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
329: PetscBool flg = PETSC_FALSE;
330: PetscViewer viewer;
333: PetscOptionsHead(PetscOptionsObject,"Pseudo-timestepping options");
334: PetscOptionsBool("-ts_monitor_pseudo","Monitor convergence","",flg,&flg,NULL);
335: if (flg) {
336: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),"stdout",&viewer);
337: TSMonitorSet(ts,TSPseudoMonitorDefault,viewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
338: }
339: flg = pseudo->increment_dt_from_initial_dt;
340: PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt","Increase dt as a ratio from original dt","TSPseudoIncrementDtFromInitialDt",flg,&flg,NULL);
341: pseudo->increment_dt_from_initial_dt = flg;
342: PetscOptionsReal("-ts_pseudo_increment","Ratio to increase dt","TSPseudoSetTimeStepIncrement",pseudo->dt_increment,&pseudo->dt_increment,NULL);
343: PetscOptionsReal("-ts_pseudo_max_dt","Maximum value for dt","TSPseudoSetMaxTimeStep",pseudo->dt_max,&pseudo->dt_max,NULL);
344: PetscOptionsReal("-ts_pseudo_fatol","Tolerance for norm of function","",pseudo->fatol,&pseudo->fatol,NULL);
345: PetscOptionsReal("-ts_pseudo_frtol","Relative tolerance for norm of function","",pseudo->frtol,&pseudo->frtol,NULL);
346: PetscOptionsTail();
347: return(0);
348: }
350: static PetscErrorCode TSView_Pseudo(TS ts,PetscViewer viewer)
351: {
353: PetscBool isascii;
356: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
357: if (isascii) {
358: TS_Pseudo *pseudo = (TS_Pseudo*) ts->data;
359: PetscViewerASCIIPrintf(viewer," frtol - relative tolerance in function value %g\n",(double)pseudo->frtol);
360: PetscViewerASCIIPrintf(viewer," fatol - absolute tolerance in function value %g\n",(double)pseudo->fatol);
361: PetscViewerASCIIPrintf(viewer," dt_initial - initial timestep %g\n",(double)pseudo->dt_initial);
362: PetscViewerASCIIPrintf(viewer," dt_increment - increase in timestep on successful step %g\n",(double)pseudo->dt_increment);
363: PetscViewerASCIIPrintf(viewer," dt_max - maximum time %g\n",(double)pseudo->dt_max);
364: }
365: return(0);
366: }
368: /* ----------------------------------------------------------------------------- */
369: /*@C
370: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
371: last timestep.
373: Logically Collective on TS
375: Input Parameters:
376: + ts - timestep context
377: . dt - user-defined function to verify timestep
378: - ctx - [optional] user-defined context for private data
379: for the timestep verification routine (may be NULL)
381: Level: advanced
383: Calling sequence of func:
384: . func (TS ts,Vec update,void *ctx,PetscReal *newdt,PetscBool *flag);
386: . update - latest solution vector
387: . ctx - [optional] timestep context
388: . newdt - the timestep to use for the next step
389: . flag - flag indicating whether the last time step was acceptable
391: Notes:
392: The routine set here will be called by TSPseudoVerifyTimeStep()
393: during the timestepping process.
395: .keywords: timestep, pseudo, set, verify
397: .seealso: TSPseudoVerifyTimeStepDefault(), TSPseudoVerifyTimeStep()
398: @*/
399: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts,PetscErrorCode (*dt)(TS,Vec,void*,PetscReal*,PetscBool*),void *ctx)
400: {
405: PetscTryMethod(ts,"TSPseudoSetVerifyTimeStep_C",(TS,PetscErrorCode (*)(TS,Vec,void*,PetscReal*,PetscBool*),void*),(ts,dt,ctx));
406: return(0);
407: }
409: /*@
410: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
411: dt when using the TSPseudoTimeStepDefault() routine.
413: Logically Collective on TS
415: Input Parameters:
416: + ts - the timestep context
417: - inc - the scaling factor >= 1.0
419: Options Database Key:
420: . -ts_pseudo_increment <increment>
422: Level: advanced
424: .keywords: timestep, pseudo, set, increment
426: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
427: @*/
428: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts,PetscReal inc)
429: {
435: PetscTryMethod(ts,"TSPseudoSetTimeStepIncrement_C",(TS,PetscReal),(ts,inc));
436: return(0);
437: }
439: /*@
440: TSPseudoSetMaxTimeStep - Sets the maximum time step
441: when using the TSPseudoTimeStepDefault() routine.
443: Logically Collective on TS
445: Input Parameters:
446: + ts - the timestep context
447: - maxdt - the maximum time step, use a non-positive value to deactivate
449: Options Database Key:
450: . -ts_pseudo_max_dt <increment>
452: Level: advanced
454: .keywords: timestep, pseudo, set
456: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
457: @*/
458: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts,PetscReal maxdt)
459: {
465: PetscTryMethod(ts,"TSPseudoSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
466: return(0);
467: }
469: /*@
470: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
471: is computed via the formula
472: $ dt = initial_dt*initial_fnorm/current_fnorm
473: rather than the default update,
474: $ dt = current_dt*previous_fnorm/current_fnorm.
476: Logically Collective on TS
478: Input Parameter:
479: . ts - the timestep context
481: Options Database Key:
482: . -ts_pseudo_increment_dt_from_initial_dt
484: Level: advanced
486: .keywords: timestep, pseudo, set, increment
488: .seealso: TSPseudoSetTimeStep(), TSPseudoTimeStepDefault()
489: @*/
490: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
491: {
496: PetscTryMethod(ts,"TSPseudoIncrementDtFromInitialDt_C",(TS),(ts));
497: return(0);
498: }
501: /*@C
502: TSPseudoSetTimeStep - Sets the user-defined routine to be
503: called at each pseudo-timestep to update the timestep.
505: Logically Collective on TS
507: Input Parameters:
508: + ts - timestep context
509: . dt - function to compute timestep
510: - ctx - [optional] user-defined context for private data
511: required by the function (may be NULL)
513: Level: intermediate
515: Calling sequence of func:
516: . func (TS ts,PetscReal *newdt,void *ctx);
518: . newdt - the newly computed timestep
519: . ctx - [optional] timestep context
521: Notes:
522: The routine set here will be called by TSPseudoComputeTimeStep()
523: during the timestepping process.
524: If not set then TSPseudoTimeStepDefault() is automatically used
526: .keywords: timestep, pseudo, set
528: .seealso: TSPseudoTimeStepDefault(), TSPseudoComputeTimeStep()
529: @*/
530: PetscErrorCode TSPseudoSetTimeStep(TS ts,PetscErrorCode (*dt)(TS,PetscReal*,void*),void *ctx)
531: {
536: PetscTryMethod(ts,"TSPseudoSetTimeStep_C",(TS,PetscErrorCode (*)(TS,PetscReal*,void*),void*),(ts,dt,ctx));
537: return(0);
538: }
540: /* ----------------------------------------------------------------------------- */
542: typedef PetscErrorCode (*FCN1)(TS,Vec,void*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/
543: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts,FCN1 dt,void *ctx)
544: {
545: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
548: pseudo->verify = dt;
549: pseudo->verifyctx = ctx;
550: return(0);
551: }
553: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts,PetscReal inc)
554: {
555: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
558: pseudo->dt_increment = inc;
559: return(0);
560: }
562: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts,PetscReal maxdt)
563: {
564: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
567: pseudo->dt_max = maxdt;
568: return(0);
569: }
571: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
572: {
573: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
576: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
577: return(0);
578: }
580: typedef PetscErrorCode (*FCN2)(TS,PetscReal*,void*); /* force argument to next function to not be extern C*/
581: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts,FCN2 dt,void *ctx)
582: {
583: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
586: pseudo->dt = dt;
587: pseudo->dtctx = ctx;
588: return(0);
589: }
591: /* ----------------------------------------------------------------------------- */
592: /*MC
593: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
595: This method solves equations of the form
597: $ F(X,Xdot) = 0
599: for steady state using the iteration
601: $ [G'] S = -F(X,0)
602: $ X += S
604: where
606: $ G(Y) = F(Y,(Y-X)/dt)
608: This is linearly-implicit Euler with the residual always evaluated "at steady
609: state". See note below.
611: Options database keys:
612: + -ts_pseudo_increment <real> - ratio of increase dt
613: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
614: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
615: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
617: Level: beginner
619: References:
620: + 1. - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
621: - 2. - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
623: Notes:
624: The residual computed by this method includes the transient term (Xdot is computed instead of
625: always being zero), but since the prediction from the last step is always the solution from the
626: last step, on the first Newton iteration we have
628: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
630: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
631: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
632: algorithm is no longer the one described in the referenced papers.
634: .seealso: TSCreate(), TS, TSSetType()
636: M*/
637: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
638: {
639: TS_Pseudo *pseudo;
641: SNES snes;
642: SNESType stype;
645: ts->ops->reset = TSReset_Pseudo;
646: ts->ops->destroy = TSDestroy_Pseudo;
647: ts->ops->view = TSView_Pseudo;
648: ts->ops->setup = TSSetUp_Pseudo;
649: ts->ops->step = TSStep_Pseudo;
650: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
651: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
652: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
653: ts->default_adapt_type = TSADAPTNONE;
654: ts->usessnes = PETSC_TRUE;
656: TSGetSNES(ts,&snes);
657: SNESGetType(snes,&stype);
658: if (!stype) {SNESSetType(snes,SNESKSPONLY);}
660: PetscNewLog(ts,&pseudo);
661: ts->data = (void*)pseudo;
663: pseudo->dt = TSPseudoTimeStepDefault;
664: pseudo->dtctx = NULL;
665: pseudo->dt_increment = 1.1;
666: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
667: pseudo->fnorm = -1;
668: pseudo->fnorm_initial = -1;
669: pseudo->fnorm_previous = -1;
670: #if defined(PETSC_USE_REAL_SINGLE)
671: pseudo->fatol = 1.e-25;
672: pseudo->frtol = 1.e-5;
673: #else
674: pseudo->fatol = 1.e-50;
675: pseudo->frtol = 1.e-12;
676: #endif
677: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetVerifyTimeStep_C",TSPseudoSetVerifyTimeStep_Pseudo);
678: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStepIncrement_C",TSPseudoSetTimeStepIncrement_Pseudo);
679: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetMaxTimeStep_C",TSPseudoSetMaxTimeStep_Pseudo);
680: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoIncrementDtFromInitialDt_C",TSPseudoIncrementDtFromInitialDt_Pseudo);
681: PetscObjectComposeFunction((PetscObject)ts,"TSPseudoSetTimeStep_C",TSPseudoSetTimeStep_Pseudo);
682: return(0);
683: }
685: /*@C
686: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping.
687: Use with TSPseudoSetTimeStep().
689: Collective on TS
691: Input Parameters:
692: . ts - the timestep context
693: . dtctx - unused timestep context
695: Output Parameter:
696: . newdt - the timestep to use for the next step
698: Level: advanced
700: .keywords: timestep, pseudo, default
702: .seealso: TSPseudoSetTimeStep(), TSPseudoComputeTimeStep()
703: @*/
704: PetscErrorCode TSPseudoTimeStepDefault(TS ts,PetscReal *newdt,void *dtctx)
705: {
706: TS_Pseudo *pseudo = (TS_Pseudo*)ts->data;
707: PetscReal inc = pseudo->dt_increment;
711: VecZeroEntries(pseudo->xdot);
712: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,pseudo->xdot,pseudo->func,PETSC_FALSE);
713: VecNorm(pseudo->func,NORM_2,&pseudo->fnorm);
714: if (pseudo->fnorm_initial < 0) {
715: /* first time through so compute initial function norm */
716: pseudo->fnorm_initial = pseudo->fnorm;
717: pseudo->fnorm_previous = pseudo->fnorm;
718: }
719: if (pseudo->fnorm == 0.0) *newdt = 1.e12*inc*ts->time_step;
720: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc*pseudo->dt_initial*pseudo->fnorm_initial/pseudo->fnorm;
721: else *newdt = inc*ts->time_step*pseudo->fnorm_previous/pseudo->fnorm;
722: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt,pseudo->dt_max);
723: pseudo->fnorm_previous = pseudo->fnorm;
724: return(0);
725: }