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: /*@C
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
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, (PetscErrorCode(*)(void **))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: /* ----------------------------------------------------------------------------- */
349: /*@C
350: TSPseudoSetVerifyTimeStep - Sets a user-defined routine to verify the quality of the
351: last timestep.
353: Logically Collective
355: Input Parameters:
356: + ts - timestep context
357: . dt - user-defined function to verify timestep
358: - ctx - [optional] user-defined context for private data for the timestep verification routine (may be `NULL`)
360: Calling sequence of `func`:
361: + ts - the time-step context
362: . update - latest solution vector
363: . ctx - [optional] user-defined timestep context
364: . newdt - the timestep to use for the next step
365: - flag - flag indicating whether the last time step was acceptable
367: Level: advanced
369: Note:
370: The routine set here will be called by `TSPseudoVerifyTimeStep()`
371: during the timestepping process.
373: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoVerifyTimeStepDefault()`, `TSPseudoVerifyTimeStep()`
374: @*/
375: PetscErrorCode TSPseudoSetVerifyTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, Vec update, void *ctx, PetscReal *newdt, PetscBool *flag), void *ctx)
376: {
377: PetscFunctionBegin;
379: PetscTryMethod(ts, "TSPseudoSetVerifyTimeStep_C", (TS, PetscErrorCode(*)(TS, Vec, void *, PetscReal *, PetscBool *), void *), (ts, dt, ctx));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@
384: TSPseudoSetTimeStepIncrement - Sets the scaling increment applied to
385: dt when using the TSPseudoTimeStepDefault() routine.
387: Logically Collective
389: Input Parameters:
390: + ts - the timestep context
391: - inc - the scaling factor >= 1.0
393: Options Database Key:
394: . -ts_pseudo_increment <increment> - set pseudo increment
396: Level: advanced
398: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
399: @*/
400: PetscErrorCode TSPseudoSetTimeStepIncrement(TS ts, PetscReal inc)
401: {
402: PetscFunctionBegin;
405: PetscTryMethod(ts, "TSPseudoSetTimeStepIncrement_C", (TS, PetscReal), (ts, inc));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*@
410: TSPseudoSetMaxTimeStep - Sets the maximum time step
411: when using the TSPseudoTimeStepDefault() routine.
413: Logically Collective
415: Input Parameters:
416: + ts - the timestep context
417: - maxdt - the maximum time step, use a non-positive value to deactivate
419: Options Database Key:
420: . -ts_pseudo_max_dt <increment> - set pseudo max dt
422: Level: advanced
424: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
425: @*/
426: PetscErrorCode TSPseudoSetMaxTimeStep(TS ts, PetscReal maxdt)
427: {
428: PetscFunctionBegin;
431: PetscTryMethod(ts, "TSPseudoSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: /*@
436: TSPseudoIncrementDtFromInitialDt - Indicates that a new timestep
437: is computed via the formula
438: $ dt = initial_dt*initial_fnorm/current_fnorm
439: rather than the default update,
440: $ dt = current_dt*previous_fnorm/current_fnorm.
442: Logically Collective
444: Input Parameter:
445: . ts - the timestep context
447: Options Database Key:
448: . -ts_pseudo_increment_dt_from_initial_dt <true,false> - use the initial dt to determine increment
450: Level: advanced
452: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoSetTimeStep()`, `TSPseudoTimeStepDefault()`
453: @*/
454: PetscErrorCode TSPseudoIncrementDtFromInitialDt(TS ts)
455: {
456: PetscFunctionBegin;
458: PetscTryMethod(ts, "TSPseudoIncrementDtFromInitialDt_C", (TS), (ts));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@C
463: TSPseudoSetTimeStep - Sets the user-defined routine to be
464: called at each pseudo-timestep to update the timestep.
466: Logically Collective
468: Input Parameters:
469: + ts - timestep context
470: . dt - function to compute timestep
471: - ctx - [optional] user-defined context for private data required by the function (may be `NULL`)
473: Calling sequence of `dt`:
474: + ts - the `TS` context
475: . newdt - the newly computed timestep
476: - ctx - [optional] user-defined context
478: Level: intermediate
480: Notes:
481: The routine set here will be called by `TSPseudoComputeTimeStep()`
482: during the timestepping process.
484: If not set then `TSPseudoTimeStepDefault()` is automatically used
486: .seealso: [](ch_ts), `TSPSEUDO`, `TSPseudoTimeStepDefault()`, `TSPseudoComputeTimeStep()`
487: @*/
488: PetscErrorCode TSPseudoSetTimeStep(TS ts, PetscErrorCode (*dt)(TS ts, PetscReal *newdt, void *ctx), void *ctx)
489: {
490: PetscFunctionBegin;
492: PetscTryMethod(ts, "TSPseudoSetTimeStep_C", (TS, PetscErrorCode(*)(TS, PetscReal *, void *), void *), (ts, dt, ctx));
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: /* ----------------------------------------------------------------------------- */
498: typedef PetscErrorCode (*FCN1)(TS, Vec, void *, PetscReal *, PetscBool *); /* force argument to next function to not be extern C*/
499: static PetscErrorCode TSPseudoSetVerifyTimeStep_Pseudo(TS ts, FCN1 dt, void *ctx)
500: {
501: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
503: PetscFunctionBegin;
504: pseudo->verify = dt;
505: pseudo->verifyctx = ctx;
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: static PetscErrorCode TSPseudoSetTimeStepIncrement_Pseudo(TS ts, PetscReal inc)
510: {
511: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
513: PetscFunctionBegin;
514: pseudo->dt_increment = inc;
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode TSPseudoSetMaxTimeStep_Pseudo(TS ts, PetscReal maxdt)
519: {
520: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
522: PetscFunctionBegin;
523: pseudo->dt_max = maxdt;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode TSPseudoIncrementDtFromInitialDt_Pseudo(TS ts)
528: {
529: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
531: PetscFunctionBegin;
532: pseudo->increment_dt_from_initial_dt = PETSC_TRUE;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: typedef PetscErrorCode (*FCN2)(TS, PetscReal *, void *); /* force argument to next function to not be extern C*/
537: static PetscErrorCode TSPseudoSetTimeStep_Pseudo(TS ts, FCN2 dt, void *ctx)
538: {
539: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
541: PetscFunctionBegin;
542: pseudo->dt = dt;
543: pseudo->dtctx = ctx;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*MC
548: TSPSEUDO - Solve steady state ODE and DAE problems with pseudo time stepping {cite}`ckk02` {cite}`kk97`
550: This method solves equations of the form
552: $$
553: F(X,Xdot) = 0
554: $$
556: for steady state using the iteration
558: $$
559: [G'] S = -F(X,0)
560: X += S
561: $$
563: where
565: $$
566: G(Y) = F(Y,(Y-X)/dt)
567: $$
569: This is linearly-implicit Euler with the residual always evaluated "at steady
570: state". See note below.
572: Options Database Keys:
573: + -ts_pseudo_increment <real> - ratio of increase dt
574: . -ts_pseudo_increment_dt_from_initial_dt <truth> - Increase dt as a ratio from original dt
575: . -ts_pseudo_fatol <atol> - stop iterating when the function norm is less than atol
576: - -ts_pseudo_frtol <rtol> - stop iterating when the function norm divided by the initial function norm is less than rtol
578: Level: beginner
580: Notes:
581: The residual computed by this method includes the transient term (Xdot is computed instead of
582: always being zero), but since the prediction from the last step is always the solution from the
583: last step, on the first Newton iteration we have
585: $$
586: Xdot = (Xpredicted - Xold)/dt = (Xold-Xold)/dt = 0
587: $$
589: Therefore, the linear system solved by the first Newton iteration is equivalent to the one
590: described above and in the papers. If the user chooses to perform multiple Newton iterations, the
591: algorithm is no longer the one described in the referenced papers.
593: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
594: M*/
595: PETSC_EXTERN PetscErrorCode TSCreate_Pseudo(TS ts)
596: {
597: TS_Pseudo *pseudo;
598: SNES snes;
599: SNESType stype;
601: PetscFunctionBegin;
602: ts->ops->reset = TSReset_Pseudo;
603: ts->ops->destroy = TSDestroy_Pseudo;
604: ts->ops->view = TSView_Pseudo;
605: ts->ops->setup = TSSetUp_Pseudo;
606: ts->ops->step = TSStep_Pseudo;
607: ts->ops->setfromoptions = TSSetFromOptions_Pseudo;
608: ts->ops->snesfunction = SNESTSFormFunction_Pseudo;
609: ts->ops->snesjacobian = SNESTSFormJacobian_Pseudo;
610: ts->default_adapt_type = TSADAPTNONE;
611: ts->usessnes = PETSC_TRUE;
613: PetscCall(TSGetSNES(ts, &snes));
614: PetscCall(SNESGetType(snes, &stype));
615: if (!stype) PetscCall(SNESSetType(snes, SNESKSPONLY));
617: PetscCall(PetscNew(&pseudo));
618: ts->data = (void *)pseudo;
620: pseudo->dt = TSPseudoTimeStepDefault;
621: pseudo->dtctx = NULL;
622: pseudo->dt_increment = 1.1;
623: pseudo->increment_dt_from_initial_dt = PETSC_FALSE;
624: pseudo->fnorm = -1;
625: pseudo->fnorm_initial = -1;
626: pseudo->fnorm_previous = -1;
627: #if defined(PETSC_USE_REAL_SINGLE)
628: pseudo->fatol = 1.e-25;
629: pseudo->frtol = 1.e-5;
630: #else
631: pseudo->fatol = 1.e-50;
632: pseudo->frtol = 1.e-12;
633: #endif
634: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetVerifyTimeStep_C", TSPseudoSetVerifyTimeStep_Pseudo));
635: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStepIncrement_C", TSPseudoSetTimeStepIncrement_Pseudo));
636: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetMaxTimeStep_C", TSPseudoSetMaxTimeStep_Pseudo));
637: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoIncrementDtFromInitialDt_C", TSPseudoIncrementDtFromInitialDt_Pseudo));
638: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSPseudoSetTimeStep_C", TSPseudoSetTimeStep_Pseudo));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /*@C
643: TSPseudoTimeStepDefault - Default code to compute pseudo-timestepping. Use with `TSPseudoSetTimeStep()`.
645: Collective
647: Input Parameters:
648: + ts - the timestep context
649: - dtctx - unused timestep context
651: Output Parameter:
652: . newdt - the timestep to use for the next step
654: Level: advanced
656: .seealso: [](ch_ts), `TSPseudoSetTimeStep()`, `TSPseudoComputeTimeStep()`, `TSPSEUDO`
657: @*/
658: PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
659: {
660: TS_Pseudo *pseudo = (TS_Pseudo *)ts->data;
661: PetscReal inc = pseudo->dt_increment;
663: PetscFunctionBegin;
664: PetscCall(VecZeroEntries(pseudo->xdot));
665: PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
666: PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
667: if (pseudo->fnorm_initial < 0) {
668: /* first time through so compute initial function norm */
669: pseudo->fnorm_initial = pseudo->fnorm;
670: pseudo->fnorm_previous = pseudo->fnorm;
671: }
672: if (pseudo->fnorm == 0.0) *newdt = 1.e12 * inc * ts->time_step;
673: else if (pseudo->increment_dt_from_initial_dt) *newdt = inc * pseudo->dt_initial * pseudo->fnorm_initial / pseudo->fnorm;
674: else *newdt = inc * ts->time_step * pseudo->fnorm_previous / pseudo->fnorm;
675: if (pseudo->dt_max > 0) *newdt = PetscMin(*newdt, pseudo->dt_max);
676: pseudo->fnorm_previous = pseudo->fnorm;
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }