Actual source code: posindep.c
1: /*
2: Code for Timestepping with implicit backwards Euler.
3: */
4: #include <petsc/private/tsimpl.h>
6: #define TSADAPTTSPSEUDO "tspseudo"
8: typedef struct {
9: Vec func; /* work vector where F(t[i],u[i]) is stored */
10: Vec xdot; /* work vector for time derivative of state */
12: /* information used for Pseudo-timestepping */
14: PetscErrorCode (*dt)(TS, PetscReal *, void *); /* compute next timestep, and related context */
15: void *dtctx;
16: PetscErrorCode (*verify)(TS, Vec, void *, PetscReal *, PetscBool *); /* verify previous timestep and related context */
17: void *verifyctx;
19: PetscReal fnorm_initial, fnorm; /* original and current norm of F(u) */
20: PetscReal fnorm_previous;
22: PetscReal dt_initial; /* initial time-step */
23: PetscReal dt_increment; /* scaling that dt is incremented each time-step */
24: PetscReal dt_max; /* maximum time step */
25: PetscBool increment_dt_from_initial_dt;
26: PetscReal fatol, frtol;
28: PetscObjectState Xstate; /* state of input vector to TSComputeIFunction() with 0 Xdot */
30: TSStepStatus status;
31: } TS_Pseudo;
33: /*@C
34: TSPseudoComputeFunction - Compute nonlinear residual for pseudo-timestepping
36: This computes the residual for $\dot U = 0$, i.e. $F(U, 0)$ for the IFunction.
38: Collective
40: Input Parameters:
41: + ts - the timestep context
42: - solution - the solution vector
44: Output Parameter:
45: + residual - the nonlinear residual
46: - fnorm - the norm of the nonlinear residual
48: Level: advanced
50: Note:
51: `TSPSEUDO` records the nonlinear residual and the `solution` vector used to generate it. If given the same `solution` vector (as determined by the vectors `PetscObjectState`), this function will return those recorded values.
53: This would be used in a custom adaptive timestepping implementation that needs access to the residual, but reuses the calculation done by `TSPSEUDO` by default.
55: To correctly get the residual reuse behavior, `solution` must be the same `Vec` that returned by `TSGetSolution()`.
57: .seealso: [](ch_ts), `TSPSEUDO`
58: @*/
59: PetscErrorCode TSPseudoComputeFunction(TS ts, Vec solution, Vec *residual, PetscReal *fnorm)
60: {
61: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
62: PetscObjectState Xstate;
64: PetscFunctionBegin;
68: if (fnorm) PetscAssertPointer(fnorm, 4);
70: PetscCall(PetscObjectStateGet((PetscObject)solution, &Xstate));
71: if (Xstate != pseudo->Xstate || pseudo->fnorm < 0) {
72: PetscCall(VecZeroEntries(pseudo->xdot));
73: PetscCall(TSComputeIFunction(ts, ts->ptime, solution, pseudo->xdot, pseudo->func, PETSC_FALSE));
74: pseudo->Xstate = Xstate;
75: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
76: }
77: if (residual) *residual = pseudo->func;
78: if (fnorm) *fnorm = pseudo->fnorm;
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode TSStep_Pseudo(TS ts)
83: {
84: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
85: PetscInt nits, lits, rejections = 0;
86: PetscBool accept;
87: PetscReal next_time_step = ts->time_step, fnorm;
88: TSAdapt adapt;
90: PetscFunctionBegin;
91: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
93: pseudo->status = TS_STEP_INCOMPLETE;
94: while (!ts->reason && pseudo->status != TS_STEP_COMPLETE) {
95: PetscCall(TSPreStage(ts, ts->ptime + ts->time_step));
96: PetscCall(SNESSolve(ts->snes, NULL, ts->vec_sol));
97: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
98: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
99: ts->snes_its += nits;
100: ts->ksp_its += lits;
101: pseudo->fnorm = -1; /* The current norm is no longer valid */
102: PetscCall(TSPostStage(ts, ts->ptime + ts->time_step, 0, &ts->vec_sol));
103: PetscCall(TSGetAdapt(ts, &adapt));
104: PetscCall(TSAdaptCheckStage(adapt, ts, ts->ptime + ts->time_step, ts->vec_sol, &accept));
105: if (!accept) goto reject_step;
107: pseudo->status = TS_STEP_PENDING;
108: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
109: pseudo->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
110: if (!accept) {
111: ts->time_step = next_time_step;
112: goto reject_step;
113: }
114: ts->ptime += ts->time_step;
115: ts->time_step = next_time_step;
116: break;
118: reject_step:
119: ts->reject++;
120: accept = PETSC_FALSE;
121: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
122: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
123: ts->reason = TS_DIVERGED_STEP_REJECTED;
124: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
127: }
129: // Check solution convergence
130: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
132: if (fnorm < pseudo->fatol) {
133: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
134: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->frtol));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
137: if (fnorm / pseudo->fnorm_initial < pseudo->frtol) {
138: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
139: PetscCall(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));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode TSReset_Pseudo(TS ts)
146: {
147: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
149: PetscFunctionBegin;
150: PetscCall(VecDestroy(&pseudo->func));
151: PetscCall(VecDestroy(&pseudo->xdot));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode TSDestroy_Pseudo(TS ts)
156: {
157: PetscFunctionBegin;
158: PetscCall(TSReset_Pseudo(ts));
159: PetscCall(PetscFree(ts->data));
160: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL));
161: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL));
162: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL));
163: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL));
164: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
169: {
170: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
172: PetscFunctionBegin;
173: *Xdot = pseudo->xdot;
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*
178: The transient residual is
180: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
182: or for ODE,
184: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
186: This is the function that must be evaluated for transient simulation and for
187: finite difference Jacobians. On the first Newton step, this algorithm uses
188: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
189: residual is actually the steady state residual. Pseudotransient
190: continuation as described in the literature is a linearly implicit
191: algorithm, it only takes this one full Newton step with the steady state
192: residual, and then advances to the next time step.
194: This routine avoids a repeated identical call to TSComputeRHSFunction() when that result
195: is already contained in pseudo->func due to the call to TSComputeIFunction() in TSStep_Pseudo()
197: */
198: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
199: {
200: TSIFunctionFn *ifunction;
201: TSRHSFunctionFn *rhsfunction;
202: void *ctx;
203: DM dm;
204: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
205: const PetscScalar mdt = 1.0 / ts->time_step;
206: PetscObjectState Xstate;
207: PetscInt snes_it;
209: PetscFunctionBegin;
215: PetscCall(TSGetDM(ts, &dm));
216: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
217: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
218: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
220: PetscCall(PetscObjectStateGet((PetscObject)X, &Xstate));
221: PetscCall(SNESGetIterationNumber(snes, &snes_it));
222: if (Xstate == pseudo->Xstate && snes_it == 1) {
223: /* reuse the TSComputeIFunction() result performed inside TSStep_Pseudo() */
224: if (ifunction) PetscCall(VecCopy(pseudo->func, Y));
225: /* note that pseudo->func contains the negation of TSComputeRHSFunction() */
226: else PetscCall(VecWAXPY(Y, 1, pseudo->func, pseudo->xdot));
227: } else {
228: PetscCall(VecAXPBYPCZ(pseudo->xdot, -mdt, mdt, 0, ts->vec_sol0, X));
229: PetscCall(TSComputeIFunction(ts, ts->ptime + ts->time_step, X, pseudo->xdot, Y, PETSC_FALSE));
230: }
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*
235: This constructs the Jacobian needed for SNES. For DAE, this is
237: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
239: and for ODE:
241: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
242: */
243: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
244: {
245: Vec Xdot;
247: PetscFunctionBegin;
248: /* Xdot has already been computed in SNESTSFormFunction_Pseudo (SNES guarantees this) */
249: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
250: PetscCall(TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: static PetscErrorCode TSSetUp_Pseudo(TS ts)
255: {
256: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
258: PetscFunctionBegin;
259: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->func));
260: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->xdot));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
265: {
266: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
267: PetscViewer viewer = (PetscViewer)dummy;
269: PetscFunctionBegin;
270: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, NULL));
271: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
272: PetscCall(PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)pseudo->fnorm));
273: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems PetscOptionsObject)
278: {
279: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
280: PetscBool flg = PETSC_FALSE;
281: PetscViewer viewer;
283: PetscFunctionBegin;
284: PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
285: PetscCall(PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL));
286: if (flg) {
287: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer));
288: PetscCall(TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscCtxDestroyFn *)PetscViewerDestroy));
289: }
290: flg = pseudo->increment_dt_from_initial_dt;
291: PetscCall(PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL));
292: pseudo->increment_dt_from_initial_dt = flg;
293: PetscCall(PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL));
294: PetscCall(PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL));
295: PetscCall(PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL));
296: PetscCall(PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL));
297: PetscOptionsHeadEnd();
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
302: {
303: PetscBool isascii;
305: PetscFunctionBegin;
306: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
307: if (isascii) {
308: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
309: PetscCall(PetscViewerASCIIPrintf(viewer, " frtol - relative tolerance in function value %g\n", (double)pseudo->frtol));
310: PetscCall(PetscViewerASCIIPrintf(viewer, " fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol));
311: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_initial - initial timestep %g\n", (double)pseudo->dt_initial));
312: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment));
313: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_max - maximum time %g\n", (double)pseudo->dt_max));
314: }
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@C
319: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
321: Collective, No Fortran Support
323: Input Parameters:
324: + ts - the timestep context
325: - dtctx - unused timestep context
327: Output Parameter:
328: . newdt - the timestep to use for the next step
330: Level: advanced
332: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
333: @*/
334: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
335: {
336: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
337: PetscReal inc = pseudo->dt_increment, fnorm;
339: PetscFunctionBegin;
340: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
341: if (pseudo->fnorm_initial < 0) {
342: /* first time through so compute initial function norm */
343: pseudo->fnorm_initial = fnorm;
344: pseudo->fnorm_previous = fnorm;
345: }
346: if (fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
347: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / fnorm;
348: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / fnorm;
349: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
350: pseudo->fnorm_previous = fnorm;
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode TSAdaptChoose_TSPseudo(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
355: {
356: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
358: PetscFunctionBegin;
359: PetscCall(PetscLogEventBegin(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
360: PetscCallBack("TSPSEUDO callback time step", (*pseudo->dt)(ts, next_h, pseudo->dtctx));
361: PetscCall(PetscLogEventEnd(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
363: *next_sc = 0;
364: *wlte = -1; /* Weighted local truncation error was not evaluated */
365: *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
366: *wlter = -1; /* Weighted relative local truncation error was not evaluated */
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode TSAdaptCheckStage_TSPseudo(TSAdapt adapt, TS ts, PetscReal t, Vec Y, PetscBool *accept)
371: {
372: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
374: PetscFunctionBegin;
375: if (pseudo->verify) {
376: PetscReal dt;
377: PetscCall(TSGetTimeStep(ts, &dt));
378: PetscCallBack("TSPSEUDO callback verify time step", (*pseudo->verify)(ts, Y, pseudo->verifyctx, &dt, accept));
379: PetscCall(TSSetTimeStep(ts, dt));
380: }
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*MC
385: TSADAPTTSPSEUDO - TSPseudo adaptive controller for time stepping
387: Level: developer
389: Note:
390: This is only meant to be used with `TSPSEUDO` time integrator.
392: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`, `TSPSEUDO`
393: M*/
394: static PetscErrorCode TSAdaptCreate_TSPseudo(TSAdapt adapt)
395: {
396: PetscFunctionBegin;
397: adapt->ops->choose = TSAdaptChoose_TSPseudo;
398: adapt->checkstage = TSAdaptCheckStage_TSPseudo;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@C
403: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
404: last timestep.
406: Logically Collective
408: Input Parameters:
409: + ts - timestep context
410: . dt - user-defined function to verify timestep
411: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
413: Calling sequence of `func`:
414: + ts - the time-step context
415: . update - latest solution vector
416: . ctx - [optional] user-defined timestep context
417: . newdt - the timestep to use for the next step
418: - flag - flag indicating whether the last time step was acceptable
420: Level: advanced
422: Note:
423: The routine set here will be called by `TSPseudoVerifyTimeStep()`
424: during the timestepping process.
426: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
427: @*/
428: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, void *ctx, PetscReal *newdt, PetscBool *flag), void *ctx)
429: {
430: PetscFunctionBegin;
432: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode (*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*@
437: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
438: dt when using the TSPseudoTimeStepDefault() routine.
440: Logically Collective
442: Input Parameters:
443: + ts - the timestep context
444: - inc - the scaling factor >= 1.0
446: Options Database Key:
447: . -ts_pseudo_increment <increment> - set pseudo increment
449: Level: advanced
451: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
452: @*/
453: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
454: {
455: PetscFunctionBegin;
458: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@
463: TSPseudoSetMaxTimeStep - Sets the maximum time step
464: when using the TSPseudoTimeStepDefault() routine.
466: Logically Collective
468: Input Parameters:
469: + ts - the timestep context
470: - maxdt - the maximum time step, use a non-positive value to deactivate
472: Options Database Key:
473: . -ts_pseudo_max_dt <increment> - set pseudo max dt
475: Level: advanced
477: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
478: @*/
479: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
480: {
481: PetscFunctionBegin;
484: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: /*@
489: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
490: is computed via the formula $ dt = initial\_dt*initial\_fnorm/current\_fnorm $ rather than the default update, $ dt = current\_dt*previous\_fnorm/current\_fnorm.$
492: Logically Collective
494: Input Parameter:
495: . ts - the timestep context
497: Options Database Key:
498: . -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial $dt$ to determine increment
500: Level: advanced
502: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
503: @*/
504: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
505: {
506: PetscFunctionBegin;
508: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*@C
513: TSPseudoSetTimeStep - Sets the user-defined routine to be
514: called at each pseudo-timestep to update the timestep.
516: Logically Collective
518: Input Parameters:
519: + ts - timestep context
520: . dt - function to compute timestep
521: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
523: Calling sequence of `dt`:
524: + ts - the `TS` context
525: . newdt - the newly computed timestep
526: - ctx - [optional] user-defined context
528: Level: intermediate
530: Notes:
531: The routine set here will be called by `TSPseudoComputeTimeStep()`
532: during the timestepping process.
534: If not set then `TSPseudoTimeStepDefault()` is automatically used
536: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
537: @*/
538: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, void *ctx), void *ctx)
539: {
540: PetscFunctionBegin;
542: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode (*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
547: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, void *ctx)
548: {
549: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
551: PetscFunctionBegin;
552: pseudo->verify = dt;
553: pseudo->verifyctx = ctx;
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
558: {
559: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
561: PetscFunctionBegin;
562: pseudo->dt_increment = inc;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
567: {
568: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
570: PetscFunctionBegin;
571: pseudo->dt_max = maxdt;
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
576: {
577: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
579: PetscFunctionBegin;
580: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
585: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, void *ctx)
586: {
587: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
589: PetscFunctionBegin;
590: pseudo->dt = dt;
591: pseudo->dtctx = ctx;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: /*MC
596: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
598: This method solves equations of the form
600: $$
601: F(X,Xdot) = 0
602: $$
604: for steady state using the iteration
606: $$
607: [G'] S = -F(X,0)
608: X += S
609: $$
611: where
613: $$
614: G(Y) = F(Y,(Y-X)/dt)
615: $$
617: This is linearly-implicit Euler with the residual always evaluated "at steady
618: state". See note below.
620: In addition to the modified solve, a dedicated adaptive timestepping scheme is implemented, mimicking the switched evolution relaxation in {cite}`ckk02`.
621: It determines the next timestep via
623: $$
624: dt_{n} = r dt_{n-1} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{n-1},0) \Vert}
625: $$
627: where $r$ is an additional growth factor (set by `-ts_pseudo_increment`).
628: An alternative formulation is also available that uses the initial timestep and function norm.
630: $$
631: dt_{n} = r dt_{0} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{0},0) \Vert}
632: $$
634: This is chosen by setting `-ts_pseudo_increment_dt_from_initial_dt`.
635: For either method, an upper limit on the timestep can be set by `-ts_pseudo_max_dt`.
637: Options Database Keys:
638: + -ts_pseudo_increment <real> - ratio of increase dt
639: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
640: . -ts_pseudo_max_dt - Maximum dt for adaptive timestepping algorithm
641: . -ts_pseudo_monitor - Monitor convergence of the function norm
642: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than `atol`
643: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than `rtol`
645: Level: beginner
647: Notes:
648: The residual computed by this method includes the transient term (Xdot is computed instead of
649: always being zero), but since the prediction from the last step is always the solution from the
650: last step, on the first Newton iteration we have
652: $$
653: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
654: $$
656: The Jacobian $ dF/dX + shift*dF/dXdot $ contains a non-zero $ dF/dXdot$ term. In the $ X' = F(X) $ case it
657: is $ dF/dX + shift*1/dt $. Hence still contains the $ 1/dt $ term so the Jacobian is not simply the
658: Jacobian of $ F $ and thus this pseudo-transient continuation is not just Newton on $F(x)=0$.
660: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
661: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
662: algorithm is no longer the one described in the referenced papers.
663: By default, the `SNESType` is set to `SNESKSPONLY` to match the algorithm from the referenced papers.
664: Setting the `SNESType` via `-snes_type` will override this default setting.
666: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
667: M*/
668: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
669: {
670: TS_Pseudo *pseudo;
671: SNES snes;
672: SNESType stype;
674: PetscFunctionBegin;
675: ts->ops->reset = TSReset_Pseudo;
676: ts->ops->destroy = TSDestroy_Pseudo;
677: ts->ops->view = TSView_Pseudo;
678: ts->ops->setup = TSSetUp_Pseudo;
679: ts->ops->step = TSStep_Pseudo;
680: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
681: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
682: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
683: ts->default_adapt_type = TSADAPTTSPSEUDO;
684: ts->usessnes = PETSC_TRUE;
686: PetscCall(TSAdaptRegister(TSADAPTTSPSEUDO, TSAdaptCreate_TSPseudo));
688: PetscCall(TSGetSNES(ts, &snes));
689: PetscCall(SNESGetType(snes, &stype));
690: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
692: PetscCall(PetscNew(&pseudo));
693: ts->data = (void *)pseudo;
695: pseudo->dt = TSPseudoTimeStepDefault;
696: pseudo->dtctx = NULL;
697: pseudo->dt_increment = 1.1;
698: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
699: pseudo->fnorm = -1;
700: pseudo->fnorm_initial = -1;
701: pseudo->fnorm_previous = -1;
702: #if defined(PETSC_USE_REAL_SINGLE)
703: pseudo->fatol = 1.e-25;
704: pseudo->frtol = 1.e-5;
705: #else
706: pseudo->fatol = 1.e-50;
707: pseudo->frtol = 1.e-12;
708: #endif
709: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
710: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
711: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
712: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
713: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }