Actual source code: posindep.c
1: /*
2: Code for Timestepping with implicit backwards Euler.
3: */
4: #include <petsc/private/tsimpl.h>
6: typedef struct {
7: Vec update; /* work vector where new solution is formed */
8: Vec func; /* work vector where F(t[i],u[i]) is stored */
9: Vec xdot; /* work vector for time derivative of state */
11: /* information used for Pseudo-timestepping */
13: PetscErrorCode (*dt)(TS, PetscReal *, void *); /* compute next timestep, and related context */
14: void *dtctx;
15: PetscErrorCode (*verify)(TS, Vec, void *, PetscReal *, PetscBool *); /* verify previous timestep and related context */
16: void *verifyctx;
18: PetscReal fnorm_initial, fnorm; /* original and current norm of F(u) */
19: PetscReal fnorm_previous;
21: PetscReal dt_initial; /* initial time-step */
22: PetscReal dt_increment; /* scaling that dt is incremented each time-step */
23: PetscReal dt_max; /* maximum time step */
24: PetscBool increment_dt_from_initial_dt;
25: PetscReal fatol, frtol;
26: } TS_Pseudo;
28: /* ------------------------------------------------------------------------------*/
30: /*@
31: TSPseudoComputeTimeStep - Computes the next timestep for a currently running
32: pseudo-timestepping process.
34: Collective
36: Input Parameter:
37: . ts - timestep context
39: Output Parameter:
40: . dt - newly computed timestep
42: Level: developer
44: Note:
45: The routine to be called here to compute the timestep should be
46: set by calling `TSPseudoSetTimeStep()`.
48: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoSetTimeStep()`
49: @*/
50: PetscErrorCode TSPseudoComputeTimeStep(TS ts, PetscReal *dt)
51: {
52: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
54: PetscFunctionBegin;
55: PetscCall(PetscLogEventBegin(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
56: PetscCall((*pseudo->dt)(ts, dt, pseudo->dtctx));
57: PetscCall(PetscLogEventEnd(TS_PseudoComputeTimeStep, ts, 0, 0, 0));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: /* ------------------------------------------------------------------------------*/
62: /*@C
63: TSPseudoVerifyTimeStepDefault - Default code to verify the quality of the last timestep.
65: Collective, No Fortran Support
67: Input Parameters:
68: + ts - the timestep context
69: . dtctx - unused timestep context
70: - update - latest solution vector
72: Output Parameters:
73: + newdt - the timestep to use for the next step
74: - flag - flag indicating whether the last time step was acceptable
76: Level: advanced
78: Note:
79: This routine always returns a flag of 1, indicating an acceptable
80: timestep.
82: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetVerifyTimeStep()`, `TSPseudoVerifyTimeStep()`
83: @*/
84: PetscErrorCode TSPseudoVerifyTimeStepDefault(TS ts, Vec update, void *dtctx, PetscReal *newdt, PetscBool *flag)
85: {
86: PetscFunctionBegin;
87: *flag = PETSC_TRUE;
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: /*@
92: TSPseudoVerifyTimeStep - Verifies whether the last timestep was acceptable.
94: Collective
96: Input Parameters:
97: + ts - timestep context
98: - update - latest solution vector
100: Output Parameters:
101: + dt - newly computed timestep (if it had to shrink)
102: - flag - indicates if current timestep was ok
104: Level: advanced
106: Notes:
107: The routine to be called here to compute the timestep should be
108: set by calling `TSPseudoSetVerifyTimeStep()`.
110: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetVerifyTimeStep()`, `TSPseudoVerifyTimeStepDefault()`
111: @*/
112: PetscErrorCode TSPseudoVerifyTimeStep(TS ts, Vec update, PetscReal *dt, PetscBool *flag)
113: {
114: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
116: PetscFunctionBegin;
117: *flag = PETSC_TRUE;
118: if (pseudo->verify) PetscCall((*pseudo->verify)(ts, update, pseudo->verifyctx, dt, flag));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /* --------------------------------------------------------------------------------*/
124: static PetscErrorCode TSStep_Pseudo(TS ts)
125: {
126: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
127: PetscInt nits, lits, reject;
128: PetscBool stepok;
129: PetscReal next_time_step = ts->time_step;
131: PetscFunctionBegin;
132: if (ts->steps == 0) pseudo->dt_initial = ts->time_step;
133: PetscCall(VecCopy(ts->vec_sol, pseudo->update));
134: PetscCall(TSPseudoComputeTimeStep(ts, &next_time_step));
135: for (reject = 0; reject < ts->max_reject; reject++, ts->reject++) {
136: ts->time_step = next_time_step;
137: PetscCall(TSPreStage(ts, ts->ptime + ts->time_step));
138: PetscCall(SNESSolve(ts->snes, NULL, pseudo->update));
139: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
140: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
141: ts->snes_its += nits;
142: ts->ksp_its += lits;
143: PetscCall(TSPostStage(ts, ts->ptime + ts->time_step, 0, &pseudo->update));
144: PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, pseudo->update, &stepok));
145: if (!stepok) {
146: next_time_step = ts->time_step;
147: continue;
148: }
149: pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
150: PetscCall(TSPseudoVerifyTimeStep(ts, pseudo->update, &next_time_step, &stepok));
151: if (stepok) break;
152: }
153: if (reject >= ts->max_reject) {
154: ts->reason = TS_DIVERGED_STEP_REJECTED;
155: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, reject));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: PetscCall(VecCopy(pseudo->update, ts->vec_sol));
160: ts->ptime += ts->time_step;
161: ts->time_step = next_time_step;
163: if (pseudo->fnorm < 0) {
164: PetscCall(VecZeroEntries(pseudo->xdot));
165: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
166: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
167: }
168: if (pseudo->fnorm < pseudo->fatol) {
169: ts->reason = TS_CONVERGED_PSEUDO_FATOL;
170: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", converged since fnorm %g < fatol %g\n", ts->steps, (double)pseudo->fnorm, (double)pseudo->frtol));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
173: if (pseudo->fnorm / pseudo->fnorm_initial < pseudo->frtol) {
174: ts->reason = TS_CONVERGED_PSEUDO_FRTOL;
175: 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));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*------------------------------------------------------------*/
182: static PetscErrorCode TSReset_Pseudo(TS ts)
183: {
184: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
186: PetscFunctionBegin;
187: PetscCall(VecDestroy(&pseudo->update));
188: PetscCall(VecDestroy(&pseudo->func));
189: PetscCall(VecDestroy(&pseudo->xdot));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode TSDestroy_Pseudo(TS ts)
194: {
195: PetscFunctionBegin;
196: PetscCall(TSReset_Pseudo(ts));
197: PetscCall(PetscFree(ts->data));
198: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", NULL));
199: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", NULL));
200: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", NULL));
201: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", NULL));
202: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", NULL));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*------------------------------------------------------------*/
208: /*
209: Compute Xdot = (X^{n+1}-X^n)/dt) = 0
210: */
211: static PetscErrorCode TSPseudoGetXdot(TS ts, Vec X, Vec *Xdot)
212: {
213: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
214: const PetscScalar mdt = 1.0 / ts->time_step, *xnp1, *xn;
215: PetscScalar *xdot;
216: PetscInt i, n;
218: PetscFunctionBegin;
219: *Xdot = NULL;
220: PetscCall(VecGetArrayRead(ts->vec_sol, &xn));
221: PetscCall(VecGetArrayRead(X, &xnp1));
222: PetscCall(VecGetArray(pseudo->xdot, &xdot));
223: PetscCall(VecGetLocalSize(X, &n));
224: for (i = 0; i < n; i++) xdot[i] = mdt * (xnp1[i] - xn[i]);
225: PetscCall(VecRestoreArrayRead(ts->vec_sol, &xn));
226: PetscCall(VecRestoreArrayRead(X, &xnp1));
227: PetscCall(VecRestoreArray(pseudo->xdot, &xdot));
228: *Xdot = pseudo->xdot;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*
233: The transient residual is
235: F(U^{n+1},(U^{n+1}-U^n)/dt) = 0
237: or for ODE,
239: (U^{n+1} - U^{n})/dt - F(U^{n+1}) = 0
241: This is the function that must be evaluated for transient simulation and for
242: finite difference Jacobians. On the first Newton step, this algorithm uses
243: a guess of U^{n+1} = U^n in which case the transient term vanishes and the
244: residual is actually the steady state residual. Pseudotransient
245: continuation as described in the literature is a linearly implicit
246: algorithm, it only takes this one Newton step with the steady state
247: residual, and then advances to the next time step.
248: */
249: static PetscErrorCode SNESTSFormFunction_Pseudo(SNES snes, Vec X, Vec Y, TS ts)
250: {
251: Vec Xdot;
253: PetscFunctionBegin;
254: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
255: PetscCall(TSComputeIFunction(ts, ts->ptime + ts->time_step, X, Xdot, Y, PETSC_FALSE));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*
260: This constructs the Jacobian needed for SNES. For DAE, this is
262: dF(X,Xdot)/dX + shift*dF(X,Xdot)/dXdot
264: and for ODE:
266: J = I/dt - J_{Frhs} where J_{Frhs} is the given Jacobian of Frhs.
267: */
268: static PetscErrorCode SNESTSFormJacobian_Pseudo(SNES snes, Vec X, Mat AA, Mat BB, TS ts)
269: {
270: Vec Xdot;
272: PetscFunctionBegin;
273: PetscCall(TSPseudoGetXdot(ts, X, &Xdot));
274: PetscCall(TSComputeIJacobian(ts, ts->ptime + ts->time_step, X, Xdot, 1. / ts->time_step, AA, BB, PETSC_FALSE));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode TSSetUp_Pseudo(TS ts)
279: {
280: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
282: PetscFunctionBegin;
283: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->update));
284: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->func));
285: PetscCall(VecDuplicate(ts->vec_sol, &pseudo->xdot));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
288: /*------------------------------------------------------------*/
290: static PetscErrorCode TSPseudoMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, void *dummy)
291: {
292: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
293: PetscViewer viewer = (PetscViewer)dummy;
295: PetscFunctionBegin;
296: if (pseudo->fnorm < 0) { /* The last computed norm is stale, recompute */
297: PetscCall(VecZeroEntries(pseudo->xdot));
298: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
299: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
300: }
301: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
302: PetscCall(PetscViewerASCIIPrintf(viewer, "TS %" PetscInt_FMT " dt %g time %g fnorm %g\n", step, (double)ts->time_step, (double)ptime, (double)pseudo->fnorm));
303: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode TSSetFromOptions_Pseudo(TS ts, PetscOptionItems *PetscOptionsObject)
308: {
309: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
310: PetscBool flg = PETSC_FALSE;
311: PetscViewer viewer;
313: PetscFunctionBegin;
314: PetscOptionsHeadBegin(PetscOptionsObject, "Pseudo-timestepping options");
315: PetscCall(PetscOptionsBool("-ts_monitor_pseudo", "Monitor convergence", "", flg, &flg, NULL));
316: if (flg) {
317: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts), "stdout", &viewer));
318: PetscCall(TSMonitorSet(ts, TSPseudoMonitorDefault, viewer, (PetscCtxDestroyFn *)PetscViewerDestroy));
319: }
320: flg = pseudo->increment_dt_from_initial_dt;
321: PetscCall(PetscOptionsBool("-ts_pseudo_increment_dt_from_initial_dt", "Increase dt as a ratio from original dt", "TSPseudoIncrementDtFromInitialDt", flg, &flg, NULL));
322: pseudo->increment_dt_from_initial_dt = flg;
323: PetscCall(PetscOptionsReal("-ts_pseudo_increment", "Ratio to increase dt", "TSPseudoSetTimeStepIncrement", pseudo->dt_increment, &pseudo->dt_increment, NULL));
324: PetscCall(PetscOptionsReal("-ts_pseudo_max_dt", "Maximum value for dt", "TSPseudoSetMaxTimeStep", pseudo->dt_max, &pseudo->dt_max, NULL));
325: PetscCall(PetscOptionsReal("-ts_pseudo_fatol", "Tolerance for norm of function", "", pseudo->fatol, &pseudo->fatol, NULL));
326: PetscCall(PetscOptionsReal("-ts_pseudo_frtol", "Relative tolerance for norm of function", "", pseudo->frtol, &pseudo->frtol, NULL));
327: PetscOptionsHeadEnd();
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode TSView_Pseudo(TS ts, PetscViewer viewer)
332: {
333: PetscBool isascii;
335: PetscFunctionBegin;
336: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
337: if (isascii) {
338: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
339: PetscCall(PetscViewerASCIIPrintf(viewer, " frtol - relative tolerance in function value %g\n", (double)pseudo->frtol));
340: PetscCall(PetscViewerASCIIPrintf(viewer, " fatol - absolute tolerance in function value %g\n", (double)pseudo->fatol));
341: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_initial - initial timestep %g\n", (double)pseudo->dt_initial));
342: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_increment - increase in timestep on successful step %g\n", (double)pseudo->dt_increment));
343: PetscCall(PetscViewerASCIIPrintf(viewer, " dt_max - maximum time %g\n", (double)pseudo->dt_max));
344: }
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@C
349: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
350: last timestep.
352: Logically Collective
354: Input Parameters:
355: + ts - timestep context
356: . dt - user-defined function to verify timestep
357: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
359: Calling sequence of `func`:
360: + ts - the time-step context
361: . update - latest solution vector
362: . ctx - [optional] user-defined timestep context
363: . newdt - the timestep to use for the next step
364: - flag - flag indicating whether the last time step was acceptable
366: Level: advanced
368: Note:
369: The routine set here will be called by `TSPseudoVerifyTimeStep()`
370: during the timestepping process.
372: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
373: @*/
374: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, void *ctx, PetscReal *newdt, PetscBool *flag), void *ctx)
375: {
376: PetscFunctionBegin;
378: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode (*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*@
383: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
384: dt when using the TSPseudoTimeStepDefault() routine.
386: Logically Collective
388: Input Parameters:
389: + ts - the timestep context
390: - inc - the scaling factor >= 1.0
392: Options Database Key:
393: . -ts_pseudo_increment <increment> - set pseudo increment
395: Level: advanced
397: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
398: @*/
399: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
400: {
401: PetscFunctionBegin;
404: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: TSPseudoSetMaxTimeStep - Sets the maximum time step
410: when using the TSPseudoTimeStepDefault() routine.
412: Logically Collective
414: Input Parameters:
415: + ts - the timestep context
416: - maxdt - the maximum time step, use a non-positive value to deactivate
418: Options Database Key:
419: . -ts_pseudo_max_dt <increment> - set pseudo max dt
421: Level: advanced
423: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
424: @*/
425: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
426: {
427: PetscFunctionBegin;
430: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
436: is computed via the formula
437: $ dt = initial_dt*initial_fnorm/current_fnorm
438: rather than the default update,
439: $ dt = current_dt*previous_fnorm/current_fnorm.
441: Logically Collective
443: Input Parameter:
444: . ts - the timestep context
446: Options Database Key:
447: . -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial dt to determine increment
449: Level: advanced
451: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
452: @*/
453: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
454: {
455: PetscFunctionBegin;
457: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /*@C
462: TSPseudoSetTimeStep - Sets the user-defined routine to be
463: called at each pseudo-timestep to update the timestep.
465: Logically Collective
467: Input Parameters:
468: + ts - timestep context
469: . dt - function to compute timestep
470: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
472: Calling sequence of `dt`:
473: + ts - the `TS` context
474: . newdt - the newly computed timestep
475: - ctx - [optional] user-defined context
477: Level: intermediate
479: Notes:
480: The routine set here will be called by `TSPseudoComputeTimeStep()`
481: during the timestepping process.
483: If not set then `TSPseudoTimeStepDefault()` is automatically used
485: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
486: @*/
487: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, void *ctx), void *ctx)
488: {
489: PetscFunctionBegin;
491: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode (*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /* ----------------------------------------------------------------------------- */
497: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
498: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, void *ctx)
499: {
500: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
502: PetscFunctionBegin;
503: pseudo->verify = dt;
504: pseudo->verifyctx = ctx;
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
509: {
510: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
512: PetscFunctionBegin;
513: pseudo->dt_increment = inc;
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
518: {
519: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
521: PetscFunctionBegin;
522: pseudo->dt_max = maxdt;
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
527: {
528: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
530: PetscFunctionBegin;
531: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
536: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, void *ctx)
537: {
538: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
540: PetscFunctionBegin;
541: pseudo->dt = dt;
542: pseudo->dtctx = ctx;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*MC
547: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
549: This method solves equations of the form
551: $$
552: F(X,Xdot) = 0
553: $$
555: for steady state using the iteration
557: $$
558: [G'] S = -F(X,0)
559: X += S
560: $$
562: where
564: $$
565: G(Y) = F(Y,(Y-X)/dt)
566: $$
568: This is linearly-implicit Euler with the residual always evaluated "at steady
569: state". See note below.
571: Options Database Keys:
572: + -ts_pseudo_increment <real> - ratio of increase dt
573: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
574: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
575: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
577: Level: beginner
579: Notes:
580: The residual computed by this method includes the transient term (Xdot is computed instead of
581: always being zero), but since the prediction from the last step is always the solution from the
582: last step, on the first Newton iteration we have
584: $$
585: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
586: $$
588: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
589: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
590: algorithm is no longer the one described in the referenced papers.
592: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
593: M*/
594: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
595: {
596: TS_Pseudo *pseudo;
597: SNES snes;
598: SNESType stype;
600: PetscFunctionBegin;
601: ts->ops->reset = TSReset_Pseudo;
602: ts->ops->destroy = TSDestroy_Pseudo;
603: ts->ops->view = TSView_Pseudo;
604: ts->ops->setup = TSSetUp_Pseudo;
605: ts->ops->step = TSStep_Pseudo;
606: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
607: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
608: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
609: ts->default_adapt_type = TSADAPTNONE;
610: ts->usessnes = PETSC_TRUE;
612: PetscCall(TSGetSNES(ts, &snes));
613: PetscCall(SNESGetType(snes, &stype));
614: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
616: PetscCall(PetscNew(&pseudo));
617: ts->data = (void *)pseudo;
619: pseudo->dt = TSPseudoTimeStepDefault;
620: pseudo->dtctx = NULL;
621: pseudo->dt_increment = 1.1;
622: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
623: pseudo->fnorm = -1;
624: pseudo->fnorm_initial = -1;
625: pseudo->fnorm_previous = -1;
626: #if defined(PETSC_USE_REAL_SINGLE)
627: pseudo->fatol = 1.e-25;
628: pseudo->frtol = 1.e-5;
629: #else
630: pseudo->fatol = 1.e-50;
631: pseudo->frtol = 1.e-12;
632: #endif
633: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
634: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
635: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
636: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
637: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@C
642: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
644: Collective, No Fortran Support
646: Input Parameters:
647: + ts - the timestep context
648: - dtctx - unused timestep context
650: Output Parameter:
651: . newdt - the timestep to use for the next step
653: Level: advanced
655: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
656: @*/
657: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
658: {
659: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
660: PetscReal inc = pseudo->dt_increment;
662: PetscFunctionBegin;
663: PetscCall(VecZeroEntries(pseudo->xdot));
664: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
665: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
666: if (pseudo->fnorm_initial < 0) {
667: /* first time through so compute initial function norm */
668: pseudo->fnorm_initial = pseudo->fnorm;
669: pseudo->fnorm_previous = pseudo->fnorm;
670: }
671: if (pseudo->fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
672: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / pseudo->fnorm;
673: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / pseudo->fnorm;
674: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
675: pseudo->fnorm_previous = pseudo->fnorm;
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }