Actual source code: eimex.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdm.h>
4: static const PetscInt TSEIMEXDefault = 3;
6: typedef struct {
7: PetscInt row_ind; /* Return the term T[row_ind][col_ind] */
8: PetscInt col_ind; /* Return the term T[row_ind][col_ind] */
9: PetscInt nstages; /* Numbers of stages in current scheme */
10: PetscInt max_rows; /* Maximum number of rows */
11: PetscInt *N; /* Harmonic sequence N[max_rows] */
12: Vec Y; /* States computed during the step, used to complete the step */
13: Vec Z; /* For shift*(Y-Z) */
14: Vec *T; /* Working table, size determined by nstages */
15: Vec YdotRHS; /* f(x) Work vector holding YdotRHS during residual evaluation */
16: Vec YdotI; /* xdot-g(x) Work vector holding YdotI = G(t,x,xdot) when xdot =0 */
17: Vec Ydot; /* f(x)+g(x) Work vector */
18: Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation) */
19: PetscReal shift;
20: PetscReal ctime;
21: PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
22: PetscBool ord_adapt; /* order adapativity */
23: TSStepStatus status;
24: } TS_EIMEX;
26: /* This function is pure */
27: static PetscInt Map(PetscInt i, PetscInt j, PetscInt s)
28: {
29: return ((2 * s - j + 1) * j / 2 + i - j);
30: }
32: static PetscErrorCode TSEvaluateStep_EIMEX(TS ts, PetscInt order, Vec X, PetscBool *done)
33: {
34: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
35: const PetscInt ns = ext->nstages;
36: PetscFunctionBegin;
37: PetscCall(VecCopy(ext->T[Map(ext->row_ind, ext->col_ind, ns)], X));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode TSStage_EIMEX(TS ts, PetscInt istage)
42: {
43: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
44: PetscReal h;
45: Vec Y = ext->Y, Z = ext->Z;
46: SNES snes;
47: TSAdapt adapt;
48: PetscInt i, its, lits;
49: PetscBool accept;
51: PetscFunctionBegin;
52: PetscCall(TSGetSNES(ts, &snes));
53: h = ts->time_step / ext->N[istage]; /* step size for the istage-th stage */
54: ext->shift = 1. / h;
55: PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again */
56: PetscCall(VecCopy(ext->VecSolPrev, Y)); /* Take the previous solution as initial step */
58: for (i = 0; i < ext->N[istage]; i++) {
59: ext->ctime = ts->ptime + h * i;
60: PetscCall(VecCopy(Y, Z)); /* Save the solution of the previous substep */
61: PetscCall(SNESSolve(snes, NULL, Y));
62: PetscCall(SNESGetIterationNumber(snes, &its));
63: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
64: ts->snes_its += its;
65: ts->ksp_its += lits;
66: PetscCall(TSGetAdapt(ts, &adapt));
67: PetscCall(TSAdaptCheckStage(adapt, ts, ext->ctime, Y, &accept));
68: }
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: static PetscErrorCode TSStep_EIMEX(TS ts)
73: {
74: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
75: const PetscInt ns = ext->nstages;
76: Vec *T = ext->T, Y = ext->Y;
77: SNES snes;
78: PetscInt i, j;
79: PetscBool accept = PETSC_FALSE;
80: PetscReal alpha, local_error, local_error_a, local_error_r;
82: PetscFunctionBegin;
83: PetscCall(TSGetSNES(ts, &snes));
84: PetscCall(SNESSetType(snes, "ksponly"));
85: ext->status = TS_STEP_INCOMPLETE;
87: PetscCall(VecCopy(ts->vec_sol, ext->VecSolPrev));
89: /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
90: for (j = 0; j < ns; j++) {
91: PetscCall(TSStage_EIMEX(ts, j));
92: PetscCall(VecCopy(Y, T[j]));
93: }
95: for (i = 1; i < ns; i++) {
96: for (j = i; j < ns; j++) {
97: alpha = -(PetscReal)ext->N[j] / ext->N[j - i];
98: PetscCall(VecAXPBYPCZ(T[Map(j, i, ns)], alpha, 1.0, 0, T[Map(j, i - 1, ns)], T[Map(j - 1, i - 1, ns)])); /* T[j][i]=alpha*T[j][i-1]+T[j-1][i-1] */
99: alpha = 1.0 / (1.0 + alpha);
100: PetscCall(VecScale(T[Map(j, i, ns)], alpha));
101: }
102: }
104: PetscCall(TSEvaluateStep(ts, ns, ts->vec_sol, NULL)); /*update ts solution */
106: if (ext->ord_adapt && ext->nstages < ext->max_rows) {
107: accept = PETSC_FALSE;
108: while (!accept && ext->nstages < ext->max_rows) {
109: PetscCall(TSErrorWeightedNorm(ts, ts->vec_sol, T[Map(ext->nstages - 1, ext->nstages - 2, ext->nstages)], ts->adapt->wnormtype, &local_error, &local_error_a, &local_error_r));
110: accept = (local_error < 1.0) ? PETSC_TRUE : PETSC_FALSE;
112: if (!accept) { /* add one more stage*/
113: PetscCall(TSStage_EIMEX(ts, ext->nstages));
114: ext->nstages++;
115: ext->row_ind++;
116: ext->col_ind++;
117: /*T table need to be recycled*/
118: PetscCall(VecDuplicateVecs(ts->vec_sol, (1 + ext->nstages) * ext->nstages / 2, &ext->T));
119: for (i = 0; i < ext->nstages - 1; i++) {
120: for (j = 0; j <= i; j++) PetscCall(VecCopy(T[Map(i, j, ext->nstages - 1)], ext->T[Map(i, j, ext->nstages)]));
121: }
122: PetscCall(VecDestroyVecs(ext->nstages * (ext->nstages - 1) / 2, &T));
123: T = ext->T; /*reset the pointer*/
124: /*recycling finished, store the new solution*/
125: PetscCall(VecCopy(Y, T[ext->nstages - 1]));
126: /*extrapolation for the newly added stage*/
127: for (i = 1; i < ext->nstages; i++) {
128: alpha = -(PetscReal)ext->N[ext->nstages - 1] / ext->N[ext->nstages - 1 - i];
129: PetscCall(VecAXPBYPCZ(T[Map(ext->nstages - 1, i, ext->nstages)], alpha, 1.0, 0, T[Map(ext->nstages - 1, i - 1, ext->nstages)], T[Map(ext->nstages - 1 - 1, i - 1, ext->nstages)])); /*T[ext->nstages-1][i]=alpha*T[ext->nstages-1][i-1]+T[ext->nstages-1-1][i-1]*/
130: alpha = 1.0 / (1.0 + alpha);
131: PetscCall(VecScale(T[Map(ext->nstages - 1, i, ext->nstages)], alpha));
132: }
133: /*update ts solution */
134: PetscCall(TSEvaluateStep(ts, ext->nstages, ts->vec_sol, NULL));
135: } /*end if !accept*/
136: } /*end while*/
138: if (ext->nstages == ext->max_rows) PetscCall(PetscInfo(ts, "Max number of rows has been used\n"));
139: } /*end if ext->ord_adapt*/
140: ts->ptime += ts->time_step;
141: ext->status = TS_STEP_COMPLETE;
143: if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /* cubic Hermit spline */
148: static PetscErrorCode TSInterpolate_EIMEX(TS ts, PetscReal itime, Vec X)
149: {
150: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
151: PetscReal t, a, b;
152: Vec Y0 = ext->VecSolPrev, Y1 = ext->Y, Ydot = ext->Ydot, YdotI = ext->YdotI;
153: const PetscReal h = ts->ptime - ts->ptime_prev;
154: PetscFunctionBegin;
155: t = (itime - ts->ptime + h) / h;
156: /* YdotI = -f(x)-g(x) */
158: PetscCall(VecZeroEntries(Ydot));
159: PetscCall(TSComputeIFunction(ts, ts->ptime - h, Y0, Ydot, YdotI, PETSC_FALSE));
161: a = 2.0 * t * t * t - 3.0 * t * t + 1.0;
162: b = -(t * t * t - 2.0 * t * t + t) * h;
163: PetscCall(VecAXPBYPCZ(X, a, b, 0.0, Y0, YdotI));
165: PetscCall(TSComputeIFunction(ts, ts->ptime, Y1, Ydot, YdotI, PETSC_FALSE));
166: a = -2.0 * t * t * t + 3.0 * t * t;
167: b = -(t * t * t - t * t) * h;
168: PetscCall(VecAXPBYPCZ(X, a, b, 1.0, Y1, YdotI));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode TSReset_EIMEX(TS ts)
174: {
175: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
176: PetscInt ns;
178: PetscFunctionBegin;
179: ns = ext->nstages;
180: PetscCall(VecDestroyVecs((1 + ns) * ns / 2, &ext->T));
181: PetscCall(VecDestroy(&ext->Y));
182: PetscCall(VecDestroy(&ext->Z));
183: PetscCall(VecDestroy(&ext->YdotRHS));
184: PetscCall(VecDestroy(&ext->YdotI));
185: PetscCall(VecDestroy(&ext->Ydot));
186: PetscCall(VecDestroy(&ext->VecSolPrev));
187: PetscCall(PetscFree(ext->N));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode TSDestroy_EIMEX(TS ts)
192: {
193: PetscFunctionBegin;
194: PetscCall(TSReset_EIMEX(ts));
195: PetscCall(PetscFree(ts->data));
196: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetMaxRows_C", NULL));
197: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetRowCol_C", NULL));
198: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetOrdAdapt_C", NULL));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode TSEIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot, Vec *YdotI, Vec *YdotRHS)
203: {
204: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
206: PetscFunctionBegin;
207: if (Z) {
208: if (dm && dm != ts->dm) {
209: PetscCall(DMGetNamedGlobalVector(dm, "TSEIMEX_Z", Z));
210: } else *Z = ext->Z;
211: }
212: if (Ydot) {
213: if (dm && dm != ts->dm) {
214: PetscCall(DMGetNamedGlobalVector(dm, "TSEIMEX_Ydot", Ydot));
215: } else *Ydot = ext->Ydot;
216: }
217: if (YdotI) {
218: if (dm && dm != ts->dm) {
219: PetscCall(DMGetNamedGlobalVector(dm, "TSEIMEX_YdotI", YdotI));
220: } else *YdotI = ext->YdotI;
221: }
222: if (YdotRHS) {
223: if (dm && dm != ts->dm) {
224: PetscCall(DMGetNamedGlobalVector(dm, "TSEIMEX_YdotRHS", YdotRHS));
225: } else *YdotRHS = ext->YdotRHS;
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode TSEIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot, Vec *YdotI, Vec *YdotRHS)
231: {
232: PetscFunctionBegin;
233: if (Z) {
234: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSEIMEX_Z", Z));
235: }
236: if (Ydot) {
237: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSEIMEX_Ydot", Ydot));
238: }
239: if (YdotI) {
240: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSEIMEX_YdotI", YdotI));
241: }
242: if (YdotRHS) {
243: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSEIMEX_YdotRHS", YdotRHS));
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*
249: This defines the nonlinear equation that is to be solved with SNES
250: Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
251: In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
252: Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
253: */
254: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes, Vec X, Vec G, TS ts)
255: {
256: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
257: Vec Ydot, Z;
258: DM dm, dmsave;
260: PetscFunctionBegin;
261: PetscCall(VecZeroEntries(G));
263: PetscCall(SNESGetDM(snes, &dm));
264: PetscCall(TSEIMEXGetVecs(ts, dm, &Z, &Ydot, NULL, NULL));
265: PetscCall(VecZeroEntries(Ydot));
266: dmsave = ts->dm;
267: ts->dm = dm;
268: PetscCall(TSComputeIFunction(ts, ext->ctime, X, Ydot, G, PETSC_FALSE));
269: /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function. */
270: PetscCall(VecCopy(G, Ydot));
271: ts->dm = dmsave;
272: PetscCall(TSEIMEXRestoreVecs(ts, dm, &Z, &Ydot, NULL, NULL));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*
278: This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
279: */
280: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts)
281: {
282: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
283: Vec Ydot;
284: DM dm, dmsave;
285: PetscFunctionBegin;
286: PetscCall(SNESGetDM(snes, &dm));
287: PetscCall(TSEIMEXGetVecs(ts, dm, NULL, &Ydot, NULL, NULL));
288: /* PetscCall(VecZeroEntries(Ydot)); */
289: /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
290: dmsave = ts->dm;
291: ts->dm = dm;
292: PetscCall(TSComputeIJacobian(ts, ts->ptime, X, Ydot, ext->shift, A, B, PETSC_TRUE));
293: ts->dm = dmsave;
294: PetscCall(TSEIMEXRestoreVecs(ts, dm, NULL, &Ydot, NULL, NULL));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine, DM coarse, void *ctx)
299: {
300: PetscFunctionBegin;
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
305: {
306: TS ts = (TS)ctx;
307: Vec Z, Z_c;
309: PetscFunctionBegin;
310: PetscCall(TSEIMEXGetVecs(ts, fine, &Z, NULL, NULL, NULL));
311: PetscCall(TSEIMEXGetVecs(ts, coarse, &Z_c, NULL, NULL, NULL));
312: PetscCall(MatRestrict(restrct, Z, Z_c));
313: PetscCall(VecPointwiseMult(Z_c, rscale, Z_c));
314: PetscCall(TSEIMEXRestoreVecs(ts, fine, &Z, NULL, NULL, NULL));
315: PetscCall(TSEIMEXRestoreVecs(ts, coarse, &Z_c, NULL, NULL, NULL));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: static PetscErrorCode TSSetUp_EIMEX(TS ts)
320: {
321: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
322: DM dm;
324: PetscFunctionBegin;
325: if (!ext->N) { /* ext->max_rows not set */
326: PetscCall(TSEIMEXSetMaxRows(ts, TSEIMEXDefault));
327: }
328: if (-1 == ext->row_ind && -1 == ext->col_ind) {
329: PetscCall(TSEIMEXSetRowCol(ts, ext->max_rows, ext->max_rows));
330: } else { /* ext->row_ind and col_ind already set */
331: if (ext->ord_adapt) PetscCall(PetscInfo(ts, "Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n"));
332: }
334: if (ext->ord_adapt) {
335: ext->nstages = 2; /* Start with the 2-stage scheme */
336: PetscCall(TSEIMEXSetRowCol(ts, ext->nstages, ext->nstages));
337: } else {
338: ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
339: }
341: PetscCall(TSGetAdapt(ts, &ts->adapt));
343: PetscCall(VecDuplicateVecs(ts->vec_sol, (1 + ext->nstages) * ext->nstages / 2, &ext->T)); /* full T table */
344: PetscCall(VecDuplicate(ts->vec_sol, &ext->YdotI));
345: PetscCall(VecDuplicate(ts->vec_sol, &ext->YdotRHS));
346: PetscCall(VecDuplicate(ts->vec_sol, &ext->Ydot));
347: PetscCall(VecDuplicate(ts->vec_sol, &ext->VecSolPrev));
348: PetscCall(VecDuplicate(ts->vec_sol, &ext->Y));
349: PetscCall(VecDuplicate(ts->vec_sol, &ext->Z));
350: PetscCall(TSGetDM(ts, &dm));
351: if (dm) PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSEIMEX, DMRestrictHook_TSEIMEX, ts));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode TSSetFromOptions_EIMEX(TS ts, PetscOptionItems *PetscOptionsObject)
356: {
357: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
358: PetscInt tindex[2];
359: PetscInt np = 2, nrows = TSEIMEXDefault;
361: PetscFunctionBegin;
362: tindex[0] = TSEIMEXDefault;
363: tindex[1] = TSEIMEXDefault;
364: PetscOptionsHeadBegin(PetscOptionsObject, "EIMEX ODE solver options");
365: {
366: PetscBool flg;
367: PetscCall(PetscOptionsInt("-ts_eimex_max_rows", "Define the maximum number of rows used", "TSEIMEXSetMaxRows", nrows, &nrows, &flg)); /* default value 3 */
368: if (flg) PetscCall(TSEIMEXSetMaxRows(ts, nrows));
369: PetscCall(PetscOptionsIntArray("-ts_eimex_row_col", "Return the specific term in the T table", "TSEIMEXSetRowCol", tindex, &np, &flg));
370: if (flg) PetscCall(TSEIMEXSetRowCol(ts, tindex[0], tindex[1]));
371: PetscCall(PetscOptionsBool("-ts_eimex_order_adapt", "Solve the problem with adaptive order", "TSEIMEXSetOrdAdapt", ext->ord_adapt, &ext->ord_adapt, NULL));
372: }
373: PetscOptionsHeadEnd();
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: static PetscErrorCode TSView_EIMEX(TS ts, PetscViewer viewer)
378: {
379: PetscFunctionBegin;
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@C
384: TSEIMEXSetMaxRows - Set the maximum number of rows for `TSEIMEX` schemes
386: Logically Collective
388: Input Parameters:
389: + ts - timestepping context
390: - nrows - maximum number of rows
392: Level: intermediate
394: .seealso: [](ch_ts), `TSEIMEXSetRowCol()`, `TSEIMEXSetOrdAdapt()`, `TSEIMEX`
395: @*/
396: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
397: {
398: PetscFunctionBegin;
400: PetscTryMethod(ts, "TSEIMEXSetMaxRows_C", (TS, PetscInt), (ts, nrows));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@C
405: TSEIMEXSetRowCol - Set the number of rows and the number of columns for the tableau that represents the T solution in the `TSEIMEX` scheme
407: Logically Collective
409: Input Parameters:
410: + ts - timestepping context
411: . row - the row
412: - col - the column
414: Level: intermediate
416: .seealso: [](ch_ts), `TSEIMEXSetMaxRows()`, `TSEIMEXSetOrdAdapt()`, `TSEIMEX`
417: @*/
418: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
419: {
420: PetscFunctionBegin;
422: PetscTryMethod(ts, "TSEIMEXSetRowCol_C", (TS, PetscInt, PetscInt), (ts, row, col));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@C
427: TSEIMEXSetOrdAdapt - Set the order adaptativity for the `TSEIMEX` schemes
429: Logically Collective
431: Input Parameters:
432: + ts - timestepping context
433: - flg - index in the T table
435: Level: intermediate
437: .seealso: [](ch_ts), `TSEIMEXSetRowCol()`, `TSEIMEX`
438: @*/
439: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
440: {
441: PetscFunctionBegin;
443: PetscTryMethod(ts, "TSEIMEXSetOrdAdapt_C", (TS, PetscBool), (ts, flg));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts, PetscInt nrows)
448: {
449: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
450: PetscInt i;
452: PetscFunctionBegin;
453: PetscCheck(nrows >= 0 && nrows <= 100, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Max number of rows (current value %" PetscInt_FMT ") should be an integer number between 1 and 100", nrows);
454: PetscCall(PetscFree(ext->N));
455: ext->max_rows = nrows;
456: PetscCall(PetscMalloc1(nrows, &ext->N));
457: for (i = 0; i < nrows; i++) ext->N[i] = i + 1;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts, PetscInt row, PetscInt col)
462: {
463: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
465: PetscFunctionBegin;
466: PetscCheck(row >= 1 && col >= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "The row or column index (current value %" PetscInt_FMT ",%" PetscInt_FMT ") should not be less than 1 ", row, col);
467: PetscCheck(row <= ext->max_rows && col <= ext->max_rows, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "The row or column index (current value %" PetscInt_FMT ",%" PetscInt_FMT ") exceeds the maximum number of rows %" PetscInt_FMT, row, col,
468: ext->max_rows);
469: PetscCheck(col <= row, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "The column index (%" PetscInt_FMT ") exceeds the row index (%" PetscInt_FMT ")", col, row);
471: ext->row_ind = row - 1;
472: ext->col_ind = col - 1; /* Array index in C starts from 0 */
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts, PetscBool flg)
477: {
478: TS_EIMEX *ext = (TS_EIMEX *)ts->data;
479: PetscFunctionBegin;
480: ext->ord_adapt = flg;
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*MC
485: TSEIMEX - Time stepping with Extrapolated IMEX methods {cite}`constantinescu_a2010a`.
487: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it
488: is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using `TSSetIFunction()` and the
489: non-stiff part with `TSSetRHSFunction()`.
491: Level: beginner
493: Notes:
494: The default is a 3-stage scheme, it can be changed with `TSEIMEXSetMaxRows()` or -ts_eimex_max_rows
496: This method currently only works with ODE, for which the stiff part $ G(t,X,Xdot) $ has the form $ Xdot + Ghat(t,X)$.
498: The general system is written as
500: $$
501: G(t,X,Xdot) = F(t,X)
502: $$
504: where G represents the stiff part and F represents the non-stiff part. The user should provide the stiff part
505: of the equation using TSSetIFunction() and the non-stiff part with `TSSetRHSFunction()`.
506: This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.
508: Another common form for the system is
510: $$
511: y'=f(x)+g(x)
512: $$
514: The relationship between F,G and f,g is
516: $$
517: G = y'-g(x), F = f(x)
518: $$
520: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSEIMEXSetMaxRows()`, `TSEIMEXSetRowCol()`, `TSEIMEXSetOrdAdapt()`, `TSType`
521: M*/
522: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
523: {
524: TS_EIMEX *ext;
526: PetscFunctionBegin;
528: ts->ops->reset = TSReset_EIMEX;
529: ts->ops->destroy = TSDestroy_EIMEX;
530: ts->ops->view = TSView_EIMEX;
531: ts->ops->setup = TSSetUp_EIMEX;
532: ts->ops->step = TSStep_EIMEX;
533: ts->ops->interpolate = TSInterpolate_EIMEX;
534: ts->ops->evaluatestep = TSEvaluateStep_EIMEX;
535: ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
536: ts->ops->snesfunction = SNESTSFormFunction_EIMEX;
537: ts->ops->snesjacobian = SNESTSFormJacobian_EIMEX;
538: ts->default_adapt_type = TSADAPTNONE;
540: ts->usessnes = PETSC_TRUE;
542: PetscCall(PetscNew(&ext));
543: ts->data = (void *)ext;
545: ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
546: ext->row_ind = -1;
547: ext->col_ind = -1;
548: ext->max_rows = TSEIMEXDefault;
549: ext->nstages = TSEIMEXDefault;
551: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX));
552: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetRowCol_C", TSEIMEXSetRowCol_EIMEX));
553: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetOrdAdapt_C", TSEIMEXSetOrdAdapt_EIMEX));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }