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 update; /* work vector where new solution is formed */
10: Vec xdot; /* work vector for time derivative of state */
12: /* information used for (adaptive) 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_previous; /* original and previous norm of F(u) */
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;
27: TSStepStatus status;
28: } TS_Pseudo;
30: #define TSPSEUDO_RESIDUAL_KEY "TSPSEUDO Residual Key"
32: /*
33: This struct is attached to any vector that has a solution in it.
34: Use and copying of this struct between Vecs allows residual information to be carried between ts->vec_sol and pseudo->update
35: */
36: typedef struct {
37: Vec func; /* work vector where F(t[i],u[i]) is stored */
38: PetscObjectState Xstate; /* state of vector given to TSComputeIFunction() with 0 Xdot to compute `residual`*/
39: } TS_Pseudo_Residual;
41: static PetscErrorCode TSPseudoResidualDestroy(PetscCtx ctx)
42: {
43: TS_Pseudo_Residual *pseudo_residual = *(TS_Pseudo_Residual **)ctx;
45: PetscFunctionBegin;
46: PetscCall(VecDestroy(&pseudo_residual->func));
47: PetscCall(PetscFree(pseudo_residual));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*@C
52: TSPseudoComputeFunction - Compute nonlinear residual for pseudo-timestepping
54: This computes the residual for $\dot U = 0$, i.e. $F(U, 0)$ for the IFunction.
56: Collective
58: Input Parameters:
59: + ts - the timestep context
60: - solution - the solution vector
62: Output Parameter:
63: + residual - the nonlinear residual
64: - fnorm - the norm of the nonlinear residual
66: Level: advanced
68: Note:
69: `TSPSEUDO` records the nonlinear residual and the `solution` vector used to generate it. If given the same `solution` vector (as determined by the vector's `PetscObjectState`), this function will return those recorded values.
71: This can be used in a custom adaptive timestepping implementation that needs access to the residual, but can reuse the calculation already done by `TSPSEUDO`.
73: To correctly get the residual reuse behavior, `solution` must be the same `Vec` that was returned by `TSGetSolution()` or the `Vec` given by `TSAdaptCheckStage()`.
75: .seealso: [](ch_ts), `TSPSEUDO`
76: @*/
77: PetscErrorCode TSPseudoComputeFunction(TS ts, Vec solution, Vec *residual, PetscReal *fnorm)
78: {
79: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
80: TS_Pseudo_Residual *pseudo_residual;
81: PetscObjectState Xstate;
83: PetscFunctionBegin;
86: if (residual) PetscAssertPointer(residual, 3);
87: if (fnorm) PetscAssertPointer(fnorm, 4);
88: PetscCheckTypeName(ts, TSPSEUDO);
90: PetscCall(PetscObjectStateGet((PetscObject)solution, &Xstate));
91: PetscCall(PetscObjectContainerQuery((PetscObject)solution, TSPSEUDO_RESIDUAL_KEY, &pseudo_residual));
92: if (!pseudo_residual) {
93: PetscCall(PetscNew(&pseudo_residual));
94: PetscCall(VecDuplicate(solution, &pseudo_residual->func));
95: pseudo_residual->Xstate = -1;
96: PetscCall(PetscObjectContainerCompose((PetscObject)solution, TSPSEUDO_RESIDUAL_KEY, pseudo_residual, TSPseudoResidualDestroy));
97: }
98: if (pseudo_residual->Xstate != Xstate) {
99: PetscCall(VecZeroEntries(pseudo->xdot));
100: PetscCall(TSComputeIFunction(ts, ts->ptime, solution, pseudo->xdot, pseudo_residual->func, PETSC_FALSE));
101: pseudo_residual->Xstate = Xstate;
102: }
104: if (residual) *residual = pseudo_residual->func;
105: if (fnorm) PetscCall(VecNorm(pseudo_residual->func, NORM_2, fnorm));
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: // Copy residual information from src to dest (if said information is valid)
110: // This allows residual information to be copied between ts->vec_sol and pseudo->update
111: static PetscErrorCode TSPseudoCopyResidualInfo(Vec src, Vec dest)
112: {
113: PetscObjectState src_state, dest_state;
114: TS_Pseudo_Residual *src_pseudo_residual = NULL, *dest_pseudo_residual = NULL;
116: PetscFunctionBegin;
117: PetscCall(PetscObjectStateGet((PetscObject)src, &src_state));
118: PetscCall(PetscObjectContainerQuery((PetscObject)src, TSPSEUDO_RESIDUAL_KEY, &src_pseudo_residual));
119: PetscCheck(src_pseudo_residual, PetscObjectComm((PetscObject)src), PETSC_ERR_ARG_WRONGSTATE, "Source vector should have TSPSEUDO residual information attached to it before copying");
121: if (src_pseudo_residual->Xstate != src_state) PetscFunctionReturn(PETSC_SUCCESS);
123: PetscCall(PetscObjectContainerQuery((PetscObject)dest, TSPSEUDO_RESIDUAL_KEY, &dest_pseudo_residual));
124: if (!dest_pseudo_residual) {
125: PetscCall(PetscNew(&dest_pseudo_residual));
126: PetscCall(VecDuplicate(dest, &dest_pseudo_residual->func));
127: PetscCall(PetscObjectContainerCompose((PetscObject)dest, TSPSEUDO_RESIDUAL_KEY, dest_pseudo_residual, TSPseudoResidualDestroy));
128: }
129: PetscCall(PetscObjectStateGet((PetscObject)dest, &dest_state));
130: PetscCall(VecCopy(src_pseudo_residual->func, dest_pseudo_residual->func));
131: dest_pseudo_residual->Xstate = dest_state;
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode TSStep_Pseudo(TS ts)
136: {
137: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
138: PetscInt nits, lits, rejections = 0, step_num;
139: PetscBool accept;
140: PetscReal next_time_step = ts->time_step, fnorm;
141: TSAdapt adapt;
143: PetscFunctionBegin;
144: PetscCall(TSGetStepNumber(ts, &step_num));
145: if (ts->start_step == step_num) {
146: pseudo->dt_initial = ts->time_step;
147: PetscCall(VecCopy(ts->vec_sol, pseudo->update));
148: }
150: pseudo->status = TS_STEP_INCOMPLETE;
151: while (!ts->reason && pseudo->status != TS_STEP_COMPLETE) {
152: PetscCall(TSPreStage(ts, ts->ptime + ts->time_step));
153: PetscCall(SNESSolve(ts->snes, NULL, pseudo->update));
154: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
155: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
156: ts->snes_its += nits;
157: ts->ksp_its += lits;
158: PetscCall(TSPostStage(ts, ts->ptime + ts->time_step, 0, &pseudo->update));
159: PetscCall(TSGetAdapt(ts, &adapt));
160: PetscCall(TSAdaptCheckStage(adapt, ts, ts->ptime + ts->time_step, pseudo->update, &accept));
161: if (!accept) goto reject_step;
163: pseudo->status = TS_STEP_PENDING;
164: PetscCall(VecCopy(pseudo->update, ts->vec_sol));
165: // Copy residual info in case TSPseudoComputeFunction() was called in TSAdaptCheckStage()
166: PetscCall(TSPseudoCopyResidualInfo(pseudo->update, ts->vec_sol));
167: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
168: pseudo->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
169: if (!accept) {
170: ts->time_step = next_time_step;
171: goto reject_step;
172: }
173: ts->ptime += ts->time_step;
174: ts->time_step = next_time_step;
175: break;
177: reject_step:
178: ts->reject++;
179: accept = PETSC_FALSE;
180: PetscCall(VecCopy(ts->vec_sol0, pseudo->update));
181: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
182: ts->reason = TS_DIVERGED_STEP_REJECTED;
183: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
186: }
188: // Check solution convergence
189: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
190: // Copy residual info back to pseudo->update. vec_sol and update have the same values, so no need for VecCopy on them
191: PetscCall(TSPseudoCopyResidualInfo(ts->vec_sol, pseudo->update));
192: if (pseudo->fnorm_initial == -1) pseudo->fnorm_initial = fnorm;
194: if (fnorm < pseudo->fatol) {
195: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
196: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)fnorm, (double)pseudo->fatol));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
199: if (fnorm / pseudo->fnorm_initial < pseudo->frtol) {
200: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
201: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g / fnorm_initial %g < frtol %g\n", ts->steps, (double)fnorm, (double)pseudo->fnorm_initial, (double)pseudo->frtol));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode TSReset_Pseudo(TS ts)
208: {
209: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
211: PetscFunctionBegin;
212: PetscCall(VecDestroy(&pseudo->update));
213: PetscCall(VecDestroy(&pseudo->xdot));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode TSDestroy_Pseudo(TS ts)
218: {
219: PetscFunctionBegin;
220: PetscCall(TSReset_Pseudo(ts));
221: PetscCall(PetscFree(ts->data));
222: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL));
223: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL));
224: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL));
225: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL));
226: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
231: {
232: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
234: PetscFunctionBegin;
235: *Xdot = pseudo->xdot;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*
240: The transient residual is
242: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
244: or for ODE,
246: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
248: This is the function that must be evaluated for transient simulation and for
249: finite difference Jacobians. On the first Newton step, this algorithm uses
250: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
251: residual is actually the steady state residual. Pseudotransient
252: continuation as described in the literature is a linearly implicit
253: algorithm, it only takes this one full Newton step with the steady state
254: residual, and then advances to the next time step.
256: This routine avoids a repeated identical call to TSComputeRHSFunction() when that result
257: is already contained in an attached TS_Pseudo_Residual struct
259: */
260: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
261: {
262: TSIFunctionFn *ifunction;
263: TSRHSFunctionFn *rhsfunction;
264: DM dm;
265: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
266: const PetscScalar mdt = 1.0 / ts->time_step;
267: PetscInt snes_it, snes_func_evals;
269: PetscFunctionBegin;
275: PetscCall(TSGetDM(ts, &dm));
276: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
277: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
278: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
280: PetscCall(SNESGetIterationNumber(snes, &snes_it));
281: PetscCall(SNESGetNumberFunctionEvals(snes, &snes_func_evals));
282: if (snes_it == 0 && snes_func_evals == 0) {
283: Vec func;
284: PetscCall(TSPseudoComputeFunction(ts, X, &func, NULL));
285: PetscCall(VecCopy(func, Y));
286: } else {
287: PetscCall(VecAXPBYPCZ(pseudo->xdot, -mdt, mdt, 0, ts->vec_sol0, X));
288: PetscCall(TSComputeIFunction(ts, ts->ptime + ts->time_step, X, pseudo->xdot, Y, PETSC_FALSE));
289: }
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*
294: This constructs the Jacobian needed for SNES. For DAE, this is
296: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
298: and for ODE:
300: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
301: */
302: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
303: {
304: Vec Xdot;
306: PetscFunctionBegin;
307: /* Xdot has already been computed in SNESTSFormFunction_Pseudo (SNES guarantees this) */
308: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
309: PetscCall(TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode TSSetUp_Pseudo(TS ts)
314: {
315: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
317: PetscFunctionBegin;
318: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->update));
319: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->xdot));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
324: {
325: PetscViewer viewer = (PetscViewer)dummy;
326: PetscReal fnorm;
328: PetscFunctionBegin;
329: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
330: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
331: PetscCall(PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)fnorm));
332: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems PetscOptionsObject)
337: {
338: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
339: PetscBool flg = PETSC_FALSE;
340: PetscViewer viewer;
342: PetscFunctionBegin;
343: PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
344: PetscCall(PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL));
345: if (flg) {
346: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer));
347: PetscCall(TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscCtxDestroyFn *)PetscViewerDestroy));
348: }
349: flg = pseudo->increment_dt_from_initial_dt;
350: PetscCall(PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL));
351: pseudo->increment_dt_from_initial_dt = flg;
352: PetscCall(PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL));
353: PetscCall(PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL));
354: PetscCall(PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL));
355: PetscCall(PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL));
356: PetscOptionsHeadEnd();
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
361: {
362: PetscBool isascii;
364: PetscFunctionBegin;
365: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
366: if (isascii) {
367: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
368: PetscCall(PetscViewerASCIIPrintf(viewer, " frtol - relative tolerance in function value %g\n", (double)pseudo->frtol));
369: PetscCall(PetscViewerASCIIPrintf(viewer, " fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol));
370: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_initial - initial timestep %g\n", (double)pseudo->dt_initial));
371: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment));
372: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_max - maximum time step %g\n", (double)pseudo->dt_max));
373: }
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /*@C
378: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
380: Collective, No Fortran Support
382: Input Parameters:
383: + ts - the timestep context
384: - dtctx - unused timestep context
386: Output Parameter:
387: . newdt - the timestep to use for the next step
389: Level: advanced
391: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
392: @*/
393: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
394: {
395: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
396: PetscReal inc = pseudo->dt_increment, fnorm;
398: PetscFunctionBegin;
399: PetscCall(TSPseudoComputeFunction(ts, ts->vec_sol, NULL, &fnorm));
400: if (pseudo->fnorm_initial < 0) {
401: /* first time through so compute initial function norm */
402: pseudo->fnorm_initial = fnorm;
403: pseudo->fnorm_previous = fnorm;
404: }
405: if (fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
406: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / fnorm;
407: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / fnorm;
408: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
409: pseudo->fnorm_previous = fnorm;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode TSAdaptChoose_TSPseudo(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter)
414: {
415: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
417: PetscFunctionBegin;
418: PetscCall(PetscLogEventBegin(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
419: PetscCallBack("TSPSEUDO callback time step", (*pseudo->dt)(ts, next_h, pseudo->dtctx));
420: PetscCall(PetscLogEventEnd(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
422: *next_sc = 0;
423: *wlte = -1; /* Weighted local truncation error was not evaluated */
424: *wltea = -1; /* Weighted absolute local truncation error was not evaluated */
425: *wlter = -1; /* Weighted relative local truncation error was not evaluated */
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: static PetscErrorCode TSAdaptCheckStage_TSPseudo(TSAdapt adapt, TS ts, PetscReal t, Vec Y, PetscBool *accept)
430: {
431: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
433: PetscFunctionBegin;
434: if (pseudo->verify) {
435: PetscReal dt;
436: PetscCall(TSGetTimeStep(ts, &dt));
437: PetscCallBack("TSPSEUDO callback verify time step", (*pseudo->verify)(ts, Y, pseudo->verifyctx, &dt, accept));
438: PetscCall(TSSetTimeStep(ts, dt));
439: }
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*MC
444: TSADAPTTSPSEUDO - TSPseudo adaptive controller for time stepping
446: Level: developer
448: Note:
449: This is only meant to be used with `TSPSEUDO` time integrator.
451: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptType`, `TSPSEUDO`
452: M*/
453: static PetscErrorCode TSAdaptCreate_TSPseudo(TSAdapt adapt)
454: {
455: PetscFunctionBegin;
456: adapt->ops->choose = TSAdaptChoose_TSPseudo;
457: adapt->checkstage = TSAdaptCheckStage_TSPseudo;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /*@C
462: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
463: last timestep.
465: Logically Collective
467: Input Parameters:
468: + ts - timestep context
469: . dt - user-defined function to verify timestep
470: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
472: Calling sequence of `func`:
473: + ts - the time-step context
474: . update - latest solution vector
475: . ctx - [optional] user-defined timestep context
476: . newdt - the timestep to use for the next step
477: - flag - flag indicating whether the last time step was acceptable
479: Level: advanced
481: Note:
482: The routine set here will be called by `TSPseudoVerifyTimeStep()`
483: during the timestepping process.
485: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
486: @*/
487: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, PetscCtx ctx, PetscReal *newdt, PetscBool *flag), PetscCtx ctx)
488: {
489: PetscFunctionBegin;
491: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode (*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
497: dt when using the TSPseudoTimeStepDefault() routine.
499: Logically Collective
501: Input Parameters:
502: + ts - the timestep context
503: - inc - the scaling factor >= 1.0
505: Options Database Key:
506: . -ts_pseudo_increment increment - set pseudo increment
508: Level: advanced
510: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
511: @*/
512: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
513: {
514: PetscFunctionBegin;
517: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@
522: TSPseudoSetMaxTimeStep - Sets the maximum time step
523: when using the TSPseudoTimeStepDefault() routine.
525: Logically Collective
527: Input Parameters:
528: + ts - the timestep context
529: - maxdt - the maximum time step, use a non-positive value to deactivate
531: Options Database Key:
532: . -ts_pseudo_max_dt increment - set pseudo max dt
534: Level: advanced
536: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
537: @*/
538: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
539: {
540: PetscFunctionBegin;
543: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
549: is computed via the formula $ dt = initial\_dt*initial\_fnorm/current\_fnorm $ rather than the default update, $ dt = current\_dt*previous\_fnorm/current\_fnorm.$
551: Logically Collective
553: Input Parameter:
554: . ts - the timestep context
556: Options Database Key:
557: . -ts_pseudo_increment_dt_from_initial_dt (true|false) - use the initial $dt$ to determine increment
559: Level: advanced
561: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
562: @*/
563: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
564: {
565: PetscFunctionBegin;
567: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*@C
572: TSPseudoSetTimeStep - Sets the user-defined routine to be
573: called at each pseudo-timestep to update the timestep.
575: Logically Collective
577: Input Parameters:
578: + ts - timestep context
579: . dt - function to compute timestep
580: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
582: Calling sequence of `dt`:
583: + ts - the `TS` context
584: . newdt - the newly computed timestep
585: - ctx - [optional] user-defined context
587: Level: intermediate
589: Notes:
590: The routine set here will be called by `TSPseudoComputeTimeStep()`
591: during the timestepping process.
593: If not set then `TSPseudoTimeStepDefault()` is automatically used
595: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
596: @*/
597: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, PetscCtx ctx), PetscCtx ctx)
598: {
599: PetscFunctionBegin;
601: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode (*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
606: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, PetscCtx ctx)
607: {
608: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
610: PetscFunctionBegin;
611: pseudo->verify = dt;
612: pseudo->verifyctx = ctx;
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
617: {
618: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
620: PetscFunctionBegin;
621: pseudo->dt_increment = inc;
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
626: {
627: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
629: PetscFunctionBegin;
630: pseudo->dt_max = maxdt;
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
635: {
636: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
638: PetscFunctionBegin;
639: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
644: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, PetscCtx ctx)
645: {
646: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
648: PetscFunctionBegin;
649: pseudo->dt = dt;
650: pseudo->dtctx = ctx;
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: /*MC
655: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
657: This method solves equations of the form
659: $$
660: F(X,Xdot) = 0
661: $$
663: for steady state using the iteration
665: $$
666: [G'] S = -F(X,0)
667: X += S
668: $$
670: where
672: $$
673: G(Y) = F(Y,(Y-X)/dt)
674: $$
676: This is linearly-implicit Euler with the residual always evaluated "at steady
677: state". See note below.
679: In addition to the modified solve, a dedicated adaptive timestepping scheme is implemented, mimicking the switched evolution relaxation in {cite}`ckk02`.
680: It determines the next timestep via
682: $$
683: dt_{n} = r dt_{n-1} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{n-1},0) \Vert}
684: $$
686: where $r$ is an additional growth factor (set by `-ts_pseudo_increment`).
687: An alternative formulation is also available that uses the initial timestep and function norm.
689: $$
690: dt_{n} = r dt_{0} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{0},0) \Vert}
691: $$
693: This is chosen by setting `-ts_pseudo_increment_dt_from_initial_dt`.
694: For either method, an upper limit on the timestep can be set by `-ts_pseudo_max_dt`.
696: Options Database Keys:
697: + -ts_pseudo_increment real - ratio of increase dt
698: . -ts_pseudo_increment_dt_from_initial_dt (true|false) - Increase dt as a ratio from original dt
699: . -ts_pseudo_max_dt - Maximum `dt` for adaptive timestepping algorithm
700: . -ts_pseudo_monitor - Monitor convergence of the function norm
701: . -ts_pseudo_fatol atol - stop iterating when the function norm is less than `atol`
702: - -ts_pseudo_frtol rtol - stop iterating when the function norm divided by the initial function norm is less than `rtol`
704: Level: beginner
706: Notes:
707: The residual computed by this method includes the transient term (Xdot is computed instead of
708: always being zero), but since the prediction from the last step is always the solution from the
709: last step, on the first Newton iteration we have
711: $$
712: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
713: $$
715: The Jacobian $ dF/dX + shift*dF/dXdot $ contains a non-zero $ dF/dXdot$ term. In the $ X' = F(X) $ case it
716: is $ dF/dX + shift*1/dt $. Hence still contains the $ 1/dt $ term so the Jacobian is not simply the
717: Jacobian of $ F $ and thus this pseudo-transient continuation is not just Newton on $F(x)=0$.
719: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
720: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
721: algorithm is no longer the one described in the referenced papers.
722: By default, the `SNESType` is set to `SNESKSPONLY` to match the algorithm from the referenced papers.
723: Setting the `SNESType` via `-snes_type` will override this default setting.
725: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
726: M*/
727: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
728: {
729: TS_Pseudo *pseudo;
730: SNES snes;
731: SNESType stype;
733: PetscFunctionBegin;
734: ts->ops->reset = TSReset_Pseudo;
735: ts->ops->destroy = TSDestroy_Pseudo;
736: ts->ops->view = TSView_Pseudo;
737: ts->ops->setup = TSSetUp_Pseudo;
738: ts->ops->step = TSStep_Pseudo;
739: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
740: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
741: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
742: ts->default_adapt_type = TSADAPTTSPSEUDO;
743: ts->usessnes = PETSC_TRUE;
745: PetscCall(TSAdaptRegister(TSADAPTTSPSEUDO, TSAdaptCreate_TSPseudo));
747: PetscCall(TSGetSNES(ts, &snes));
748: PetscCall(SNESGetType(snes, &stype));
749: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
751: PetscCall(PetscNew(&pseudo));
752: ts->data = (void *)pseudo;
754: pseudo->dt = TSPseudoTimeStepDefault;
755: pseudo->dtctx = NULL;
756: pseudo->dt_increment = 1.1;
757: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
758: pseudo->fnorm_initial = -1;
759: pseudo->fnorm_previous = -1;
760: #if defined(PETSC_USE_REAL_SINGLE)
761: pseudo->fatol = 1.e-25;
762: pseudo->frtol = 1.e-5;
763: #else
764: pseudo->fatol = 1.e-50;
765: pseudo->frtol = 1.e-12;
766: #endif
767: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
768: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
769: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
770: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
771: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }