Actual source code: arkimex.c

  1: /*
  2:   Code for timestepping with additive Runge-Kutta IMEX method

  4:   Notes:
  5:   The general system is written as

  7:   F(t,U,Udot) = G(t,U)

  9:   where F represents the stiff part of the physics and G represents the non-stiff part.

 11: */
 12: #include <petsc/private/tsimpl.h>
 13: #include <petscdm.h>

 15: static TSARKIMEXType  TSARKIMEXDefault = TSARKIMEX3;
 16: static PetscBool      TSARKIMEXRegisterAllCalled;
 17: static PetscBool      TSARKIMEXPackageInitialized;
 18: static PetscErrorCode TSExtrapolate_ARKIMEX(TS, PetscReal, Vec);

 20: typedef struct _ARKTableau *ARKTableau;
 21: struct _ARKTableau {
 22:   char      *name;
 23:   PetscInt   order;                /* Classical approximation order of the method */
 24:   PetscInt   s;                    /* Number of stages */
 25:   PetscBool  stiffly_accurate;     /* The implicit part is stiffly accurate*/
 26:   PetscBool  FSAL_implicit;        /* The implicit part is FSAL*/
 27:   PetscBool  explicit_first_stage; /* The implicit part has an explicit first stage*/
 28:   PetscInt   pinterp;              /* Interpolation order */
 29:   PetscReal *At, *bt, *ct;         /* Stiff tableau */
 30:   PetscReal *A, *b, *c;            /* Non-stiff tableau */
 31:   PetscReal *bembedt, *bembed;     /* Embedded formula of order one less (order-1) */
 32:   PetscReal *binterpt, *binterp;   /* Dense output formula */
 33:   PetscReal  ccfl;                 /* Placeholder for CFL coefficient relative to forward Euler */
 34: };
 35: typedef struct _ARKTableauLink *ARKTableauLink;
 36: struct _ARKTableauLink {
 37:   struct _ARKTableau tab;
 38:   ARKTableauLink     next;
 39: };
 40: static ARKTableauLink ARKTableauList;

 42: typedef struct {
 43:   ARKTableau   tableau;
 44:   Vec         *Y;            /* States computed during the step */
 45:   Vec         *YdotI;        /* Time derivatives for the stiff part */
 46:   Vec         *YdotRHS;      /* Function evaluations for the non-stiff part */
 47:   Vec         *Y_prev;       /* States computed during the previous time step */
 48:   Vec         *YdotI_prev;   /* Time derivatives for the stiff part for the previous time step*/
 49:   Vec         *YdotRHS_prev; /* Function evaluations for the non-stiff part for the previous time step*/
 50:   Vec          Ydot0;        /* Holds the slope from the previous step in FSAL case */
 51:   Vec          Ydot;         /* Work vector holding Ydot during residual evaluation */
 52:   Vec          Z;            /* Ydot = shift(Y-Z) */
 53:   PetscScalar *work;         /* Scalar work */
 54:   PetscReal    scoeff;       /* shift = scoeff/dt */
 55:   PetscReal    stage_time;
 56:   PetscBool    imex;
 57:   PetscBool    extrapolate; /* Extrapolate initial guess from previous time-step stage values */
 58:   TSStepStatus status;

 60:   /* context for sensitivity analysis */
 61:   Vec *VecsDeltaLam;   /* Increment of the adjoint sensitivity w.r.t IC at stage */
 62:   Vec *VecsSensiTemp;  /* Vectors to be multiplied with Jacobian transpose */
 63:   Vec *VecsSensiPTemp; /* Temporary Vectors to store JacobianP-transpose-vector product */
 64: } TS_ARKIMEX;
 65: /*MC
 66:      TSARKIMEXARS122 - Second order ARK IMEX scheme.

 68:      This method has one explicit stage and one implicit stage.

 70:      Options Database Key:
 71: .      -ts_arkimex_type ars122 - set arkimex type to ars122

 73:      Level: advanced

 75:      References:
 76: .    * -  U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997).

 78: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
 79: M*/
 80: /*MC
 81:      TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part.

 83:      This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu.

 85:      Options Database Key:
 86: .      -ts_arkimex_type a2 - set arkimex type to a2

 88:      Level: advanced

 90: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
 91: M*/
 92: /*MC
 93:      TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part.

 95:      This method has two implicit stages, and L-stable implicit scheme.

 97:      Options Database Key:
 98: .      -ts_arkimex_type l2 - set arkimex type to l2

100:      Level: advanced

102:     References:
103: .   * -  L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005.

105: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
106: M*/
107: /*MC
108:      TSARKIMEX1BEE - First order backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method.

110:      This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used.

112:      Options Database Key:
113: .      -ts_arkimex_type 1bee - set arkimex type to 1bee

115:      Level: advanced

117: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
118: M*/
119: /*MC
120:      TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part.

122:      This method has one explicit stage and two implicit stages. The implicit part is the same as in TSARKIMEX2D and TSARKIMEX2E, but the explicit part has a larger stability region on the negative real axis. This method was provided by Emil Constantinescu.

124:      Options Database Key:
125: .      -ts_arkimex_type 2c - set arkimex type to 2c

127:      Level: advanced

129: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
130: M*/
131: /*MC
132:      TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part.

134:      This method has one explicit stage and two implicit stages. The stability function is independent of the explicit part in the infinity limit of the implicit component. This method was provided by Emil Constantinescu.

136:      Options Database Key:
137: .      -ts_arkimex_type 2d - set arkimex type to 2d

139:      Level: advanced

141: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
142: M*/
143: /*MC
144:      TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part.

146:      This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu.

148:      Options Database Key:
149: .      -ts_arkimex_type 2e - set arkimex type to 2e

151:     Level: advanced

153: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
154: M*/
155: /*MC
156:      TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme.

158:      This method has three implicit stages.

160:      References:
161: .    * -  L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005.

163:      This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375

165:      Options Database Key:
166: .      -ts_arkimex_type prssp2 - set arkimex type to prssp2

168:      Level: advanced

170: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
171: M*/
172: /*MC
173:      TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part.

175:      This method has one explicit stage and three implicit stages.

177:      Options Database Key:
178: .      -ts_arkimex_type 3 - set arkimex type to 3

180:      Level: advanced

182:      References:
183: .    * -  Kennedy and Carpenter 2003.

185: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
186: M*/
187: /*MC
188:      TSARKIMEXARS443 - Third order ARK IMEX scheme.

190:      This method has one explicit stage and four implicit stages.

192:      Options Database Key:
193: .      -ts_arkimex_type ars443 - set arkimex type to ars443

195:      Level: advanced

197:      References:
198: +    * -  U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997).
199: -    * -  This method is referred to as ARS(4,4,3) in https://arxiv.org/abs/1110.4375

201: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
202: M*/
203: /*MC
204:      TSARKIMEXBPR3 - Third order ARK IMEX scheme.

206:      This method has one explicit stage and four implicit stages.

208:      Options Database Key:
209: .      -ts_arkimex_type bpr3 - set arkimex type to bpr3

211:      Level: advanced

213:      References:
214: .    * - This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375

216: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
217: M*/
218: /*MC
219:      TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part.

221:      This method has one explicit stage and four implicit stages.

223:      Options Database Key:
224: .      -ts_arkimex_type 4 - set arkimex type to4

226:      Level: advanced

228:      References:
229: .    * - Kennedy and Carpenter 2003.

231: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
232: M*/
233: /*MC
234:      TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part.

236:      This method has one explicit stage and five implicit stages.

238:      Options Database Key:
239: .      -ts_arkimex_type 5 - set arkimex type to 5

241:      Level: advanced

243:      References:
244: .    * - Kennedy and Carpenter 2003.

246: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEXSetType()`
247: M*/

249: /*@C
250:   TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in `TSARKIMEX`

252:   Not Collective, but should be called by all processes which will need the schemes to be registered

254:   Level: advanced

256: .seealso: [](ch_ts), `TS`, `TSARKIMEX`, `TSARKIMEXRegisterDestroy()`
257: @*/
258: PetscErrorCode TSARKIMEXRegisterAll(void)
259: {
260:   PetscFunctionBegin;
261:   if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
262:   TSARKIMEXRegisterAllCalled = PETSC_TRUE;

264:   {
265:     const PetscReal A[3][3] =
266:       {
267:         {0.0, 0.0, 0.0},
268:         {0.0, 0.0, 0.0},
269:         {0.0, 0.5, 0.0}
270:     },
271:                     At[3][3] = {{1.0, 0.0, 0.0}, {0.0, 0.5, 0.0}, {0.0, 0.5, 0.5}}, b[3] = {0.0, 0.5, 0.5}, bembedt[3] = {1.0, 0.0, 0.0};
272:     PetscCall(TSARKIMEXRegister(TSARKIMEX1BEE, 2, 3, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
273:   }
274:   {
275:     const PetscReal A[2][2] =
276:       {
277:         {0.0, 0.0},
278:         {0.5, 0.0}
279:     },
280:                     At[2][2] = {{0.0, 0.0}, {0.0, 0.5}}, b[2] = {0.0, 1.0}, bembedt[2] = {0.5, 0.5};
281:     /* binterpt[2][2] = {{1.0,-1.0},{0.0,1.0}};  second order dense output has poor stability properties and hence it is not currently in use*/
282:     PetscCall(TSARKIMEXRegister(TSARKIMEXARS122, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
283:   }
284:   {
285:     const PetscReal A[2][2] =
286:       {
287:         {0.0, 0.0},
288:         {1.0, 0.0}
289:     },
290:                     At[2][2] = {{0.0, 0.0}, {0.5, 0.5}}, b[2] = {0.5, 0.5}, bembedt[2] = {0.0, 1.0};
291:     /* binterpt[2][2] = {{1.0,-0.5},{0.0,0.5}}  second order dense output has poor stability properties and hence it is not currently in use*/
292:     PetscCall(TSARKIMEXRegister(TSARKIMEXA2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 1, b, NULL));
293:   }
294:   {
295:     /* const PetscReal us2 = 1.0-1.0/PetscSqrtReal((PetscReal)2.0);    Direct evaluation: 0.2928932188134524755992. Used below to ensure all values are available at compile time   */
296:     const PetscReal A[2][2] =
297:       {
298:         {0.0, 0.0},
299:         {1.0, 0.0}
300:     },
301:                     At[2][2] = {{0.2928932188134524755992, 0.0}, {1.0 - 2.0 * 0.2928932188134524755992, 0.2928932188134524755992}}, b[2] = {0.5, 0.5}, bembedt[2] = {0.0, 1.0}, binterpt[2][2] = {{(0.2928932188134524755992 - 1.0) / (2.0 * 0.2928932188134524755992 - 1.0), -1 / (2.0 * (1.0 - 2.0 * 0.2928932188134524755992))}, {1 - (0.2928932188134524755992 - 1.0) / (2.0 * 0.2928932188134524755992 - 1.0), -1 / (2.0 * (1.0 - 2.0 * 0.2928932188134524755992))}}, binterp[2][2] = {{1.0, -0.5}, {0.0, 0.5}};
302:     PetscCall(TSARKIMEXRegister(TSARKIMEXL2, 2, 2, &At[0][0], b, NULL, &A[0][0], b, NULL, bembedt, bembedt, 2, binterpt[0], binterp[0]));
303:   }
304:   {
305:     /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0),  Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time   */
306:     const PetscReal A[3][3] =
307:       {
308:         {0,                           0,   0},
309:         {2 - 1.414213562373095048802, 0,   0},
310:         {0.5,                         0.5, 0}
311:     },
312:                     At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}};
313:     PetscCall(TSARKIMEXRegister(TSARKIMEX2C, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
314:   }
315:   {
316:     /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0),  Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time   */
317:     const PetscReal A[3][3] =
318:       {
319:         {0,                           0,    0},
320:         {2 - 1.414213562373095048802, 0,    0},
321:         {0.75,                        0.25, 0}
322:     },
323:                     At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}};
324:     PetscCall(TSARKIMEXRegister(TSARKIMEX2D, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
325:   }
326:   { /* Optimal for linear implicit part */
327:     /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0),  Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time   */
328:     const PetscReal A[3][3] =
329:       {
330:         {0,                                     0,                                     0},
331:         {2 - 1.414213562373095048802,           0,                                     0},
332:         {(3 - 2 * 1.414213562373095048802) / 6, (3 + 2 * 1.414213562373095048802) / 6, 0}
333:     },
334:                     At[3][3] = {{0, 0, 0}, {1 - 1 / 1.414213562373095048802, 1 - 1 / 1.414213562373095048802, 0}, {1 / (2 * 1.414213562373095048802), 1 / (2 * 1.414213562373095048802), 1 - 1 / 1.414213562373095048802}}, bembedt[3] = {(4. - 1.414213562373095048802) / 8., (4. - 1.414213562373095048802) / 8., 1 / (2. * 1.414213562373095048802)}, binterpt[3][2] = {{1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 / 1.414213562373095048802, -1.0 / (2.0 * 1.414213562373095048802)}, {1.0 - 1.414213562373095048802, 1.0 / 1.414213562373095048802}};
335:     PetscCall(TSARKIMEXRegister(TSARKIMEX2E, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
336:   }
337:   { /* Optimal for linear implicit part */
338:     const PetscReal A[3][3] =
339:       {
340:         {0,   0,   0},
341:         {0.5, 0,   0},
342:         {0.5, 0.5, 0}
343:     },
344:                     At[3][3] = {{0.25, 0, 0}, {0, 0.25, 0}, {1. / 3, 1. / 3, 1. / 3}};
345:     PetscCall(TSARKIMEXRegister(TSARKIMEXPRSSP2, 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, NULL, NULL, 0, NULL, NULL));
346:   }
347:   {
348:     const PetscReal A[4][4] =
349:       {
350:         {0,                                0,                                0,                                 0},
351:         {1767732205903. / 2027836641118.,  0,                                0,                                 0},
352:         {5535828885825. / 10492691773637., 788022342437. / 10882634858940.,  0,                                 0},
353:         {6485989280629. / 16251701735622., -4246266847089. / 9704473918619., 10755448449292. / 10357097424841., 0}
354:     },
355:                     At[4][4] = {{0, 0, 0, 0}, {1767732205903. / 4055673282236., 1767732205903. / 4055673282236., 0, 0}, {2746238789719. / 10658868560708., -640167445237. / 6845629431997., 1767732205903. / 4055673282236., 0}, {1471266399579. / 7840856788654., -4482444167858. / 7529755066697., 11266239266428. / 11593286722821., 1767732205903. / 4055673282236.}}, bembedt[4] = {2756255671327. / 12835298489170., -10771552573575. / 22201958757719., 9247589265047. / 10645013368117., 2193209047091. / 5459859503100.}, binterpt[4][2] = {{4655552711362. / 22874653954995., -215264564351. / 13552729205753.}, {-18682724506714. / 9892148508045., 17870216137069. / 13817060693119.}, {34259539580243. / 13192909600954., -28141676662227. / 17317692491321.}, {584795268549. / 6622622206610., 2508943948391. / 7218656332882.}};
356:     PetscCall(TSARKIMEXRegister(TSARKIMEX3, 3, 4, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 2, binterpt[0], NULL));
357:   }
358:   {
359:     const PetscReal A[5][5] =
360:       {
361:         {0,        0,       0,      0,       0},
362:         {1. / 2,   0,       0,      0,       0},
363:         {11. / 18, 1. / 18, 0,      0,       0},
364:         {5. / 6,   -5. / 6, .5,     0,       0},
365:         {1. / 4,   7. / 4,  3. / 4, -7. / 4, 0}
366:     },
367:                     At[5][5] = {{0, 0, 0, 0, 0}, {0, 1. / 2, 0, 0, 0}, {0, 1. / 6, 1. / 2, 0, 0}, {0, -1. / 2, 1. / 2, 1. / 2, 0}, {0, 3. / 2, -3. / 2, 1. / 2, 1. / 2}}, *bembedt = NULL;
368:     PetscCall(TSARKIMEXRegister(TSARKIMEXARS443, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL));
369:   }
370:   {
371:     const PetscReal A[5][5] =
372:       {
373:         {0,      0,      0,      0, 0},
374:         {1,      0,      0,      0, 0},
375:         {4. / 9, 2. / 9, 0,      0, 0},
376:         {1. / 4, 0,      3. / 4, 0, 0},
377:         {1. / 4, 0,      3. / 5, 0, 0}
378:     },
379:                     At[5][5] = {{0, 0, 0, 0, 0}, {.5, .5, 0, 0, 0}, {5. / 18, -1. / 9, .5, 0, 0}, {.5, 0, 0, .5, 0}, {.25, 0, .75, -.5, .5}}, *bembedt = NULL;
380:     PetscCall(TSARKIMEXRegister(TSARKIMEXBPR3, 3, 5, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 0, NULL, NULL));
381:   }
382:   {
383:     const PetscReal A[6][6] =
384:       {
385:         {0,                               0,                                 0,                                 0,                                0,              0},
386:         {1. / 2,                          0,                                 0,                                 0,                                0,              0},
387:         {13861. / 62500.,                 6889. / 62500.,                    0,                                 0,                                0,              0},
388:         {-116923316275. / 2393684061468., -2731218467317. / 15368042101831., 9408046702089. / 11113171139209.,  0,                                0,              0},
389:         {-451086348788. / 2902428689909., -2682348792572. / 7519795681897.,  12662868775082. / 11960479115383., 3355817975965. / 11060851509271., 0,              0},
390:         {647845179188. / 3216320057751.,  73281519250. / 8382639484533.,     552539513391. / 3454668386233.,    3354512671639. / 8306763924573.,  4040. / 17871., 0}
391:     },
392:                     At[6][6] = {{0, 0, 0, 0, 0, 0}, {1. / 4, 1. / 4, 0, 0, 0, 0}, {8611. / 62500., -1743. / 31250., 1. / 4, 0, 0, 0}, {5012029. / 34652500., -654441. / 2922500., 174375. / 388108., 1. / 4, 0, 0}, {15267082809. / 155376265600., -71443401. / 120774400., 730878875. / 902184768., 2285395. / 8070912., 1. / 4, 0}, {82889. / 524892., 0, 15625. / 83664., 69875. / 102672., -2260. / 8211, 1. / 4}}, bembedt[6] = {4586570599. / 29645900160., 0, 178811875. / 945068544., 814220225. / 1159782912., -3700637. / 11593932., 61727. / 225920.}, binterpt[6][3] = {{6943876665148. / 7220017795957., -54480133. / 30881146., 6818779379841. / 7100303317025.},  {0, 0, 0},
393:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     {7640104374378. / 9702883013639., -11436875. / 14766696., 2173542590792. / 12501825683035.}, {-20649996744609. / 7521556579894., 174696575. / 18121608., -31592104683404. / 5083833661969.},
394:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     {8854892464581. / 2390941311638., -12120380. / 966161., 61146701046299. / 7138195549469.},   {-11397109935349. / 6675773540249., 3843. / 706., -17219254887155. / 4939391667607.}};
395:     PetscCall(TSARKIMEXRegister(TSARKIMEX4, 4, 6, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL));
396:   }
397:   {
398:     const PetscReal A[8][8] =
399:       {
400:         {0,                                  0,                              0,                                 0,                                  0,                               0,                                 0,                               0},
401:         {41. / 100,                          0,                              0,                                 0,                                  0,                               0,                                 0,                               0},
402:         {367902744464. / 2072280473677.,     677623207551. / 8224143866563., 0,                                 0,                                  0,                               0,                                 0,                               0},
403:         {1268023523408. / 10340822734521.,   0,                              1029933939417. / 13636558850479.,  0,                                  0,                               0,                                 0,                               0},
404:         {14463281900351. / 6315353703477.,   0,                              66114435211212. / 5879490589093.,  -54053170152839. / 4284798021562.,  0,                               0,                                 0,                               0},
405:         {14090043504691. / 34967701212078.,  0,                              15191511035443. / 11219624916014., -18461159152457. / 12425892160975., -281667163811. / 9011619295870., 0,                                 0,                               0},
406:         {19230459214898. / 13134317526959.,  0,                              21275331358303. / 2942455364971.,  -38145345988419. / 4862620318723.,  -1. / 8,                         -1. / 8,                           0,                               0},
407:         {-19977161125411. / 11928030595625., 0,                              -40795976796054. / 6384907823539., 177454434618887. / 12078138498510., 782672205425. / 8267701900261.,  -69563011059811. / 9646580694205., 7356628210526. / 4942186776405., 0}
408:     },
409:                     At[8][8] = {{0, 0, 0, 0, 0, 0, 0, 0}, {41. / 200., 41. / 200., 0, 0, 0, 0, 0, 0}, {41. / 400., -567603406766. / 11931857230679., 41. / 200., 0, 0, 0, 0, 0}, {683785636431. / 9252920307686., 0, -110385047103. / 1367015193373., 41. / 200., 0, 0, 0, 0}, {3016520224154. / 10081342136671., 0, 30586259806659. / 12414158314087., -22760509404356. / 11113319521817., 41. / 200., 0, 0, 0}, {218866479029. / 1489978393911., 0, 638256894668. / 5436446318841., -1179710474555. / 5321154724896., -60928119172. / 8023461067671., 41. / 200., 0, 0}, {1020004230633. / 5715676835656., 0, 25762820946817. / 25263940353407., -2161375909145. / 9755907335909., -211217309593. / 5846859502534., -4269925059573. / 7827059040749., 41. / 200, 0}, {-872700587467. / 9133579230613., 0, 0, 22348218063261. / 9555858737531., -1143369518992. / 8141816002931., -39379526789629. / 19018526304540., 32727382324388. / 42900044865799., 41. / 200.}}, bembedt[8] = {-975461918565. / 9796059967033., 0, 0, 78070527104295. / 32432590147079., -548382580838. / 3424219808633., -33438840321285. / 15594753105479., 3629800801594. / 4656183773603., 4035322873751. / 18575991585200.}, binterpt[8][3] = {{-17674230611817. / 10670229744614., 43486358583215. / 12773830924787., -9257016797708. / 5021505065439.}, {0, 0, 0}, {0, 0, 0}, {65168852399939. / 7868540260826., -91478233927265. / 11067650958493., 26096422576131. / 11239449250142.}, {15494834004392. / 5936557850923., -79368583304911. / 10890268929626., 92396832856987. / 20362823103730.}, {-99329723586156. / 26959484932159., -12239297817655. / 9152339842473., 30029262896817. / 10175596800299.}, {-19024464361622. / 5461577185407., 115839755401235. / 10719374521269., -26136350496073. / 3983972220547.}, {-6511271360970. / 6095937251113., 5843115559534. / 2180450260947., -5289405421727. / 3760307252460.}};
410:     PetscCall(TSARKIMEXRegister(TSARKIMEX5, 5, 8, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembedt, 3, binterpt[0], NULL));
411:   }
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: /*@C
416:    TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by `TSARKIMEXRegister()`.

418:    Not Collective

420:    Level: advanced

422: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXRegister()`, `TSARKIMEXRegisterAll()`
423: @*/
424: PetscErrorCode TSARKIMEXRegisterDestroy(void)
425: {
426:   ARKTableauLink link;

428:   PetscFunctionBegin;
429:   while ((link = ARKTableauList)) {
430:     ARKTableau t   = &link->tab;
431:     ARKTableauList = link->next;
432:     PetscCall(PetscFree6(t->At, t->bt, t->ct, t->A, t->b, t->c));
433:     PetscCall(PetscFree2(t->bembedt, t->bembed));
434:     PetscCall(PetscFree2(t->binterpt, t->binterp));
435:     PetscCall(PetscFree(t->name));
436:     PetscCall(PetscFree(link));
437:   }
438:   TSARKIMEXRegisterAllCalled = PETSC_FALSE;
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: /*@C
443:   TSARKIMEXInitializePackage - This function initializes everything in the `TSARKIMEX` package. It is called
444:   from `TSInitializePackage()`.

446:   Level: developer

448: .seealso: [](ch_ts), `PetscInitialize()`, `TSARKIMEXFinalizePackage()`
449: @*/
450: PetscErrorCode TSARKIMEXInitializePackage(void)
451: {
452:   PetscFunctionBegin;
453:   if (TSARKIMEXPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
454:   TSARKIMEXPackageInitialized = PETSC_TRUE;
455:   PetscCall(TSARKIMEXRegisterAll());
456:   PetscCall(PetscRegisterFinalize(TSARKIMEXFinalizePackage));
457:   PetscFunctionReturn(PETSC_SUCCESS);
458: }

460: /*@C
461:   TSARKIMEXFinalizePackage - This function destroys everything in the `TSARKIMEX` package. It is
462:   called from `PetscFinalize()`.

464:   Level: developer

466: .seealso: [](ch_ts), `PetscFinalize()`, `TSARKIMEXInitializePackage()`
467: @*/
468: PetscErrorCode TSARKIMEXFinalizePackage(void)
469: {
470:   PetscFunctionBegin;
471:   TSARKIMEXPackageInitialized = PETSC_FALSE;
472:   PetscCall(TSARKIMEXRegisterDestroy());
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: /*@C
477:    TSARKIMEXRegister - register a `TSARKIMEX` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation

479:    Not Collective, but the same schemes should be registered on all processes on which they will be used

481:    Input Parameters:
482: +  name - identifier for method
483: .  order - approximation order of method
484: .  s - number of stages, this is the dimension of the matrices below
485: .  At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major)
486: .  bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At)
487: .  ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At)
488: .  A - Non-stiff stage coefficients (dimension s*s, row-major)
489: .  b - Non-stiff step completion table (dimension s; NULL to use last row of At)
490: .  c - Non-stiff abscissa (dimension s; NULL to use row sums of A)
491: .  bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available)
492: .  bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided)
493: .  pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp
494: .  binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp)
495: -  binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt)

497:    Level: advanced

499:    Note:
500:    Several `TSARKIMEX` methods are provided, this function is only needed to create new methods.

502: .seealso: [](ch_ts), `TSARKIMEX`, `TSType`, `TS`
503: @*/
504: PetscErrorCode TSARKIMEXRegister(TSARKIMEXType name, PetscInt order, PetscInt s, const PetscReal At[], const PetscReal bt[], const PetscReal ct[], const PetscReal A[], const PetscReal b[], const PetscReal c[], const PetscReal bembedt[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[], const PetscReal binterp[])
505: {
506:   ARKTableauLink link;
507:   ARKTableau     t;
508:   PetscInt       i, j;

510:   PetscFunctionBegin;
511:   PetscCall(TSARKIMEXInitializePackage());
512:   PetscCall(PetscNew(&link));
513:   t = &link->tab;
514:   PetscCall(PetscStrallocpy(name, &t->name));
515:   t->order = order;
516:   t->s     = s;
517:   PetscCall(PetscMalloc6(s * s, &t->At, s, &t->bt, s, &t->ct, s * s, &t->A, s, &t->b, s, &t->c));
518:   PetscCall(PetscArraycpy(t->At, At, s * s));
519:   PetscCall(PetscArraycpy(t->A, A, s * s));
520:   if (bt) PetscCall(PetscArraycpy(t->bt, bt, s));
521:   else
522:     for (i = 0; i < s; i++) t->bt[i] = At[(s - 1) * s + i];
523:   if (b) PetscCall(PetscArraycpy(t->b, b, s));
524:   else
525:     for (i = 0; i < s; i++) t->b[i] = t->bt[i];
526:   if (ct) PetscCall(PetscArraycpy(t->ct, ct, s));
527:   else
528:     for (i = 0; i < s; i++)
529:       for (j = 0, t->ct[i] = 0; j < s; j++) t->ct[i] += At[i * s + j];
530:   if (c) PetscCall(PetscArraycpy(t->c, c, s));
531:   else
532:     for (i = 0; i < s; i++)
533:       for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j];
534:   t->stiffly_accurate = PETSC_TRUE;
535:   for (i = 0; i < s; i++)
536:     if (t->At[(s - 1) * s + i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE;
537:   t->explicit_first_stage = PETSC_TRUE;
538:   for (i = 0; i < s; i++)
539:     if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE;
540:   /*def of FSAL can be made more precise*/
541:   t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate);
542:   if (bembedt) {
543:     PetscCall(PetscMalloc2(s, &t->bembedt, s, &t->bembed));
544:     PetscCall(PetscArraycpy(t->bembedt, bembedt, s));
545:     PetscCall(PetscArraycpy(t->bembed, bembed ? bembed : bembedt, s));
546:   }

548:   t->pinterp = pinterp;
549:   PetscCall(PetscMalloc2(s * pinterp, &t->binterpt, s * pinterp, &t->binterp));
550:   PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp));
551:   PetscCall(PetscArraycpy(t->binterp, binterp ? binterp : binterpt, s * pinterp));
552:   link->next     = ARKTableauList;
553:   ARKTableauList = link;
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*
558:  The step completion formula is

560:  x1 = x0 - h bt^T YdotI + h b^T YdotRHS

562:  This function can be called before or after ts->vec_sol has been updated.
563:  Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order.
564:  We can write

566:  x1e = x0 - h bet^T YdotI + h be^T YdotRHS
567:      = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS
568:      = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS

570:  so we can evaluate the method with different order even after the step has been optimistically completed.
571: */
572: static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts, PetscInt order, Vec X, PetscBool *done)
573: {
574:   TS_ARKIMEX  *ark = (TS_ARKIMEX *)ts->data;
575:   ARKTableau   tab = ark->tableau;
576:   PetscScalar *w   = ark->work;
577:   PetscReal    h;
578:   PetscInt     s = tab->s, j;

580:   PetscFunctionBegin;
581:   switch (ark->status) {
582:   case TS_STEP_INCOMPLETE:
583:   case TS_STEP_PENDING:
584:     h = ts->time_step;
585:     break;
586:   case TS_STEP_COMPLETE:
587:     h = ts->ptime - ts->ptime_prev;
588:     break;
589:   default:
590:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
591:   }
592:   if (order == tab->order) {
593:     if (ark->status == TS_STEP_INCOMPLETE) {
594:       if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */
595:         PetscCall(VecCopy(ark->Y[s - 1], X));
596:       } else { /* Use the standard completion formula (bt,b) */
597:         PetscCall(VecCopy(ts->vec_sol, X));
598:         for (j = 0; j < s; j++) w[j] = h * tab->bt[j];
599:         PetscCall(VecMAXPY(X, s, w, ark->YdotI));
600:         if (ark->imex) { /* Method is IMEX, complete the explicit formula */
601:           for (j = 0; j < s; j++) w[j] = h * tab->b[j];
602:           PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
603:         }
604:       }
605:     } else PetscCall(VecCopy(ts->vec_sol, X));
606:     if (done) *done = PETSC_TRUE;
607:     PetscFunctionReturn(PETSC_SUCCESS);
608:   } else if (order == tab->order - 1) {
609:     if (!tab->bembedt) goto unavailable;
610:     if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */
611:       PetscCall(VecCopy(ts->vec_sol, X));
612:       for (j = 0; j < s; j++) w[j] = h * tab->bembedt[j];
613:       PetscCall(VecMAXPY(X, s, w, ark->YdotI));
614:       for (j = 0; j < s; j++) w[j] = h * tab->bembed[j];
615:       PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
616:     } else { /* Rollback and re-complete using (bet-be,be-b) */
617:       PetscCall(VecCopy(ts->vec_sol, X));
618:       for (j = 0; j < s; j++) w[j] = h * (tab->bembedt[j] - tab->bt[j]);
619:       PetscCall(VecMAXPY(X, tab->s, w, ark->YdotI));
620:       for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]);
621:       PetscCall(VecMAXPY(X, s, w, ark->YdotRHS));
622:     }
623:     if (done) *done = PETSC_TRUE;
624:     PetscFunctionReturn(PETSC_SUCCESS);
625:   }
626: unavailable:
627:   if (done) *done = PETSC_FALSE;
628:   else
629:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "ARKIMEX '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", tab->name,
630:             tab->order, order);
631:   PetscFunctionReturn(PETSC_SUCCESS);
632: }

634: static PetscErrorCode TSARKIMEXTestMassIdentity(TS ts, PetscBool *id)
635: {
636:   Vec         Udot, Y1, Y2;
637:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
638:   PetscReal   norm;

640:   PetscFunctionBegin;
641:   PetscCall(VecDuplicate(ts->vec_sol, &Udot));
642:   PetscCall(VecDuplicate(ts->vec_sol, &Y1));
643:   PetscCall(VecDuplicate(ts->vec_sol, &Y2));
644:   PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y1, ark->imex));
645:   PetscCall(VecSetRandom(Udot, NULL));
646:   PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, Udot, Y2, ark->imex));
647:   PetscCall(VecAXPY(Y2, -1.0, Y1));
648:   PetscCall(VecAXPY(Y2, -1.0, Udot));
649:   PetscCall(VecNorm(Y2, NORM_2, &norm));
650:   if (norm < 100.0 * PETSC_MACHINE_EPSILON) {
651:     *id = PETSC_TRUE;
652:   } else {
653:     *id = PETSC_FALSE;
654:     PetscCall(PetscInfo((PetscObject)ts, "IFunction(Udot = random) - IFunction(Udot = 0) is not near Udot, %g, suspect mass matrix implied in IFunction() is not the identity as required\n", (double)norm));
655:   }
656:   PetscCall(VecDestroy(&Udot));
657:   PetscCall(VecDestroy(&Y1));
658:   PetscCall(VecDestroy(&Y2));
659:   PetscFunctionReturn(PETSC_SUCCESS);
660: }

662: static PetscErrorCode TSRollBack_ARKIMEX(TS ts)
663: {
664:   TS_ARKIMEX      *ark = (TS_ARKIMEX *)ts->data;
665:   ARKTableau       tab = ark->tableau;
666:   const PetscInt   s   = tab->s;
667:   const PetscReal *bt = tab->bt, *b = tab->b;
668:   PetscScalar     *w     = ark->work;
669:   Vec             *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS;
670:   PetscInt         j;
671:   PetscReal        h;

673:   PetscFunctionBegin;
674:   switch (ark->status) {
675:   case TS_STEP_INCOMPLETE:
676:   case TS_STEP_PENDING:
677:     h = ts->time_step;
678:     break;
679:   case TS_STEP_COMPLETE:
680:     h = ts->ptime - ts->ptime_prev;
681:     break;
682:   default:
683:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
684:   }
685:   for (j = 0; j < s; j++) w[j] = -h * bt[j];
686:   PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotI));
687:   for (j = 0; j < s; j++) w[j] = -h * b[j];
688:   PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS));
689:   PetscFunctionReturn(PETSC_SUCCESS);
690: }

692: static PetscErrorCode TSStep_ARKIMEX(TS ts)
693: {
694:   TS_ARKIMEX      *ark = (TS_ARKIMEX *)ts->data;
695:   ARKTableau       tab = ark->tableau;
696:   const PetscInt   s   = tab->s;
697:   const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c;
698:   PetscScalar     *w = ark->work;
699:   Vec             *Y = ark->Y, *YdotI = ark->YdotI, *YdotRHS = ark->YdotRHS, Ydot = ark->Ydot, Ydot0 = ark->Ydot0, Z = ark->Z;
700:   PetscBool        extrapolate = ark->extrapolate;
701:   TSAdapt          adapt;
702:   SNES             snes;
703:   PetscInt         i, j, its, lits;
704:   PetscInt         rejections = 0;
705:   PetscBool        stageok, accept = PETSC_TRUE;
706:   PetscReal        next_time_step = ts->time_step;

708:   PetscFunctionBegin;
709:   if (ark->extrapolate && !ark->Y_prev) {
710:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
711:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
712:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
713:   }

715:   if (!ts->steprollback) {
716:     if (ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */
717:       PetscCall(VecCopy(YdotI[s - 1], Ydot0));
718:     }
719:     if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */
720:       for (i = 0; i < s; i++) {
721:         PetscCall(VecCopy(Y[i], ark->Y_prev[i]));
722:         PetscCall(VecCopy(YdotRHS[i], ark->YdotRHS_prev[i]));
723:         PetscCall(VecCopy(YdotI[i], ark->YdotI_prev[i]));
724:       }
725:     }
726:   }

728:   if (ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) {
729:     TS ts_start;
730:     if (PetscDefined(USE_DEBUG)) {
731:       PetscBool id = PETSC_FALSE;
732:       PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
733:       PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "This scheme requires an identity mass matrix, however the TSIFunction you provide does not utilize an identity mass matrix");
734:     }
735:     PetscCall(TSClone(ts, &ts_start));
736:     PetscCall(TSSetSolution(ts_start, ts->vec_sol));
737:     PetscCall(TSSetTime(ts_start, ts->ptime));
738:     PetscCall(TSSetMaxSteps(ts_start, ts->steps + 1));
739:     PetscCall(TSSetMaxTime(ts_start, ts->ptime + ts->time_step));
740:     PetscCall(TSSetExactFinalTime(ts_start, TS_EXACTFINALTIME_STEPOVER));
741:     PetscCall(TSSetTimeStep(ts_start, ts->time_step));
742:     PetscCall(TSSetType(ts_start, TSARKIMEX));
743:     PetscCall(TSARKIMEXSetFullyImplicit(ts_start, PETSC_TRUE));
744:     PetscCall(TSARKIMEXSetType(ts_start, TSARKIMEX1BEE));

746:     PetscCall(TSRestartStep(ts_start));
747:     PetscCall(TSSolve(ts_start, ts->vec_sol));
748:     PetscCall(TSGetTime(ts_start, &ts->ptime));
749:     PetscCall(TSGetTimeStep(ts_start, &ts->time_step));

751:     { /* Save the initial slope for the next step */
752:       TS_ARKIMEX *ark_start = (TS_ARKIMEX *)ts_start->data;
753:       PetscCall(VecCopy(ark_start->YdotI[ark_start->tableau->s - 1], Ydot0));
754:     }
755:     ts->steps++;

757:     /* Set the correct TS in SNES */
758:     /* We'll try to bypass this by changing the method on the fly */
759:     {
760:       PetscCall(TSGetSNES(ts, &snes));
761:       PetscCall(TSSetSNES(ts, snes));
762:     }
763:     PetscCall(TSDestroy(&ts_start));
764:   }

766:   ark->status = TS_STEP_INCOMPLETE;
767:   while (!ts->reason && ark->status != TS_STEP_COMPLETE) {
768:     PetscReal t = ts->ptime;
769:     PetscReal h = ts->time_step;
770:     for (i = 0; i < s; i++) {
771:       ark->stage_time = t + h * ct[i];
772:       PetscCall(TSPreStage(ts, ark->stage_time));
773:       if (At[i * s + i] == 0) { /* This stage is explicit */
774:         PetscCheck(i == 0 || ts->equation_type < TS_EQ_IMPLICIT, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Explicit stages other than the first one are not supported for implicit problems");
775:         PetscCall(VecCopy(ts->vec_sol, Y[i]));
776:         for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
777:         PetscCall(VecMAXPY(Y[i], i, w, YdotI));
778:         for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
779:         PetscCall(VecMAXPY(Y[i], i, w, YdotRHS));
780:       } else {
781:         ark->scoeff = 1. / At[i * s + i];
782:         /* Ydot = shift*(Y-Z) */
783:         PetscCall(VecCopy(ts->vec_sol, Z));
784:         for (j = 0; j < i; j++) w[j] = h * At[i * s + j];
785:         PetscCall(VecMAXPY(Z, i, w, YdotI));
786:         for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
787:         PetscCall(VecMAXPY(Z, i, w, YdotRHS));
788:         if (extrapolate && !ts->steprestart) {
789:           /* Initial guess extrapolated from previous time step stage values */
790:           PetscCall(TSExtrapolate_ARKIMEX(ts, c[i], Y[i]));
791:         } else {
792:           /* Initial guess taken from last stage */
793:           PetscCall(VecCopy(i > 0 ? Y[i - 1] : ts->vec_sol, Y[i]));
794:         }
795:         PetscCall(TSGetSNES(ts, &snes));
796:         PetscCall(SNESSolve(snes, NULL, Y[i]));
797:         PetscCall(SNESGetIterationNumber(snes, &its));
798:         PetscCall(SNESGetLinearSolveIterations(snes, &lits));
799:         ts->snes_its += its;
800:         ts->ksp_its += lits;
801:         PetscCall(TSGetAdapt(ts, &adapt));
802:         PetscCall(TSAdaptCheckStage(adapt, ts, ark->stage_time, Y[i], &stageok));
803:         if (!stageok) {
804:           /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to
805:            * use extrapolation to initialize the solves on the next attempt. */
806:           extrapolate = PETSC_FALSE;
807:           goto reject_step;
808:         }
809:       }
810:       if (ts->equation_type >= TS_EQ_IMPLICIT) {
811:         if (i == 0 && tab->explicit_first_stage) {
812:           PetscCheck(tab->stiffly_accurate, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s is not stiffly accurate and therefore explicit-first stage methods cannot be used if the equation is implicit because the slope cannot be evaluated",
813:                      ark->tableau->name);
814:           PetscCall(VecCopy(Ydot0, YdotI[0])); /* YdotI = YdotI(tn-1) */
815:         } else {
816:           PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
817:         }
818:       } else {
819:         if (i == 0 && tab->explicit_first_stage) {
820:           PetscCall(VecZeroEntries(Ydot));
821:           PetscCall(TSComputeIFunction(ts, t + h * ct[i], Y[i], Ydot, YdotI[i], ark->imex)); /* YdotI = -G(t,Y,0)   */
822:           PetscCall(VecScale(YdotI[i], -1.0));
823:         } else {
824:           PetscCall(VecAXPBYPCZ(YdotI[i], -ark->scoeff / h, ark->scoeff / h, 0, Z, Y[i])); /* YdotI = shift*(X-Z) */
825:         }
826:         if (ark->imex) {
827:           PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i]));
828:         } else {
829:           PetscCall(VecZeroEntries(YdotRHS[i]));
830:         }
831:       }
832:       PetscCall(TSPostStage(ts, ark->stage_time, i, Y));
833:     }

835:     ark->status = TS_STEP_INCOMPLETE;
836:     PetscCall(TSEvaluateStep_ARKIMEX(ts, tab->order, ts->vec_sol, NULL));
837:     ark->status = TS_STEP_PENDING;
838:     PetscCall(TSGetAdapt(ts, &adapt));
839:     PetscCall(TSAdaptCandidatesClear(adapt));
840:     PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
841:     PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
842:     ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
843:     if (!accept) { /* Roll back the current step */
844:       PetscCall(TSRollBack_ARKIMEX(ts));
845:       ts->time_step = next_time_step;
846:       goto reject_step;
847:     }

849:     ts->ptime += ts->time_step;
850:     ts->time_step = next_time_step;
851:     break;

853:   reject_step:
854:     ts->reject++;
855:     accept = PETSC_FALSE;
856:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
857:       ts->reason = TS_DIVERGED_STEP_REJECTED;
858:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
859:     }
860:   }
861:   PetscFunctionReturn(PETSC_SUCCESS);
862: }
863: /*
864:   This adjoint step function assumes the partitioned ODE system has an identity mass matrix and thus can be represented in the form
865:   Udot = H(t,U) + G(t,U)
866:   This corresponds to F(t,U,Udot) = Udot-H(t,U).

868:   The complete adjoint equations are
869:   (shift*I - dHdu) lambda_s[i]   = 1/at[i][i] (
870:     (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j]
871:     + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j]),  i = s-1,...,0
872:   lambda_n = lambda_{n+1} + \sum_{j=1}^s lambda_s[j]
873:   mu_{n+1}[i]   = h (at[i][i] dHdP lambda_s[i]
874:     + (b_i dGdP + bt[i] dHdP) lambda_{n+1} + dGdP \sum_{j=i+1}^s a[j][i] lambda_s[j]
875:     + dHdP \sum_{j=i+1}^s at[j][i] lambda_s[j]), i = s-1,...,0
876:   mu_n = mu_{n+1} + \sum_{j=1}^s mu_{n+1}[j]
877:   where shift = 1/(at[i][i]*h)

879:   If at[i][i] is 0, the first equation falls back to
880:   lambda_s[i] = h (
881:     (b_i dGdU + bt[i] dHdU) lambda_{n+1} + dGdU \sum_{j=i+1}^s a[j][i] lambda_s[j]
882:     + dHdU \sum_{j=i+1}^s at[j][i] lambda_s[j])

884: */
885: static PetscErrorCode TSAdjointStep_ARKIMEX(TS ts)
886: {
887:   TS_ARKIMEX      *ark = (TS_ARKIMEX *)ts->data;
888:   ARKTableau       tab = ark->tableau;
889:   const PetscInt   s   = tab->s;
890:   const PetscReal *At = tab->At, *A = tab->A, *ct = tab->ct, *c = tab->c, *b = tab->b, *bt = tab->bt;
891:   PetscScalar     *w = ark->work;
892:   Vec             *Y = ark->Y, Ydot = ark->Ydot, *VecsDeltaLam = ark->VecsDeltaLam, *VecsSensiTemp = ark->VecsSensiTemp, *VecsSensiPTemp = ark->VecsSensiPTemp;
893:   Mat              Jex, Jim, Jimpre;
894:   PetscInt         i, j, nadj;
895:   PetscReal        t                 = ts->ptime, stage_time_ex;
896:   PetscReal        adjoint_time_step = -ts->time_step; /* always positive since ts->time_step is negative */
897:   KSP              ksp;

899:   PetscFunctionBegin;
900:   ark->status = TS_STEP_INCOMPLETE;
901:   PetscCall(SNESGetKSP(ts->snes, &ksp));
902:   PetscCall(TSGetRHSMats_Private(ts, &Jex, NULL));
903:   PetscCall(TSGetIJacobian(ts, &Jim, &Jimpre, NULL, NULL));

905:   for (i = s - 1; i >= 0; i--) {
906:     ark->stage_time = t - adjoint_time_step * (1.0 - ct[i]);
907:     stage_time_ex   = t - adjoint_time_step * (1.0 - c[i]);
908:     if (At[i * s + i] == 0) { // This stage is explicit
909:       ark->scoeff = 0.;
910:     } else {
911:       ark->scoeff = -1. / At[i * s + i]; // this makes shift=ark->scoeff/ts->time_step positive since ts->time_step is negative
912:     }
913:     PetscCall(TSComputeSNESJacobian(ts, Y[i], Jim, Jimpre));
914:     PetscCall(TSComputeRHSJacobian(ts, stage_time_ex, Y[i], Jex, Jex));
915:     if (ts->vecs_sensip) {
916:       PetscCall(TSComputeIJacobianP(ts, ark->stage_time, Y[i], Ydot, ark->scoeff / adjoint_time_step, ts->Jacp, PETSC_TRUE)); // get dFdP (-dHdP), Ydot not really used since mass matrix is identity
917:       PetscCall(TSComputeRHSJacobianP(ts, stage_time_ex, Y[i], ts->Jacprhs));                                                 // get dGdP
918:     }
919:     for (nadj = 0; nadj < ts->numcost; nadj++) {
920:       /* Build RHS (stored in VecsDeltaLam) for first-order adjoint */
921:       if (s - i - 1 > 0) {
922:         /* Temp = -\sum_{j=i+1}^s at[j][i] lambda_s[j] */
923:         for (j = i + 1; j < s; j++) w[j - i - 1] = -At[j * s + i];
924:         PetscCall(VecSet(VecsSensiTemp[nadj], 0));
925:         PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
926:         /* (shift I - dHdU) Temp */
927:         PetscCall(MatMultTranspose(Jim, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i]));
928:         /* cancel out shift Temp where shift=-scoeff/h */
929:         PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], ark->scoeff / adjoint_time_step, VecsSensiTemp[nadj]));
930:         if (ts->vecs_sensip) {
931:           /* - dHdP Temp */
932:           PetscCall(MatMultTranspose(ts->Jacp, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
933:           /* mu_n += h dHdP Temp */
934:           PetscCall(VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsSensiPTemp[nadj]));
935:         }
936:         /* Temp = \sum_{j=i+1}^s a[j][i] lambda_s[j] */
937:         for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i];
938:         PetscCall(VecSet(VecsSensiTemp[nadj], 0));
939:         PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
940:         /* dGdU Temp */
941:         PetscCall(MatMultTransposeAdd(Jex, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
942:         if (ts->vecs_sensip) {
943:           /* dGdP Temp */
944:           PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecsSensiPTemp[nadj]));
945:           /* mu_n += h dGdP Temp */
946:           PetscCall(VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, VecsSensiPTemp[nadj]));
947:         }
948:       } else {
949:         PetscCall(VecSet(VecsDeltaLam[nadj * s + i], 0)); // make sure it is initialized
950:       }
951:       if (bt[i]) { // bt[i] dHdU lambda_{n+1}
952:         /* (shift I - dHdU)^T lambda_{n+1} */
953:         PetscCall(MatMultTranspose(Jim, ts->vecs_sensi[nadj], VecsSensiTemp[nadj]));
954:         /* -bt[i] (shift I - dHdU)^T lambda_{n+1} */
955:         PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i], VecsSensiTemp[nadj]));
956:         /* cancel out -bt[i] shift lambda_{n+1} */
957:         PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -bt[i] * ark->scoeff / adjoint_time_step, ts->vecs_sensi[nadj]));
958:         if (ts->vecs_sensip) {
959:           /* -dHdP lambda_{n+1} */
960:           PetscCall(MatMultTranspose(ts->Jacp, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj]));
961:           /* mu_n += h bt[i] dHdP lambda_{n+1} */
962:           PetscCall(VecAXPY(ts->vecs_sensip[nadj], -bt[i] * adjoint_time_step, VecsSensiPTemp[nadj]));
963:         }
964:       }
965:       if (b[i]) { // b[i] dGdU lambda_{n+1}
966:         /* dGdU lambda_{n+1} */
967:         PetscCall(MatMultTranspose(Jex, ts->vecs_sensi[nadj], VecsSensiTemp[nadj]));
968:         /* b[i] dGdU lambda_{n+1} */
969:         PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], b[i], VecsSensiTemp[nadj]));
970:         if (ts->vecs_sensip) {
971:           /* dGdP lambda_{n+1} */
972:           PetscCall(MatMultTranspose(ts->Jacprhs, ts->vecs_sensi[nadj], VecsSensiPTemp[nadj]));
973:           /* mu_n += h b[i] dGdP lambda_{n+1} */
974:           PetscCall(VecAXPY(ts->vecs_sensip[nadj], b[i] * adjoint_time_step, VecsSensiPTemp[nadj]));
975:         }
976:       }
977:       /* Build LHS for first-order adjoint */
978:       if (At[i * s + i] == 0) { // This stage is explicit
979:         PetscCall(VecScale(VecsDeltaLam[nadj * s + i], adjoint_time_step));
980:       } else {
981:         KSP                ksp;
982:         KSPConvergedReason kspreason;
983:         PetscCall(SNESGetKSP(ts->snes, &ksp));
984:         PetscCall(KSPSetOperators(ksp, Jim, Jimpre));
985:         PetscCall(VecScale(VecsDeltaLam[nadj * s + i], 1. / At[i * s + i]));
986:         PetscCall(KSPSolveTranspose(ksp, VecsDeltaLam[nadj * s + i], VecsDeltaLam[nadj * s + i]));
987:         PetscCall(KSPGetConvergedReason(ksp, &kspreason));
988:         if (kspreason < 0) {
989:           ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
990:           PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj));
991:         }
992:         if (ts->vecs_sensip) {
993:           /* -dHdP lambda_s[i] */
994:           PetscCall(MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj * s + i], VecsSensiPTemp[nadj]));
995:           /* mu_n += h at[i][i] dHdP lambda_s[i] */
996:           PetscCall(VecAXPY(ts->vecs_sensip[nadj], -At[i * s + i] * adjoint_time_step, VecsSensiPTemp[nadj]));
997:         }
998:       }
999:     }
1000:   }
1001:   for (j = 0; j < s; j++) w[j] = 1.0;
1002:   for (nadj = 0; nadj < ts->numcost; nadj++) // no need to do this for mu's
1003:     PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s]));
1004:   ark->status = TS_STEP_COMPLETE;
1005:   PetscFunctionReturn(PETSC_SUCCESS);
1006: }

1008: static PetscErrorCode TSInterpolate_ARKIMEX(TS ts, PetscReal itime, Vec X)
1009: {
1010:   TS_ARKIMEX      *ark = (TS_ARKIMEX *)ts->data;
1011:   PetscInt         s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j;
1012:   PetscReal        h;
1013:   PetscReal        tt, t;
1014:   PetscScalar     *bt, *b;
1015:   const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp;

1017:   PetscFunctionBegin;
1018:   PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name);
1019:   switch (ark->status) {
1020:   case TS_STEP_INCOMPLETE:
1021:   case TS_STEP_PENDING:
1022:     h = ts->time_step;
1023:     t = (itime - ts->ptime) / h;
1024:     break;
1025:   case TS_STEP_COMPLETE:
1026:     h = ts->ptime - ts->ptime_prev;
1027:     t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1028:     break;
1029:   default:
1030:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1031:   }
1032:   PetscCall(PetscMalloc2(s, &bt, s, &b));
1033:   for (i = 0; i < s; i++) bt[i] = b[i] = 0;
1034:   for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1035:     for (i = 0; i < s; i++) {
1036:       bt[i] += h * Bt[i * pinterp + j] * tt;
1037:       b[i] += h * B[i * pinterp + j] * tt;
1038:     }
1039:   }
1040:   PetscCall(VecCopy(ark->Y[0], X));
1041:   PetscCall(VecMAXPY(X, s, bt, ark->YdotI));
1042:   PetscCall(VecMAXPY(X, s, b, ark->YdotRHS));
1043:   PetscCall(PetscFree2(bt, b));
1044:   PetscFunctionReturn(PETSC_SUCCESS);
1045: }

1047: static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts, PetscReal c, Vec X)
1048: {
1049:   TS_ARKIMEX      *ark = (TS_ARKIMEX *)ts->data;
1050:   PetscInt         s = ark->tableau->s, pinterp = ark->tableau->pinterp, i, j;
1051:   PetscReal        h, h_prev, t, tt;
1052:   PetscScalar     *bt, *b;
1053:   const PetscReal *Bt = ark->tableau->binterpt, *B = ark->tableau->binterp;

1055:   PetscFunctionBegin;
1056:   PetscCheck(Bt && B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSARKIMEX %s does not have an interpolation formula", ark->tableau->name);
1057:   PetscCall(PetscCalloc2(s, &bt, s, &b));
1058:   h      = ts->time_step;
1059:   h_prev = ts->ptime - ts->ptime_prev;
1060:   t      = 1 + h / h_prev * c;
1061:   for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1062:     for (i = 0; i < s; i++) {
1063:       bt[i] += h * Bt[i * pinterp + j] * tt;
1064:       b[i] += h * B[i * pinterp + j] * tt;
1065:     }
1066:   }
1067:   PetscCheck(ark->Y_prev, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Stages from previous step have not been stored");
1068:   PetscCall(VecCopy(ark->Y_prev[0], X));
1069:   PetscCall(VecMAXPY(X, s, bt, ark->YdotI_prev));
1070:   PetscCall(VecMAXPY(X, s, b, ark->YdotRHS_prev));
1071:   PetscCall(PetscFree2(bt, b));
1072:   PetscFunctionReturn(PETSC_SUCCESS);
1073: }

1075: /*------------------------------------------------------------*/

1077: static PetscErrorCode TSARKIMEXTableauReset(TS ts)
1078: {
1079:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1080:   ARKTableau  tab = ark->tableau;

1082:   PetscFunctionBegin;
1083:   if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1084:   PetscCall(PetscFree(ark->work));
1085:   PetscCall(VecDestroyVecs(tab->s, &ark->Y));
1086:   PetscCall(VecDestroyVecs(tab->s, &ark->YdotI));
1087:   PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS));
1088:   PetscCall(VecDestroyVecs(tab->s, &ark->Y_prev));
1089:   PetscCall(VecDestroyVecs(tab->s, &ark->YdotI_prev));
1090:   PetscCall(VecDestroyVecs(tab->s, &ark->YdotRHS_prev));
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

1094: static PetscErrorCode TSReset_ARKIMEX(TS ts)
1095: {
1096:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1098:   PetscFunctionBegin;
1099:   PetscCall(TSARKIMEXTableauReset(ts));
1100:   PetscCall(VecDestroy(&ark->Ydot));
1101:   PetscCall(VecDestroy(&ark->Ydot0));
1102:   PetscCall(VecDestroy(&ark->Z));
1103:   PetscFunctionReturn(PETSC_SUCCESS);
1104: }

1106: static PetscErrorCode TSAdjointReset_ARKIMEX(TS ts)
1107: {
1108:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1109:   ARKTableau  tab = ark->tableau;

1111:   PetscFunctionBegin;
1112:   PetscCall(VecDestroyVecs(tab->s * ts->numcost, &ark->VecsDeltaLam));
1113:   PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiTemp));
1114:   PetscCall(VecDestroyVecs(ts->numcost, &ark->VecsSensiPTemp));
1115:   PetscFunctionReturn(PETSC_SUCCESS);
1116: }

1118: static PetscErrorCode TSARKIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1119: {
1120:   TS_ARKIMEX *ax = (TS_ARKIMEX *)ts->data;

1122:   PetscFunctionBegin;
1123:   if (Z) {
1124:     if (dm && dm != ts->dm) {
1125:       PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1126:     } else *Z = ax->Z;
1127:   }
1128:   if (Ydot) {
1129:     if (dm && dm != ts->dm) {
1130:       PetscCall(DMGetNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1131:     } else *Ydot = ax->Ydot;
1132:   }
1133:   PetscFunctionReturn(PETSC_SUCCESS);
1134: }

1136: static PetscErrorCode TSARKIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot)
1137: {
1138:   PetscFunctionBegin;
1139:   if (Z) {
1140:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Z", Z));
1141:   }
1142:   if (Ydot) {
1143:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSARKIMEX_Ydot", Ydot));
1144:   }
1145:   PetscFunctionReturn(PETSC_SUCCESS);
1146: }

1148: /*
1149:   This defines the nonlinear equation that is to be solved with SNES
1150:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
1151: */
1152: static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes, Vec X, Vec F, TS ts)
1153: {
1154:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1155:   DM          dm, dmsave;
1156:   Vec         Z, Ydot;
1157:   PetscReal   shift = ark->scoeff / ts->time_step;

1159:   PetscFunctionBegin;
1160:   PetscCall(SNESGetDM(snes, &dm));
1161:   PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, &Ydot));
1162:   PetscCall(VecAXPBYPCZ(Ydot, -shift, shift, 0, Z, X)); /* Ydot = shift*(X-Z) */
1163:   dmsave = ts->dm;
1164:   ts->dm = dm;

1166:   PetscCall(TSComputeIFunction(ts, ark->stage_time, X, Ydot, F, ark->imex));

1168:   ts->dm = dmsave;
1169:   PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, &Ydot));
1170:   PetscFunctionReturn(PETSC_SUCCESS);
1171: }

1173: static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts)
1174: {
1175:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1176:   DM          dm, dmsave;
1177:   Vec         Ydot;
1178:   PetscReal   shift = ark->scoeff / ts->time_step;

1180:   PetscFunctionBegin;
1181:   PetscCall(SNESGetDM(snes, &dm));
1182:   PetscCall(TSARKIMEXGetVecs(ts, dm, NULL, &Ydot));
1183:   /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */
1184:   dmsave = ts->dm;
1185:   ts->dm = dm;

1187:   PetscCall(TSComputeIJacobian(ts, ark->stage_time, X, Ydot, shift, A, B, ark->imex));

1189:   ts->dm = dmsave;
1190:   PetscCall(TSARKIMEXRestoreVecs(ts, dm, NULL, &Ydot));
1191:   PetscFunctionReturn(PETSC_SUCCESS);
1192: }

1194: static PetscErrorCode TSGetStages_ARKIMEX(TS ts, PetscInt *ns, Vec *Y[])
1195: {
1196:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1198:   PetscFunctionBegin;
1199:   if (ns) *ns = ark->tableau->s;
1200:   if (Y) *Y = ark->Y;
1201:   PetscFunctionReturn(PETSC_SUCCESS);
1202: }

1204: static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine, DM coarse, void *ctx)
1205: {
1206:   PetscFunctionBegin;
1207:   PetscFunctionReturn(PETSC_SUCCESS);
1208: }

1210: static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1211: {
1212:   TS  ts = (TS)ctx;
1213:   Vec Z, Z_c;

1215:   PetscFunctionBegin;
1216:   PetscCall(TSARKIMEXGetVecs(ts, fine, &Z, NULL));
1217:   PetscCall(TSARKIMEXGetVecs(ts, coarse, &Z_c, NULL));
1218:   PetscCall(MatRestrict(restrct, Z, Z_c));
1219:   PetscCall(VecPointwiseMult(Z_c, rscale, Z_c));
1220:   PetscCall(TSARKIMEXRestoreVecs(ts, fine, &Z, NULL));
1221:   PetscCall(TSARKIMEXRestoreVecs(ts, coarse, &Z_c, NULL));
1222:   PetscFunctionReturn(PETSC_SUCCESS);
1223: }

1225: static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm, DM subdm, void *ctx)
1226: {
1227:   PetscFunctionBegin;
1228:   PetscFunctionReturn(PETSC_SUCCESS);
1229: }

1231: static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1232: {
1233:   TS  ts = (TS)ctx;
1234:   Vec Z, Z_c;

1236:   PetscFunctionBegin;
1237:   PetscCall(TSARKIMEXGetVecs(ts, dm, &Z, NULL));
1238:   PetscCall(TSARKIMEXGetVecs(ts, subdm, &Z_c, NULL));

1240:   PetscCall(VecScatterBegin(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));
1241:   PetscCall(VecScatterEnd(gscat, Z, Z_c, INSERT_VALUES, SCATTER_FORWARD));

1243:   PetscCall(TSARKIMEXRestoreVecs(ts, dm, &Z, NULL));
1244:   PetscCall(TSARKIMEXRestoreVecs(ts, subdm, &Z_c, NULL));
1245:   PetscFunctionReturn(PETSC_SUCCESS);
1246: }

1248: static PetscErrorCode TSARKIMEXTableauSetUp(TS ts)
1249: {
1250:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1251:   ARKTableau  tab = ark->tableau;

1253:   PetscFunctionBegin;
1254:   PetscCall(PetscMalloc1(tab->s, &ark->work));
1255:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y));
1256:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI));
1257:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS));
1258:   if (ark->extrapolate) {
1259:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->Y_prev));
1260:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotI_prev));
1261:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ark->YdotRHS_prev));
1262:   }
1263:   PetscFunctionReturn(PETSC_SUCCESS);
1264: }

1266: static PetscErrorCode TSSetUp_ARKIMEX(TS ts)
1267: {
1268:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1269:   DM          dm;
1270:   SNES        snes;

1272:   PetscFunctionBegin;
1273:   PetscCall(TSARKIMEXTableauSetUp(ts));
1274:   PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot));
1275:   PetscCall(VecDuplicate(ts->vec_sol, &ark->Ydot0));
1276:   PetscCall(VecDuplicate(ts->vec_sol, &ark->Z));
1277:   PetscCall(TSGetDM(ts, &dm));
1278:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
1279:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
1280:   PetscCall(TSGetSNES(ts, &snes));
1281:   PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1283: /*------------------------------------------------------------*/

1285: static PetscErrorCode TSAdjointSetUp_ARKIMEX(TS ts)
1286: {
1287:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1288:   ARKTableau  tab = ark->tableau;

1290:   PetscFunctionBegin;
1291:   PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], tab->s * ts->numcost, &ark->VecsDeltaLam));
1292:   PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &ark->VecsSensiTemp));
1293:   if (ts->vecs_sensip) { PetscCall(VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &ark->VecsSensiPTemp)); }
1294:   if (PetscDefined(USE_DEBUG)) {
1295:     PetscBool id = PETSC_FALSE;
1296:     PetscCall(TSARKIMEXTestMassIdentity(ts, &id));
1297:     PetscCheck(id, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Adjoint ARKIMEX requires an identity mass matrix, however the TSIFunction you provide does not utilize an identity mass matrix");
1298:   }
1299:   PetscFunctionReturn(PETSC_SUCCESS);
1300: }

1302: static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts, PetscOptionItems *PetscOptionsObject)
1303: {
1304:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1306:   PetscFunctionBegin;
1307:   PetscOptionsHeadBegin(PetscOptionsObject, "ARKIMEX ODE solver options");
1308:   {
1309:     ARKTableauLink link;
1310:     PetscInt       count, choice;
1311:     PetscBool      flg;
1312:     const char   **namelist;
1313:     for (link = ARKTableauList, count = 0; link; link = link->next, count++)
1314:       ;
1315:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
1316:     for (link = ARKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1317:     PetscCall(PetscOptionsEList("-ts_arkimex_type", "Family of ARK IMEX method", "TSARKIMEXSetType", (const char *const *)namelist, count, ark->tableau->name, &choice, &flg));
1318:     if (flg) PetscCall(TSARKIMEXSetType(ts, namelist[choice]));
1319:     PetscCall(PetscFree(namelist));

1321:     flg = (PetscBool)!ark->imex;
1322:     PetscCall(PetscOptionsBool("-ts_arkimex_fully_implicit", "Solve the problem fully implicitly", "TSARKIMEXSetFullyImplicit", flg, &flg, NULL));
1323:     ark->imex = (PetscBool)!flg;
1324:     PetscCall(PetscOptionsBool("-ts_arkimex_initial_guess_extrapolate", "Extrapolate the initial guess for the stage solution from stage values of the previous time step", "", ark->extrapolate, &ark->extrapolate, NULL));
1325:   }
1326:   PetscOptionsHeadEnd();
1327:   PetscFunctionReturn(PETSC_SUCCESS);
1328: }

1330: static PetscErrorCode TSView_ARKIMEX(TS ts, PetscViewer viewer)
1331: {
1332:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;
1333:   PetscBool   iascii;

1335:   PetscFunctionBegin;
1336:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1337:   if (iascii) {
1338:     ARKTableau    tab = ark->tableau;
1339:     TSARKIMEXType arktype;
1340:     char          buf[512];
1341:     PetscBool     flg;

1343:     PetscCall(TSARKIMEXGetType(ts, &arktype));
1344:     PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg));
1345:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ARK IMEX %s\n", arktype));
1346:     PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ct));
1347:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Stiff abscissa       ct = %s\n", buf));
1348:     PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c));
1349:     PetscCall(PetscViewerASCIIPrintf(viewer, "Fully implicit: %s\n", flg ? "yes" : "no"));
1350:     PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s\n", tab->stiffly_accurate ? "yes" : "no"));
1351:     PetscCall(PetscViewerASCIIPrintf(viewer, "Explicit first stage: %s\n", tab->explicit_first_stage ? "yes" : "no"));
1352:     PetscCall(PetscViewerASCIIPrintf(viewer, "FSAL property: %s\n", tab->FSAL_implicit ? "yes" : "no"));
1353:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Nonstiff abscissa     c = %s\n", buf));
1354:   }
1355:   PetscFunctionReturn(PETSC_SUCCESS);
1356: }

1358: static PetscErrorCode TSLoad_ARKIMEX(TS ts, PetscViewer viewer)
1359: {
1360:   SNES    snes;
1361:   TSAdapt adapt;

1363:   PetscFunctionBegin;
1364:   PetscCall(TSGetAdapt(ts, &adapt));
1365:   PetscCall(TSAdaptLoad(adapt, viewer));
1366:   PetscCall(TSGetSNES(ts, &snes));
1367:   PetscCall(SNESLoad(snes, viewer));
1368:   /* function and Jacobian context for SNES when used with TS is always ts object */
1369:   PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
1370:   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
1371:   PetscFunctionReturn(PETSC_SUCCESS);
1372: }

1374: /*@C
1375:   TSARKIMEXSetType - Set the type of `TSARKIMEX` scheme

1377:   Logically Collective

1379:   Input Parameters:
1380: +  ts - timestepping context
1381: -  arktype - type of `TSARKIMEX` scheme

1383:   Options Database Key:
1384: .  -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5> - set `TSARKIMEX` scheme type

1386:   Level: intermediate

1388: .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEX`, `TSARKIMEXType`, `TSARKIMEX1BEE`, `TSARKIMEXA2`, `TSARKIMEXL2`, `TSARKIMEXARS122`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEXPRSSP2`,
1389:           `TSARKIMEX3`, `TSARKIMEXBPR3`, `TSARKIMEXARS443`, `TSARKIMEX4`, `TSARKIMEX5`
1390: @*/
1391: PetscErrorCode TSARKIMEXSetType(TS ts, TSARKIMEXType arktype)
1392: {
1393:   PetscFunctionBegin;
1396:   PetscTryMethod(ts, "TSARKIMEXSetType_C", (TS, TSARKIMEXType), (ts, arktype));
1397:   PetscFunctionReturn(PETSC_SUCCESS);
1398: }

1400: /*@C
1401:   TSARKIMEXGetType - Get the type of `TSARKIMEX` scheme

1403:   Logically Collective

1405:   Input Parameter:
1406: .  ts - timestepping context

1408:   Output Parameter:
1409: .  arktype - type of `TSARKIMEX` scheme

1411:   Level: intermediate

1413: .seealso: [](ch_ts), `TSARKIMEX`c, `TSARKIMEXGetType()`
1414: @*/
1415: PetscErrorCode TSARKIMEXGetType(TS ts, TSARKIMEXType *arktype)
1416: {
1417:   PetscFunctionBegin;
1419:   PetscUseMethod(ts, "TSARKIMEXGetType_C", (TS, TSARKIMEXType *), (ts, arktype));
1420:   PetscFunctionReturn(PETSC_SUCCESS);
1421: }

1423: /*@
1424:   TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly, including the part that is normally solved explicitly

1426:   Logically Collective

1428:   Input Parameters:
1429: +  ts - timestepping context
1430: -  flg - `PETSC_TRUE` for fully implicit

1432:   Level: intermediate

1434: .seealso: [](ch_ts), `TSARKIMEX`, `TSARKIMEXGetType()`, `TSARKIMEXGetFullyImplicit()`
1435: @*/
1436: PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts, PetscBool flg)
1437: {
1438:   PetscFunctionBegin;
1441:   PetscTryMethod(ts, "TSARKIMEXSetFullyImplicit_C", (TS, PetscBool), (ts, flg));
1442:   PetscFunctionReturn(PETSC_SUCCESS);
1443: }

1445: /*@
1446:   TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly

1448:   Logically Collective

1450:   Input Parameter:
1451: .  ts - timestepping context

1453:   Output Parameter:
1454: .  flg - `PETSC_TRUE` for fully implicit

1456:   Level: intermediate

1458: .seealso: [](ch_ts), `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`
1459: @*/
1460: PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts, PetscBool *flg)
1461: {
1462:   PetscFunctionBegin;
1465:   PetscUseMethod(ts, "TSARKIMEXGetFullyImplicit_C", (TS, PetscBool *), (ts, flg));
1466:   PetscFunctionReturn(PETSC_SUCCESS);
1467: }

1469: static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts, TSARKIMEXType *arktype)
1470: {
1471:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1473:   PetscFunctionBegin;
1474:   *arktype = ark->tableau->name;
1475:   PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1477: static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts, TSARKIMEXType arktype)
1478: {
1479:   TS_ARKIMEX    *ark = (TS_ARKIMEX *)ts->data;
1480:   PetscBool      match;
1481:   ARKTableauLink link;

1483:   PetscFunctionBegin;
1484:   if (ark->tableau) {
1485:     PetscCall(PetscStrcmp(ark->tableau->name, arktype, &match));
1486:     if (match) PetscFunctionReturn(PETSC_SUCCESS);
1487:   }
1488:   for (link = ARKTableauList; link; link = link->next) {
1489:     PetscCall(PetscStrcmp(link->tab.name, arktype, &match));
1490:     if (match) {
1491:       if (ts->setupcalled) PetscCall(TSARKIMEXTableauReset(ts));
1492:       ark->tableau = &link->tab;
1493:       if (ts->setupcalled) PetscCall(TSARKIMEXTableauSetUp(ts));
1494:       ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1495:       PetscFunctionReturn(PETSC_SUCCESS);
1496:     }
1497:   }
1498:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", arktype);
1499: }

1501: static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts, PetscBool flg)
1502: {
1503:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1505:   PetscFunctionBegin;
1506:   ark->imex = (PetscBool)!flg;
1507:   PetscFunctionReturn(PETSC_SUCCESS);
1508: }

1510: static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts, PetscBool *flg)
1511: {
1512:   TS_ARKIMEX *ark = (TS_ARKIMEX *)ts->data;

1514:   PetscFunctionBegin;
1515:   *flg = (PetscBool)!ark->imex;
1516:   PetscFunctionReturn(PETSC_SUCCESS);
1517: }

1519: static PetscErrorCode TSDestroy_ARKIMEX(TS ts)
1520: {
1521:   PetscFunctionBegin;
1522:   PetscCall(TSReset_ARKIMEX(ts));
1523:   if (ts->dm) {
1524:     PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSARKIMEX, DMRestrictHook_TSARKIMEX, ts));
1525:     PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSARKIMEX, DMSubDomainRestrictHook_TSARKIMEX, ts));
1526:   }
1527:   PetscCall(PetscFree(ts->data));
1528:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", NULL));
1529:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", NULL));
1530:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", NULL));
1531:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", NULL));
1532:   PetscFunctionReturn(PETSC_SUCCESS);
1533: }

1535: /* ------------------------------------------------------------ */
1536: /*MC
1537:       TSARKIMEX - ODE and DAE solver using additive Runge-Kutta IMEX schemes

1539:   These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1540:   nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1541:   of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.

1543:   Level: beginner

1545:   Notes:
1546:   The default is `TSARKIMEX3`, it can be changed with `TSARKIMEXSetType()` or -ts_arkimex_type

1548:   If the equation is implicit or a DAE, then `TSSetEquationType()` needs to be set accordingly. Refer to the manual for further information.

1550:   Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).

1552:   Consider trying `TSROSW` if the stiff part is linear or weakly nonlinear.

1554: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSARKIMEXSetType()`, `TSARKIMEXGetType()`, `TSARKIMEXSetFullyImplicit()`, `TSARKIMEXGetFullyImplicit()`,
1555:           `TSARKIMEX1BEE`, `TSARKIMEX2C`, `TSARKIMEX2D`, `TSARKIMEX2E`, `TSARKIMEX3`, `TSARKIMEXL2`, `TSARKIMEXA2`, `TSARKIMEXARS122`,
1556:           `TSARKIMEX4`, `TSARKIMEX5`, `TSARKIMEXPRSSP2`, `TSARKIMEXARS443`, `TSARKIMEXBPR3`, `TSARKIMEXType`, `TSARKIMEXRegister()`, `TSType`
1557: M*/
1558: PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts)
1559: {
1560:   TS_ARKIMEX *ark;

1562:   PetscFunctionBegin;
1563:   PetscCall(TSARKIMEXInitializePackage());

1565:   ts->ops->reset          = TSReset_ARKIMEX;
1566:   ts->ops->adjointreset   = TSAdjointReset_ARKIMEX;
1567:   ts->ops->destroy        = TSDestroy_ARKIMEX;
1568:   ts->ops->view           = TSView_ARKIMEX;
1569:   ts->ops->load           = TSLoad_ARKIMEX;
1570:   ts->ops->setup          = TSSetUp_ARKIMEX;
1571:   ts->ops->adjointsetup   = TSAdjointSetUp_ARKIMEX;
1572:   ts->ops->step           = TSStep_ARKIMEX;
1573:   ts->ops->interpolate    = TSInterpolate_ARKIMEX;
1574:   ts->ops->evaluatestep   = TSEvaluateStep_ARKIMEX;
1575:   ts->ops->rollback       = TSRollBack_ARKIMEX;
1576:   ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX;
1577:   ts->ops->snesfunction   = SNESTSFormFunction_ARKIMEX;
1578:   ts->ops->snesjacobian   = SNESTSFormJacobian_ARKIMEX;
1579:   ts->ops->getstages      = TSGetStages_ARKIMEX;
1580:   ts->ops->adjointstep    = TSAdjointStep_ARKIMEX;

1582:   ts->usessnes = PETSC_TRUE;

1584:   PetscCall(PetscNew(&ark));
1585:   ts->data  = (void *)ark;
1586:   ark->imex = PETSC_TRUE;

1588:   ark->VecsDeltaLam   = NULL;
1589:   ark->VecsSensiTemp  = NULL;
1590:   ark->VecsSensiPTemp = NULL;

1592:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetType_C", TSARKIMEXGetType_ARKIMEX));
1593:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetType_C", TSARKIMEXSetType_ARKIMEX));
1594:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXSetFullyImplicit_C", TSARKIMEXSetFullyImplicit_ARKIMEX));
1595:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSARKIMEXGetFullyImplicit_C", TSARKIMEXGetFullyImplicit_ARKIMEX));

1597:   PetscCall(TSARKIMEXSetType(ts, TSARKIMEXDefault));
1598:   PetscFunctionReturn(PETSC_SUCCESS);
1599: }