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;
27: PetscObjectState Xstate; /* state of input vector to TSComputeIFunction() with 0 Xdot */
28: } TS_Pseudo;
30: /* ------------------------------------------------------------------------------*/
32: /*@
33: TSPseudoComputeTimeStep - Computes the next timestep for a currently running
34: pseudo-timestepping process.
36: Collective
38: Input Parameter:
39: . ts - timestep context
41: Output Parameter:
42: . dt - newly computed timestep
44: Level: developer
46: Note:
47: The routine to be called here to compute the timestep should be
48: set by calling `TSPseudoSetTimeStep()`.
50: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoSetTimeStep()`
51: @*/
52: PetscErrorCode TSPseudoComputeTimeStep(TS ts, PetscReal *dt)
53: {
54: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
56: PetscFunctionBegin;
57: PetscCall(PetscLogEventBegin(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
58: PetscCall((*pseudo->dt)(ts, dt, pseudo->dtctx));
59: PetscCall(PetscLogEventEnd(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: /* ------------------------------------------------------------------------------*/
64: /*@C
65: TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
67: Collective, No Fortran Support
69: Input Parameters:
70: + ts - the timestep context
71: . dtctx - unused timestep context
72: - update - latest solution vector
74: Output Parameters:
75: + newdt - the timestep to use for the next step
76: - flag - flag indicating whether the last time step was acceptable
78: Level: advanced
80: Note:
81: This routine always returns a flag of 1, indicating an acceptable
82: timestep.
84: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetVerifyTimeStep()`, `TSPseudoVerifyTimeStep()`
85: @*/
86: PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts, Vec update, void *dtctx, PetscReal *newdt, PetscBool *flag)
87: {
88: PetscFunctionBegin;
89: *flag = PETSC_TRUE;
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*@
94: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
96: Collective
98: Input Parameters:
99: + ts - timestep context
100: - update - latest solution vector
102: Output Parameters:
103: + dt - newly computed timestep (if it had to shrink)
104: - flag - indicates if current timestep was ok
106: Level: advanced
108: Notes:
109: The routine to be called here to compute the timestep should be
110: set by calling `TSPseudoSetVerifyTimeStep()`.
112: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetVerifyTimeStep()`, `TSPseudoVerifyTimeStepDefault()`
113: @*/
114: PetscErrorCode TSPseudoVerifyTimeStep(TS ts, Vec update, PetscReal *dt, PetscBool *flag)
115: {
116: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
118: PetscFunctionBegin;
119: *flag = PETSC_TRUE;
120: if (pseudo->verify) PetscCall((*pseudo->verify)(ts, update, pseudo->verifyctx, dt, flag));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /* --------------------------------------------------------------------------------*/
126: static PetscErrorCode TSStep_Pseudo(TS ts)
127: {
128: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
129: PetscInt nits, lits, reject;
130: PetscBool stepok;
131: PetscReal next_time_step = ts->time_step;
133: PetscFunctionBegin;
134: if (ts->steps == 0) {
135: pseudo->dt_initial = ts->time_step;
136: PetscCall(VecCopy(ts->vec_sol, pseudo->update)); /* in all future updates pseudo->update already contains the current time solution */
137: }
138: PetscCall(TSPseudoComputeTimeStep(ts, &next_time_step));
139: for (reject = 0; reject < ts->max_reject; reject++, ts->reject++) {
140: ts->time_step = next_time_step;
141: if (reject) PetscCall(VecCopy(ts->vec_sol, pseudo->update)); /* need to copy the solution because we got a rejection and pseudo->update contains intermediate computation */
142: PetscCall(TSPreStage(ts, ts->ptime + ts->time_step));
143: PetscCall(SNESSolve(ts->snes, NULL, pseudo->update));
144: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
145: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
146: ts->snes_its += nits;
147: ts->ksp_its += lits;
148: PetscCall(TSPostStage(ts, ts->ptime + ts->time_step, 0, &pseudo->update));
149: PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, pseudo->update, &stepok));
150: if (!stepok) {
151: next_time_step = ts->time_step;
152: continue;
153: }
154: pseudo->fnorm = -1; /* The current norm is no longer valid */
155: PetscCall(TSPseudoVerifyTimeStep(ts, pseudo->update, &next_time_step, &stepok));
156: if (stepok) break;
157: }
158: if (reject >= ts->max_reject) {
159: ts->reason = TS_DIVERGED_STEP_REJECTED;
160: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, reject));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: PetscCall(VecCopy(pseudo->update, ts->vec_sol));
165: ts->ptime += ts->time_step;
166: ts->time_step = next_time_step;
168: PetscCall(VecZeroEntries(pseudo->xdot));
169: PetscCall(TSComputeIFunction(ts, ts->ptime, pseudo->update, pseudo->xdot, pseudo->func, PETSC_FALSE));
170: PetscCall(PetscObjectStateGet((PetscObject)pseudo->update, &pseudo->Xstate));
171: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
173: if (pseudo->fnorm < pseudo->fatol) {
174: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
175: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->frtol));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
178: if (pseudo->fnorm / pseudo->fnorm_initial < pseudo->frtol) {
179: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
180: 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));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*------------------------------------------------------------*/
188: static PetscErrorCode TSReset_Pseudo(TS ts)
189: {
190: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
192: PetscFunctionBegin;
193: PetscCall(VecDestroy(&pseudo->update));
194: PetscCall(VecDestroy(&pseudo->func));
195: PetscCall(VecDestroy(&pseudo->xdot));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode TSDestroy_Pseudo(TS ts)
200: {
201: PetscFunctionBegin;
202: PetscCall(TSReset_Pseudo(ts));
203: PetscCall(PetscFree(ts->data));
204: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL));
205: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL));
206: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL));
207: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL));
208: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*------------------------------------------------------------*/
214: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
215: {
216: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
218: PetscFunctionBegin;
219: *Xdot = pseudo->xdot;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: The transient residual is
226: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
228: or for ODE,
230: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
232: This is the function that must be evaluated for transient simulation and for
233: finite difference Jacobians. On the first Newton step, this algorithm uses
234: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
235: residual is actually the steady state residual. Pseudotransient
236: continuation as described in the literature is a linearly implicit
237: algorithm, it only takes this one full Newton step with the steady state
238: residual, and then advances to the next time step.
240: This routine avoids a repeated identical call to TSComputeRHSFunction() when that result
241: is already contained in pseudo->func due to the call to TSComputeIFunction() in TSStep_Pseudo()
243: */
244: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
245: {
246: TSIFunctionFn *ifunction;
247: TSRHSFunctionFn *rhsfunction;
248: void *ctx;
249: DM dm;
250: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
251: const PetscScalar mdt = 1.0 / ts->time_step;
252: PetscBool KSPSNES;
253: PetscObjectState Xstate;
255: PetscFunctionBegin;
261: PetscCall(TSGetDM(ts, &dm));
262: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
263: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
264: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
266: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESKSPONLY, &KSPSNES));
267: PetscCall(PetscObjectStateGet((PetscObject)X, &Xstate));
268: if (Xstate != pseudo->Xstate || ifunction || !KSPSNES) {
269: PetscCall(VecAXPBYPCZ(pseudo->xdot, -mdt, mdt, 0, ts->vec_sol, X));
270: PetscCall(TSComputeIFunction(ts, ts->ptime + ts->time_step, X, pseudo->xdot, Y, PETSC_FALSE));
271: } else {
272: /* reuse the TSComputeIFunction() result performed inside TSStep_Pseudo() */
273: /* note that pseudo->func contains the negation of TSComputeRHSFunction() */
274: PetscCall(VecWAXPY(Y, 1, pseudo->func, pseudo->xdot));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*
280: This constructs the Jacobian needed for SNES. For DAE, this is
282: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
284: and for ODE:
286: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
287: */
288: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
289: {
290: Vec Xdot;
292: PetscFunctionBegin;
293: /* Xdot has already been computed in SNESTSFormFunction_Pseudo (SNES guarantees this) */
294: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
295: PetscCall(TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode TSSetUp_Pseudo(TS ts)
300: {
301: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
303: PetscFunctionBegin;
304: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->update));
305: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->func));
306: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->xdot));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*------------------------------------------------------------*/
312: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
313: {
314: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
315: PetscViewer viewer = (PetscViewer)dummy;
317: PetscFunctionBegin;
318: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
319: PetscCall(VecZeroEntries(pseudo->xdot));
320: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
321: PetscCall(PetscObjectStateGet((PetscObject)ts->vec_sol, &pseudo->Xstate));
322: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
323: }
324: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
325: PetscCall(PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)pseudo->fnorm));
326: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems PetscOptionsObject)
331: {
332: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
333: PetscBool flg = PETSC_FALSE;
334: PetscViewer viewer;
336: PetscFunctionBegin;
337: PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
338: PetscCall(PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL));
339: if (flg) {
340: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer));
341: PetscCall(TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscCtxDestroyFn *)PetscViewerDestroy));
342: }
343: flg = pseudo->increment_dt_from_initial_dt;
344: PetscCall(PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL));
345: pseudo->increment_dt_from_initial_dt = flg;
346: PetscCall(PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL));
347: PetscCall(PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL));
348: PetscCall(PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL));
349: PetscCall(PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL));
350: PetscOptionsHeadEnd();
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
355: {
356: PetscBool isascii;
358: PetscFunctionBegin;
359: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
360: if (isascii) {
361: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
362: PetscCall(PetscViewerASCIIPrintf(viewer, " frtol - relative tolerance in function value %g\n", (double)pseudo->frtol));
363: PetscCall(PetscViewerASCIIPrintf(viewer, " fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol));
364: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_initial - initial timestep %g\n", (double)pseudo->dt_initial));
365: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment));
366: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_max - maximum time %g\n", (double)pseudo->dt_max));
367: }
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: /*@C
372: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
373: last timestep.
375: Logically Collective
377: Input Parameters:
378: + ts - timestep context
379: . dt - user-defined function to verify timestep
380: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
382: Calling sequence of `func`:
383: + ts - the time-step context
384: . update - latest solution vector
385: . ctx - [optional] user-defined timestep context
386: . newdt - the timestep to use for the next step
387: - flag - flag indicating whether the last time step was acceptable
389: Level: advanced
391: Note:
392: The routine set here will be called by `TSPseudoVerifyTimeStep()`
393: during the timestepping process.
395: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
396: @*/
397: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, void *ctx, PetscReal *newdt, PetscBool *flag), void *ctx)
398: {
399: PetscFunctionBegin;
401: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode (*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: /*@
406: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
407: dt when using the TSPseudoTimeStepDefault() routine.
409: Logically Collective
411: Input Parameters:
412: + ts - the timestep context
413: - inc - the scaling factor >= 1.0
415: Options Database Key:
416: . -ts_pseudo_increment <increment> - set pseudo increment
418: Level: advanced
420: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
421: @*/
422: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
423: {
424: PetscFunctionBegin;
427: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: /*@
432: TSPseudoSetMaxTimeStep - Sets the maximum time step
433: when using the TSPseudoTimeStepDefault() routine.
435: Logically Collective
437: Input Parameters:
438: + ts - the timestep context
439: - maxdt - the maximum time step, use a non-positive value to deactivate
441: Options Database Key:
442: . -ts_pseudo_max_dt <increment> - set pseudo max dt
444: Level: advanced
446: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
447: @*/
448: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
449: {
450: PetscFunctionBegin;
453: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
459: is computed via the formula $ dt = initial\_dt*initial\_fnorm/current\_fnorm $ rather than the default update, $ dt = current\_dt*previous\_fnorm/current\_fnorm.$
461: Logically Collective
463: Input Parameter:
464: . ts - the timestep context
466: Options Database Key:
467: . -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial $dt$ to determine increment
469: Level: advanced
471: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
472: @*/
473: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
474: {
475: PetscFunctionBegin;
477: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@C
482: TSPseudoSetTimeStep - Sets the user-defined routine to be
483: called at each pseudo-timestep to update the timestep.
485: Logically Collective
487: Input Parameters:
488: + ts - timestep context
489: . dt - function to compute timestep
490: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
492: Calling sequence of `dt`:
493: + ts - the `TS` context
494: . newdt - the newly computed timestep
495: - ctx - [optional] user-defined context
497: Level: intermediate
499: Notes:
500: The routine set here will be called by `TSPseudoComputeTimeStep()`
501: during the timestepping process.
503: If not set then `TSPseudoTimeStepDefault()` is automatically used
505: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
506: @*/
507: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, void *ctx), void *ctx)
508: {
509: PetscFunctionBegin;
511: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode (*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: /* ----------------------------------------------------------------------------- */
517: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
518: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, void *ctx)
519: {
520: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
522: PetscFunctionBegin;
523: pseudo->verify = dt;
524: pseudo->verifyctx = ctx;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
529: {
530: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
532: PetscFunctionBegin;
533: pseudo->dt_increment = inc;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
538: {
539: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
541: PetscFunctionBegin;
542: pseudo->dt_max = maxdt;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
547: {
548: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
550: PetscFunctionBegin;
551: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
556: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, void *ctx)
557: {
558: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
560: PetscFunctionBegin;
561: pseudo->dt = dt;
562: pseudo->dtctx = ctx;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*MC
567: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
569: This method solves equations of the form
571: $$
572: F(X,Xdot) = 0
573: $$
575: for steady state using the iteration
577: $$
578: [G'] S = -F(X,0)
579: X += S
580: $$
582: where
584: $$
585: G(Y) = F(Y,(Y-X)/dt)
586: $$
588: This is linearly-implicit Euler with the residual always evaluated "at steady
589: state". See note below.
591: In addition to the modified solve, a dedicated adaptive timestepping scheme is implemented, mimicking the switched evolution relaxation in {cite}`ckk02`.
592: It determines the next timestep via
594: $$
595: dt_{n} = r dt_{n-1} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{n-1},0) \Vert}
596: $$
598: where $r$ is an additional growth factor (set by `-ts_pseudo_increment`).
599: An alternative formulation is also available that uses the initial timestep and function norm.
601: $$
602: dt_{n} = r dt_{0} \frac{\Vert F(X_{n},0) \Vert}{\Vert F(X_{0},0) \Vert}
603: $$
605: This is chosen by setting `-ts_pseudo_increment_dt_from_initial_dt`.
606: For either method, an upper limit on the timestep can be set by `-ts_pseudo_max_dt`.
608: Options Database Keys:
609: + -ts_pseudo_increment <real> - ratio of increase dt
610: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
611: . -ts_pseudo_max_dt - Maximum dt for adaptive timestepping algorithm
612: . -ts_pseudo_monitor - Monitor convergence of the function norm
613: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than `atol`
614: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than `rtol`
616: Level: beginner
618: Notes:
619: The residual computed by this method includes the transient term (Xdot is computed instead of
620: always being zero), but since the prediction from the last step is always the solution from the
621: last step, on the first Newton iteration we have
623: $$
624: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
625: $$
627: The Jacobian $ dF/dX + shift*dF/dXdot $ contains a non-zero $ dF/dXdot$ term. In the $ X' = F(X) $ case it
628: is $ dF/dX + shift*1/dt $. Hence still contains the $ 1/dt $ term so the Jacobian is not simply the
629: Jacobian of $ F $ and thus this pseudo-transient continuation is not just Newton on $F(x)=0$.
631: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
632: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
633: algorithm is no longer the one described in the referenced papers.
634: By default, the `SNESType` is set to `SNESKSPONLY` to match the algorithm from the referenced papers.
635: Setting the `SNESType` via `-snes_type` will override this default setting.
637: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
638: M*/
639: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
640: {
641: TS_Pseudo *pseudo;
642: SNES snes;
643: SNESType stype;
645: PetscFunctionBegin;
646: ts->ops->reset = TSReset_Pseudo;
647: ts->ops->destroy = TSDestroy_Pseudo;
648: ts->ops->view = TSView_Pseudo;
649: ts->ops->setup = TSSetUp_Pseudo;
650: ts->ops->step = TSStep_Pseudo;
651: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
652: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
653: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
654: ts->default_adapt_type = TSADAPTNONE;
655: ts->usessnes = PETSC_TRUE;
657: PetscCall(TSGetSNES(ts, &snes));
658: PetscCall(SNESGetType(snes, &stype));
659: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
661: PetscCall(PetscNew(&pseudo));
662: ts->data = (void *)pseudo;
664: pseudo->dt = TSPseudoTimeStepDefault;
665: pseudo->dtctx = NULL;
666: pseudo->dt_increment = 1.1;
667: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
668: pseudo->fnorm = -1;
669: pseudo->fnorm_initial = -1;
670: pseudo->fnorm_previous = -1;
671: #if defined(PETSC_USE_REAL_SINGLE)
672: pseudo->fatol = 1.e-25;
673: pseudo->frtol = 1.e-5;
674: #else
675: pseudo->fatol = 1.e-50;
676: pseudo->frtol = 1.e-12;
677: #endif
678: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
679: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
680: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
681: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
682: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /*@C
687: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
689: Collective, No Fortran Support
691: Input Parameters:
692: + ts - the timestep context
693: - dtctx - unused timestep context
695: Output Parameter:
696: . newdt - the timestep to use for the next step
698: Level: advanced
700: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
701: @*/
702: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
703: {
704: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
705: PetscReal inc = pseudo->dt_increment;
707: PetscFunctionBegin;
708: if (pseudo->fnorm < 0.0) {
709: PetscCall(VecZeroEntries(pseudo->xdot));
710: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
711: PetscCall(PetscObjectStateGet((PetscObject)ts->vec_sol, &pseudo->Xstate));
712: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
713: }
714: if (pseudo->fnorm_initial < 0) {
715: /* first time through so compute initial function norm */
716: pseudo->fnorm_initial = pseudo->fnorm;
717: pseudo->fnorm_previous = pseudo->fnorm;
718: }
719: if (pseudo->fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
720: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / pseudo->fnorm;
721: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / pseudo->fnorm;
722: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
723: pseudo->fnorm_previous = pseudo->fnorm;
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }