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