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;
67: if (residual) PetscAssertPointer(residual, 3);
68: if (fnorm) PetscAssertPointer(fnorm, 4);
69: PetscCheckTypeName(ts, TSPSEUDO);
71: PetscCall(PetscObjectStateGet((PetscObject)solution, &Xstate));
72: if (Xstate != pseudo->Xstate || pseudo->fnorm < 0) {
73: PetscCall(VecZeroEntries(pseudo->xdot));
74: PetscCall(TSComputeIFunction(ts, ts->ptime, solution, pseudo->xdot, pseudo->func, PETSC_FALSE));
75: pseudo->Xstate = Xstate;
76: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
77: }
78: if (residual) *residual = pseudo->func;
79: if (fnorm) *fnorm = pseudo->fnorm;
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode TSStep_Pseudo(TS ts)
84: {
85: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
86: PetscInt nits, lits, rejections = 0;
87: PetscBool accept;
88: PetscReal next_time_step = ts->time_step, fnorm;
89: TSAdapt adapt;
91: PetscFunctionBegin;
92: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
94: pseudo->status = TS_STEP_INCOMPLETE;
95: while (!ts->reason && pseudo->status != TS_STEP_COMPLETE) {
96: PetscCall(TSPreStage(ts, ts->ptime + ts->time_step));
97: PetscCall(SNESSolve(ts->snes, NULL, ts->vec_sol));
98: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
99: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
100: ts->snes_its += nits;
101: ts->ksp_its += lits;
102: pseudo->fnorm = -1; /* The current norm is no longer valid */
103: PetscCall(TSPostStage(ts, ts->ptime + ts->time_step, 0, &ts->vec_sol));
104: PetscCall(TSGetAdapt(ts, &adapt));
105: PetscCall(TSAdaptCheckStage(adapt, ts, ts->ptime + ts->time_step, ts->vec_sol, &accept));
106: if (!accept) goto reject_step;
108: pseudo->status = TS_STEP_PENDING;
109: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
110: pseudo->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
111: if (!accept) {
112: ts->time_step = next_time_step;
113: goto reject_step;
114: }
115: ts->ptime += ts->time_step;
116: ts->time_step = next_time_step;
117: break;
119: reject_step:
120: ts->reject++;
121: accept = PETSC_FALSE;
122: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
123: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
124: ts->reason = TS_DIVERGED_STEP_REJECTED;
125: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
128: }
130: // Check solution convergence
131: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
132: if (pseudo->fnorm_initial == -1) pseudo->fnorm_initial = fnorm;
134: if (fnorm < pseudo->fatol) {
135: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
136: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->frtol));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
139: if (fnorm / pseudo->fnorm_initial < pseudo->frtol) {
140: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
141: 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));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: static PetscErrorCode TSReset_Pseudo(TS ts)
148: {
149: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
151: PetscFunctionBegin;
152: PetscCall(VecDestroy(&pseudo->func));
153: PetscCall(VecDestroy(&pseudo->xdot));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode TSDestroy_Pseudo(TS ts)
158: {
159: PetscFunctionBegin;
160: PetscCall(TSReset_Pseudo(ts));
161: PetscCall(PetscFree(ts->data));
162: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL));
163: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL));
164: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL));
165: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL));
166: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
171: {
172: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
174: PetscFunctionBegin;
175: *Xdot = pseudo->xdot;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*
180: The transient residual is
182: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
184: or for ODE,
186: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
188: This is the function that must be evaluated for transient simulation and for
189: finite difference Jacobians. On the first Newton step, this algorithm uses
190: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
191: residual is actually the steady state residual. Pseudotransient
192: continuation as described in the literature is a linearly implicit
193: algorithm, it only takes this one full Newton step with the steady state
194: residual, and then advances to the next time step.
196: This routine avoids a repeated identical call to TSComputeRHSFunction() when that result
197: is already contained in pseudo->func due to the call to TSComputeIFunction() in TSStep_Pseudo()
199: */
200: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
201: {
202: TSIFunctionFn *ifunction;
203: TSRHSFunctionFn *rhsfunction;
204: void *ctx;
205: DM dm;
206: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
207: const PetscScalar mdt = 1.0 / ts->time_step;
208: PetscObjectState Xstate;
209: PetscInt snes_it;
211: PetscFunctionBegin;
217: PetscCall(TSGetDM(ts, &dm));
218: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
219: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
220: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
222: PetscCall(PetscObjectStateGet((PetscObject)X, &Xstate));
223: PetscCall(SNESGetIterationNumber(snes, &snes_it));
224: if (Xstate == pseudo->Xstate && snes_it == 1) {
225: /* reuse the TSComputeIFunction() result performed inside TSStep_Pseudo() */
226: if (ifunction) PetscCall(VecCopy(pseudo->func, Y));
227: /* note that pseudo->func contains the negation of TSComputeRHSFunction() */
228: else PetscCall(VecWAXPY(Y, 1, pseudo->func, pseudo->xdot));
229: } else {
230: PetscCall(VecAXPBYPCZ(pseudo->xdot, -mdt, mdt, 0, ts->vec_sol0, X));
231: PetscCall(TSComputeIFunction(ts, ts->ptime + ts->time_step, X, pseudo->xdot, Y, PETSC_FALSE));
232: }
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: /*
237: This constructs the Jacobian needed for SNES. For DAE, this is
239: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
241: and for ODE:
243: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
244: */
245: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
246: {
247: Vec Xdot;
249: PetscFunctionBegin;
250: /* Xdot has already been computed in SNESTSFormFunction_Pseudo (SNES guarantees this) */
251: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
252: PetscCall(TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode TSSetUp_Pseudo(TS ts)
257: {
258: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
260: PetscFunctionBegin;
261: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->func));
262: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->xdot));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
267: {
268: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
269: PetscViewer viewer = (PetscViewer)dummy;
271: PetscFunctionBegin;
272: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, NULL));
273: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
274: PetscCall(PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)pseudo->fnorm));
275: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems PetscOptionsObject)
280: {
281: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
282: PetscBool flg = PETSC_FALSE;
283: PetscViewer viewer;
285: PetscFunctionBegin;
286: PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
287: PetscCall(PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL));
288: if (flg) {
289: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer));
290: PetscCall(TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscCtxDestroyFn *)PetscViewerDestroy));
291: }
292: flg = pseudo->increment_dt_from_initial_dt;
293: PetscCall(PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL));
294: pseudo->increment_dt_from_initial_dt = flg;
295: PetscCall(PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL));
296: PetscCall(PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL));
297: PetscCall(PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL));
298: PetscCall(PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL));
299: PetscOptionsHeadEnd();
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
304: {
305: PetscBool isascii;
307: PetscFunctionBegin;
308: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
309: if (isascii) {
310: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
311: PetscCall(PetscViewerASCIIPrintf(viewer, " frtol - relative tolerance in function value %g\n", (double)pseudo->frtol));
312: PetscCall(PetscViewerASCIIPrintf(viewer, " fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol));
313: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_initial - initial timestep %g\n", (double)pseudo->dt_initial));
314: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment));
315: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_max - maximum time %g\n", (double)pseudo->dt_max));
316: }
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@C
321: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
323: Collective, No Fortran Support
325: Input Parameters:
326: + ts - the timestep context
327: - dtctx - unused timestep context
329: Output Parameter:
330: . newdt - the timestep to use for the next step
332: Level: advanced
334: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
335: @*/
336: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
337: {
338: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
339: PetscReal inc = pseudo->dt_increment, fnorm;
341: PetscFunctionBegin;
342: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
343: if (pseudo->fnorm_initial < 0) {
344: /* first time through so compute initial function norm */
345: pseudo->fnorm_initial = fnorm;
346: pseudo->fnorm_previous = fnorm;
347: }
348: if (fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
349: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / fnorm;
350: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / fnorm;
351: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
352: pseudo->fnorm_previous = fnorm;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode TSAdaptChoose_TSPseudo(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
357: {
358: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
360: PetscFunctionBegin;
361: PetscCall(PetscLogEventBegin(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
362: PetscCallBack("TSPSEUDO callback time step", (*pseudo->dt)(ts, next_h, pseudo->dtctx));
363: PetscCall(PetscLogEventEnd(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
365: *next_sc = 0;
366: *wlte = -1; /* Weighted local truncation error was not evaluated */
367: *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
368: *wlter = -1; /* Weighted relative local truncation error was not evaluated */
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode TSAdaptCheckStage_TSPseudo(TSAdapt adapt, TS ts, PetscReal t, Vec Y, PetscBool *accept)
373: {
374: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
376: PetscFunctionBegin;
377: if (pseudo->verify) {
378: PetscReal dt;
379: PetscCall(TSGetTimeStep(ts, &dt));
380: PetscCallBack("TSPSEUDO callback verify time step", (*pseudo->verify)(ts, Y, pseudo->verifyctx, &dt, accept));
381: PetscCall(TSSetTimeStep(ts, dt));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*MC
387: TSADAPTTSPSEUDO - TSPseudo adaptive controller for time stepping
389: Level: developer
391: Note:
392: This is only meant to be used with `TSPSEUDO` time integrator.
394: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`, `TSPSEUDO`
395: M*/
396: static PetscErrorCode TSAdaptCreate_TSPseudo(TSAdapt adapt)
397: {
398: PetscFunctionBegin;
399: adapt->ops->choose = TSAdaptChoose_TSPseudo;
400: adapt->checkstage = TSAdaptCheckStage_TSPseudo;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@C
405: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
406: last timestep.
408: Logically Collective
410: Input Parameters:
411: + ts - timestep context
412: . dt - user-defined function to verify timestep
413: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
415: Calling sequence of `func`:
416: + ts - the time-step context
417: . update - latest solution vector
418: . ctx - [optional] user-defined timestep context
419: . newdt - the timestep to use for the next step
420: - flag - flag indicating whether the last time step was acceptable
422: Level: advanced
424: Note:
425: The routine set here will be called by `TSPseudoVerifyTimeStep()`
426: during the timestepping process.
428: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
429: @*/
430: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, PetscCtx ctx, PetscReal *newdt, PetscBool *flag), PetscCtx ctx)
431: {
432: PetscFunctionBegin;
434: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode (*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
440: dt when using the TSPseudoTimeStepDefault() routine.
442: Logically Collective
444: Input Parameters:
445: + ts - the timestep context
446: - inc - the scaling factor >= 1.0
448: Options Database Key:
449: . -ts_pseudo_increment increment - set pseudo increment
451: Level: advanced
453: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
454: @*/
455: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
456: {
457: PetscFunctionBegin;
460: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: TSPseudoSetMaxTimeStep - Sets the maximum time step
466: when using the TSPseudoTimeStepDefault() routine.
468: Logically Collective
470: Input Parameters:
471: + ts - the timestep context
472: - maxdt - the maximum time step, use a non-positive value to deactivate
474: Options Database Key:
475: . -ts_pseudo_max_dt increment - set pseudo max dt
477: Level: advanced
479: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
480: @*/
481: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
482: {
483: PetscFunctionBegin;
486: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
492: is computed via the formula $ dt = initial\_dt*initial\_fnorm/current\_fnorm $ rather than the default update, $ dt = current\_dt*previous\_fnorm/current\_fnorm.$
494: Logically Collective
496: Input Parameter:
497: . ts - the timestep context
499: Options Database Key:
500: . -ts_pseudo_increment_dt_from_initial_dt (true|false) - use the initial $dt$ to determine increment
502: Level: advanced
504: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
505: @*/
506: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
507: {
508: PetscFunctionBegin;
510: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*@C
515: TSPseudoSetTimeStep - Sets the user-defined routine to be
516: called at each pseudo-timestep to update the timestep.
518: Logically Collective
520: Input Parameters:
521: + ts - timestep context
522: . dt - function to compute timestep
523: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
525: Calling sequence of `dt`:
526: + ts - the `TS` context
527: . newdt - the newly computed timestep
528: - ctx - [optional] user-defined context
530: Level: intermediate
532: Notes:
533: The routine set here will be called by `TSPseudoComputeTimeStep()`
534: during the timestepping process.
536: If not set then `TSPseudoTimeStepDefault()` is automatically used
538: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
539: @*/
540: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, PetscCtx ctx), PetscCtx ctx)
541: {
542: PetscFunctionBegin;
544: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode (*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
549: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, PetscCtx ctx)
550: {
551: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
553: PetscFunctionBegin;
554: pseudo->verify = dt;
555: pseudo->verifyctx = ctx;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
560: {
561: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
563: PetscFunctionBegin;
564: pseudo->dt_increment = inc;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
569: {
570: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
572: PetscFunctionBegin;
573: pseudo->dt_max = maxdt;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
578: {
579: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
581: PetscFunctionBegin;
582: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
587: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, PetscCtx ctx)
588: {
589: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
591: PetscFunctionBegin;
592: pseudo->dt = dt;
593: pseudo->dtctx = ctx;
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*MC
598: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
600: This method solves equations of the form
602: $$
603: F(X,Xdot) = 0
604: $$
606: for steady state using the iteration
608: $$
609: [G'] S = -F(X,0)
610: X += S
611: $$
613: where
615: $$
616: G(Y) = F(Y,(Y-X)/dt)
617: $$
619: This is linearly-implicit Euler with the residual always evaluated "at steady
620: state". See note below.
622: In addition to the modified solve, a dedicated adaptive timestepping scheme is implemented, mimicking the switched evolution relaxation in {cite}`ckk02`.
623: It determines the next timestep via
625: $$
626: dt_{n} = r dt_{n-1} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{n-1},0) \Vert}
627: $$
629: where $r$ is an additional growth factor (set by `-ts_pseudo_increment`).
630: An alternative formulation is also available that uses the initial timestep and function norm.
632: $$
633: dt_{n} = r dt_{0} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{0},0) \Vert}
634: $$
636: This is chosen by setting `-ts_pseudo_increment_dt_from_initial_dt`.
637: For either method, an upper limit on the timestep can be set by `-ts_pseudo_max_dt`.
639: Options Database Keys:
640: + -ts_pseudo_increment real - ratio of increase dt
641: . -ts_pseudo_increment_dt_from_initial_dt (true|false) - Increase dt as a ratio from original dt
642: . -ts_pseudo_max_dt - Maximum `dt` for adaptive timestepping algorithm
643: . -ts_pseudo_monitor - Monitor convergence of the function norm
644: . -ts_pseudo_fatol atol - stop iterating when the function norm is less than `atol`
645: - -ts_pseudo_frtol rtol - stop iterating when the function norm divided by the initial function norm is less than `rtol`
647: Level: beginner
649: Notes:
650: The residual computed by this method includes the transient term (Xdot is computed instead of
651: always being zero), but since the prediction from the last step is always the solution from the
652: last step, on the first Newton iteration we have
654: $$
655: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
656: $$
658: The Jacobian $ dF/dX + shift*dF/dXdot $ contains a non-zero $ dF/dXdot$ term. In the $ X' = F(X) $ case it
659: is $ dF/dX + shift*1/dt $. Hence still contains the $ 1/dt $ term so the Jacobian is not simply the
660: Jacobian of $ F $ and thus this pseudo-transient continuation is not just Newton on $F(x)=0$.
662: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
663: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
664: algorithm is no longer the one described in the referenced papers.
665: By default, the `SNESType` is set to `SNESKSPONLY` to match the algorithm from the referenced papers.
666: Setting the `SNESType` via `-snes_type` will override this default setting.
668: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
669: M*/
670: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
671: {
672: TS_Pseudo *pseudo;
673: SNES snes;
674: SNESType stype;
676: PetscFunctionBegin;
677: ts->ops->reset = TSReset_Pseudo;
678: ts->ops->destroy = TSDestroy_Pseudo;
679: ts->ops->view = TSView_Pseudo;
680: ts->ops->setup = TSSetUp_Pseudo;
681: ts->ops->step = TSStep_Pseudo;
682: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
683: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
684: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
685: ts->default_adapt_type = TSADAPTTSPSEUDO;
686: ts->usessnes = PETSC_TRUE;
688: PetscCall(TSAdaptRegister(TSADAPTTSPSEUDO, TSAdaptCreate_TSPseudo));
690: PetscCall(TSGetSNES(ts, &snes));
691: PetscCall(SNESGetType(snes, &stype));
692: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
694: PetscCall(PetscNew(&pseudo));
695: ts->data = (void *)pseudo;
697: pseudo->dt = TSPseudoTimeStepDefault;
698: pseudo->dtctx = NULL;
699: pseudo->dt_increment = 1.1;
700: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
701: pseudo->fnorm = -1;
702: pseudo->fnorm_initial = -1;
703: pseudo->fnorm_previous = -1;
704: #if defined(PETSC_USE_REAL_SINGLE)
705: pseudo->fatol = 1.e-25;
706: pseudo->frtol = 1.e-5;
707: #else
708: pseudo->fatol = 1.e-50;
709: pseudo->frtol = 1.e-12;
710: #endif
711: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
712: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
713: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
714: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
715: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }