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
36: Input Parameter:
37: . ts - timestep context
39: Output Parameter:
40: . dt - newly computed timestep
42: Level: developer
44: Note:
45: The routine to be called here to compute the timestep should be
46: set by calling `TSPseudoSetTimeStep()`.
48: .seealso: [](chapter_ts), `TSPSEUDO`, `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
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: [](chapter_ts), `TSPSEUDO`, `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
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: [](chapter_ts), `TSPSEUDO`, `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
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: 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: Level: advanced
358: Note:
359: The routine set here will be called by `TSPseudoVerifyTimeStep()`
360: during the timestepping process.
362: .seealso: [](chapter_ts), `TSPSEUDO`, `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
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: [](chapter_ts), `TSPSEUDO`, `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
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: [](chapter_ts), `TSPSEUDO`, `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
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: [](chapter_ts), `TSPSEUDO`, `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
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: Calling sequence of func:
460: $ func (TS ts,PetscReal *newdt,void *ctx);
462: + newdt - the newly computed timestep
463: - ctx - [optional] timestep context
465: Level: intermediate
467: Notes:
468: The routine set here will be called by `TSPseudoComputeTimeStep()`
469: during the timestepping process.
471: If not set then `TSPseudoTimeStepDefault()` is automatically used
473: .seealso: [](chapter_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
474: @*/
475: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS, PetscReal *, void *), void *ctx)
476: {
478: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode(*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
479: return 0;
480: }
482: /* ----------------------------------------------------------------------------- */
484: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
485: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, void *ctx)
486: {
487: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
489: pseudo->verify = dt;
490: pseudo->verifyctx = ctx;
491: return 0;
492: }
494: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
495: {
496: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
498: pseudo->dt_increment = inc;
499: return 0;
500: }
502: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
503: {
504: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
506: pseudo->dt_max = maxdt;
507: return 0;
508: }
510: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
511: {
512: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
514: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
515: return 0;
516: }
518: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
519: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, void *ctx)
520: {
521: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
523: pseudo->dt = dt;
524: pseudo->dtctx = ctx;
525: return 0;
526: }
528: /* ----------------------------------------------------------------------------- */
529: /*MC
530: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping
532: This method solves equations of the form
534: $ F(X,Xdot) = 0
536: for steady state using the iteration
538: $ [G'] S = -F(X,0)
539: $ X += S
541: where
543: $ G(Y) = F(Y,(Y-X)/dt)
545: This is linearly-implicit Euler with the residual always evaluated "at steady
546: state". See note below.
548: Options Database Keys:
549: + -ts_pseudo_increment <real> - ratio of increase dt
550: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
551: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
552: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
554: Level: beginner
556: Notes:
557: The residual computed by this method includes the transient term (Xdot is computed instead of
558: always being zero), but since the prediction from the last step is always the solution from the
559: last step, on the first Newton iteration we have
561: $ Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
563: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
564: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
565: algorithm is no longer the one described in the referenced papers.
567: References:
568: + * - Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient Continuation and Differential Algebraic Equations, 2003.
569: - * - C. T. Kelley and David E. Keyes, Convergence analysis of Pseudotransient Continuation, 1998.
571: .seealso: [](chapter_ts), `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. Use with `TSPseudoSetTimeStep()`.
622: Collective
624: Input Parameters:
625: + ts - the timestep context
626: - dtctx - unused timestep context
628: Output Parameter:
629: . newdt - the timestep to use for the next step
631: Level: advanced
633: .seealso: [](chapter_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
634: @*/
635: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
636: {
637: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
638: PetscReal inc = pseudo->dt_increment;
640: VecZeroEntries(pseudo->xdot);
641: TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE);
642: VecNorm(pseudo->func, NORM_2, &pseudo->fnorm);
643: if (pseudo->fnorm_initial < 0) {
644: /* first time through so compute initial function norm */
645: pseudo->fnorm_initial = pseudo->fnorm;
646: pseudo->fnorm_previous = pseudo->fnorm;
647: }
648: if (pseudo->fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
649: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / pseudo->fnorm;
650: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / pseudo->fnorm;
651: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
652: pseudo->fnorm_previous = pseudo->fnorm;
653: return 0;
654: }