Actual source code: rosw.c

  1: /*
  2:   Code for timestepping with Rosenbrock W methods

  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.
 10:   This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian.

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

 16: #include <petsc/private/kernels/blockinvert.h>

 18: static TSRosWType TSRosWDefault = TSROSWRA34PW2;
 19: static PetscBool  TSRosWRegisterAllCalled;
 20: static PetscBool  TSRosWPackageInitialized;

 22: typedef struct _RosWTableau *RosWTableau;
 23: struct _RosWTableau {
 24:   char      *name;
 25:   PetscInt   order;             /* Classical approximation order of the method */
 26:   PetscInt   s;                 /* Number of stages */
 27:   PetscInt   pinterp;           /* Interpolation order */
 28:   PetscReal *A;                 /* Propagation table, strictly lower triangular */
 29:   PetscReal *Gamma;             /* Stage table, lower triangular with nonzero diagonal */
 30:   PetscBool *GammaZeroDiag;     /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */
 31:   PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/
 32:   PetscReal *b;                 /* Step completion table */
 33:   PetscReal *bembed;            /* Step completion table for embedded method of order one less */
 34:   PetscReal *ASum;              /* Row sum of A */
 35:   PetscReal *GammaSum;          /* Row sum of Gamma, only needed for non-autonomous systems */
 36:   PetscReal *At;                /* Propagation table in transformed variables */
 37:   PetscReal *bt;                /* Step completion table in transformed variables */
 38:   PetscReal *bembedt;           /* Step completion table of order one less in transformed variables */
 39:   PetscReal *GammaInv;          /* Inverse of Gamma, used for transformed variables */
 40:   PetscReal  ccfl;              /* Placeholder for CFL coefficient relative to forward Euler */
 41:   PetscReal *binterpt;          /* Dense output formula */
 42: };
 43: typedef struct _RosWTableauLink *RosWTableauLink;
 44: struct _RosWTableauLink {
 45:   struct _RosWTableau tab;
 46:   RosWTableauLink     next;
 47: };
 48: static RosWTableauLink RosWTableauList;

 50: typedef struct {
 51:   RosWTableau  tableau;
 52:   Vec         *Y;            /* States computed during the step, used to complete the step */
 53:   Vec          Ydot;         /* Work vector holding Ydot during residual evaluation */
 54:   Vec          Ystage;       /* Work vector for the state value at each stage */
 55:   Vec          Zdot;         /* Ydot = Zdot + shift*Y */
 56:   Vec          Zstage;       /* Y = Zstage + Y */
 57:   Vec          vec_sol_prev; /* Solution from the previous step (used for interpolation and rollback)*/
 58:   PetscScalar *work;         /* Scalar work space of length number of stages, used to prepare VecMAXPY() */
 59:   PetscReal    scoeff;       /* shift = scoeff/dt */
 60:   PetscReal    stage_time;
 61:   PetscReal    stage_explicit;     /* Flag indicates that the current stage is explicit */
 62:   PetscBool    recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
 63:   TSStepStatus status;
 64: } TS_RosW;

 66: /*MC
 67:      TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method).

 69:      Only an approximate Jacobian is needed.

 71:      Level: intermediate

 73: .seealso: [](ch_ts), `TSROSW`
 74: M*/

 76: /*MC
 77:      TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method).

 79:      Only an approximate Jacobian is needed.

 81:      Level: intermediate

 83: .seealso: [](ch_ts), `TSROSW`
 84: M*/

 86: /*MC
 87:      TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme.

 89:      Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P.

 91:      Level: intermediate

 93: .seealso: [](ch_ts), `TSROSW`
 94: M*/

 96: /*MC
 97:      TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme.

 99:      Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M.

101:      Level: intermediate

103: .seealso: [](ch_ts), `TSROSW`
104: M*/

106: /*MC
107:      TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`

109:      Only an approximate Jacobian is needed. By default, it is only recomputed once per step.

111:      This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73.

113:      Level: intermediate

115: .seealso: [](ch_ts), `TSROSW`
116: M*/

118: /*MC
119:      TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`.

121:      Only an approximate Jacobian is needed. By default, it is only recomputed once per step.

123:      This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48.

125:      Level: intermediate

127: .seealso: [](ch_ts), `TSROSW`
128: M*/

130: /*MC
131:      TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`

133:      By default, the Jacobian is only recomputed once per step.

135:      Both the third order and embedded second order methods are stiffly accurate and L-stable.

137:      Level: intermediate

139: .seealso: [](ch_ts), `TSROSW`, `TSROSWSANDU3`
140: M*/

142: /*MC
143:      TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`

145:      By default, the Jacobian is only recomputed once per step.

147:      The third order method is L-stable, but not stiffly accurate.
148:      The second order embedded method is strongly A-stable with R(infty) = 0.5.
149:      The internal stages are L-stable.
150:      This method is called ROS3 in {cite}`sandu_1997`.

152:      Level: intermediate

154: .seealso: [](ch_ts), `TSROSW`, `TSROSWRODAS3`
155: M*/

157: /*MC
158:      TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages

160:      By default, the Jacobian is only recomputed once per step.

162:      A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3)

164:      Level: intermediate

166: .seealso: [](ch_ts), `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP`
167: M*/

169: /*MC
170:      TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages

172:      By default, the Jacobian is only recomputed once per step.

174:      L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)

176:      Level: intermediate

178: .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP`
179: M*/

181: /*MC
182:      TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages

184:      By default, the Jacobian is only recomputed once per step.

186:      L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)

188:      Level: intermediate

190: .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP`
191: M*/

193: /*MC
194:      TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop {cite}`kaps1979generalized`

196:      By default, the Jacobian is only recomputed once per step.

198:      A(89.3 degrees)-stable, |R(infty)| = 0.454.

200:      This method does not provide a dense output formula.

202:      Level: intermediate

204:      Note:
205:      See Section 4 Table 7.2 in {cite}`wanner1996solving`

207: .seealso: [](ch_ts), `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`
208: M*/

210: /*MC
211:      TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine {cite}`shampine1982implementation`

213:      By default, the Jacobian is only recomputed once per step.

215:      A-stable, |R(infty)| = 1/3.

217:      This method does not provide a dense output formula.

219:      Level: intermediate

221:      Note:
222:      See Section 4 Table 7.2 in {cite}`wanner1996solving`

224: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L`
225: M*/

227: /*MC
228:      TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen {cite}`veldhuizen1984d`

230:      By default, the Jacobian is only recomputed once per step.

232:      A(89.5 degrees)-stable, |R(infty)| = 0.24.

234:      This method does not provide a dense output formula.

236:      Level: intermediate

238:      Note:
239:      See Section 4 Table 7.2 in {cite}`wanner1996solving`

241: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
242: M*/

244: /*MC
245:      TSROSW4L - four stage, fourth order Rosenbrock (not W) method

247:      By default, the Jacobian is only recomputed once per step.

249:      A-stable and L-stable

251:      This method does not provide a dense output formula.

253:      Level: intermediate

255:      Note:
256:      See Section 4 Table 7.2 in {cite}`wanner1996solving`

258: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
259: M*/

261: /*@C
262:   TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in `TSROSW`

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

266:   Level: advanced

268: .seealso: [](ch_ts), `TSROSW`, `TSRosWRegisterDestroy()`
269: @*/
270: PetscErrorCode TSRosWRegisterAll(void)
271: {
272:   PetscFunctionBegin;
273:   if (TSRosWRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
274:   TSRosWRegisterAllCalled = PETSC_TRUE;

276:   {
277:     const PetscReal A        = 0;
278:     const PetscReal Gamma    = 1;
279:     const PetscReal b        = 1;
280:     const PetscReal binterpt = 1;

282:     PetscCall(TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
283:   }

285:   {
286:     const PetscReal A        = 0;
287:     const PetscReal Gamma    = 0.5;
288:     const PetscReal b        = 1;
289:     const PetscReal binterpt = 1;

291:     PetscCall(TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
292:   }

294:   {
295:     /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0);   Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */
296:     const PetscReal A[2][2] = {
297:       {0,  0},
298:       {1., 0}
299:     };
300:     const PetscReal Gamma[2][2] = {
301:       {1.707106781186547524401,       0                      },
302:       {-2. * 1.707106781186547524401, 1.707106781186547524401}
303:     };
304:     const PetscReal b[2]  = {0.5, 0.5};
305:     const PetscReal b1[2] = {1.0, 0.0};
306:     PetscReal       binterpt[2][2];
307:     binterpt[0][0] = 1.707106781186547524401 - 1.0;
308:     binterpt[1][0] = 2.0 - 1.707106781186547524401;
309:     binterpt[0][1] = 1.707106781186547524401 - 1.5;
310:     binterpt[1][1] = 1.5 - 1.707106781186547524401;

312:     PetscCall(TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]));
313:   }
314:   {
315:     /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0);   Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */
316:     const PetscReal A[2][2] = {
317:       {0,  0},
318:       {1., 0}
319:     };
320:     const PetscReal Gamma[2][2] = {
321:       {0.2928932188134524755992,       0                       },
322:       {-2. * 0.2928932188134524755992, 0.2928932188134524755992}
323:     };
324:     const PetscReal b[2]  = {0.5, 0.5};
325:     const PetscReal b1[2] = {1.0, 0.0};
326:     PetscReal       binterpt[2][2];
327:     binterpt[0][0] = 0.2928932188134524755992 - 1.0;
328:     binterpt[1][0] = 2.0 - 0.2928932188134524755992;
329:     binterpt[0][1] = 0.2928932188134524755992 - 1.5;
330:     binterpt[1][1] = 1.5 - 0.2928932188134524755992;

332:     PetscCall(TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]));
333:   }
334:   {
335:     /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */
336:     PetscReal       binterpt[3][2];
337:     const PetscReal A[3][3] = {
338:       {0,                      0, 0},
339:       {1.5773502691896257e+00, 0, 0},
340:       {0.5,                    0, 0}
341:     };
342:     const PetscReal Gamma[3][3] = {
343:       {7.8867513459481287e-01,  0,                       0                     },
344:       {-1.5773502691896257e+00, 7.8867513459481287e-01,  0                     },
345:       {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01}
346:     };
347:     const PetscReal b[3]  = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01};
348:     const PetscReal b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01};

350:     binterpt[0][0] = -0.8094010767585034;
351:     binterpt[1][0] = -0.5;
352:     binterpt[2][0] = 2.3094010767585034;
353:     binterpt[0][1] = 0.9641016151377548;
354:     binterpt[1][1] = 0.5;
355:     binterpt[2][1] = -1.4641016151377548;

357:     PetscCall(TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
358:   }
359:   {
360:     PetscReal binterpt[4][3];
361:     /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
362:     const PetscReal A[4][4] = {
363:       {0,                      0,                       0,  0},
364:       {8.7173304301691801e-01, 0,                       0,  0},
365:       {8.4457060015369423e-01, -1.1299064236484185e-01, 0,  0},
366:       {0,                      0,                       1., 0}
367:     };
368:     const PetscReal Gamma[4][4] = {
369:       {4.3586652150845900e-01,  0,                       0,                      0                     },
370:       {-8.7173304301691801e-01, 4.3586652150845900e-01,  0,                      0                     },
371:       {-9.0338057013044082e-01, 5.4180672388095326e-02,  4.3586652150845900e-01, 0                     },
372:       {2.4212380706095346e-01,  -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01}
373:     };
374:     const PetscReal b[4]  = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01};
375:     const PetscReal b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01};

377:     binterpt[0][0] = 1.0564298455794094;
378:     binterpt[1][0] = 2.296429974281067;
379:     binterpt[2][0] = -1.307599564525376;
380:     binterpt[3][0] = -1.045260255335102;
381:     binterpt[0][1] = -1.3864882699759573;
382:     binterpt[1][1] = -8.262611700275677;
383:     binterpt[2][1] = 7.250979895056055;
384:     binterpt[3][1] = 2.398120075195581;
385:     binterpt[0][2] = 0.5721822314575016;
386:     binterpt[1][2] = 4.742931142090097;
387:     binterpt[2][2] = -4.398120075195578;
388:     binterpt[3][2] = -0.9169932983520199;

390:     PetscCall(TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
391:   }
392:   {
393:     /* const PetscReal g = 0.5;       Directly written in-place below */
394:     const PetscReal A[4][4] = {
395:       {0,    0,     0,   0},
396:       {0,    0,     0,   0},
397:       {1.,   0,     0,   0},
398:       {0.75, -0.25, 0.5, 0}
399:     };
400:     const PetscReal Gamma[4][4] = {
401:       {0.5,     0,       0,       0  },
402:       {1.,      0.5,     0,       0  },
403:       {-0.25,   -0.25,   0.5,     0  },
404:       {1. / 12, 1. / 12, -2. / 3, 0.5}
405:     };
406:     const PetscReal b[4]  = {5. / 6, -1. / 6, -1. / 6, 0.5};
407:     const PetscReal b2[4] = {0.75, -0.25, 0.5, 0};

409:     PetscCall(TSRosWRegister(TSROSWRODAS3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
410:   }
411:   {
412:     /*const PetscReal g = 0.43586652150845899941601945119356;       Directly written in-place below */
413:     const PetscReal A[3][3] = {
414:       {0,                                  0, 0},
415:       {0.43586652150845899941601945119356, 0, 0},
416:       {0.43586652150845899941601945119356, 0, 0}
417:     };
418:     const PetscReal Gamma[3][3] = {
419:       {0.43586652150845899941601945119356,  0,                                  0                                 },
420:       {-0.19294655696029095575009695436041, 0.43586652150845899941601945119356, 0                                 },
421:       {0,                                   1.74927148125794685173529749738960, 0.43586652150845899941601945119356}
422:     };
423:     const PetscReal b[3]  = {-0.75457412385404315829818998646589, 1.94100407061964420292840123379419, -0.18642994676560104463021124732829};
424:     const PetscReal b2[3] = {-1.53358745784149585370766523913002, 2.81745131148625772213931745457622, -0.28386385364476186843165221544619};

426:     PetscReal binterpt[3][2];
427:     binterpt[0][0] = 3.793692883777660870425141387941;
428:     binterpt[1][0] = -2.918692883777660870425141387941;
429:     binterpt[2][0] = 0.125;
430:     binterpt[0][1] = -0.725741064379812106687651020584;
431:     binterpt[1][1] = 0.559074397713145440020984353917;
432:     binterpt[2][1] = 0.16666666666666666666666666666667;

434:     PetscCall(TSRosWRegister(TSROSWSANDU3, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
435:   }
436:   {
437:     /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0;
438:      * Direct evaluation: s3 = 1.732050807568877293527;
439:      *                     g = 0.7886751345948128822546;
440:      * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */
441:     const PetscReal A[3][3] = {
442:       {0,    0,    0},
443:       {1,    0,    0},
444:       {0.25, 0.25, 0}
445:     };
446:     const PetscReal Gamma[3][3] = {
447:       {0,                                       0,                                      0                       },
448:       {(-3.0 - 1.732050807568877293527) / 6.0,  0.7886751345948128822546,               0                       },
449:       {(-3.0 - 1.732050807568877293527) / 24.0, (-3.0 - 1.732050807568877293527) / 8.0, 0.7886751345948128822546}
450:     };
451:     const PetscReal b[3]  = {1. / 6., 1. / 6., 2. / 3.};
452:     const PetscReal b2[3] = {1. / 4., 1. / 4., 1. / 2.};
453:     PetscReal       binterpt[3][2];

455:     binterpt[0][0] = 0.089316397477040902157517886164709;
456:     binterpt[1][0] = -0.91068360252295909784248211383529;
457:     binterpt[2][0] = 1.8213672050459181956849642276706;
458:     binterpt[0][1] = 0.077350269189625764509148780501957;
459:     binterpt[1][1] = 1.077350269189625764509148780502;
460:     binterpt[2][1] = -1.1547005383792515290182975610039;

462:     PetscCall(TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
463:   }

465:   {
466:     const PetscReal A[4][4] = {
467:       {0,       0,       0,       0},
468:       {1. / 2., 0,       0,       0},
469:       {1. / 2., 1. / 2., 0,       0},
470:       {1. / 6., 1. / 6., 1. / 6., 0}
471:     };
472:     const PetscReal Gamma[4][4] = {
473:       {1. / 2., 0,        0,       0},
474:       {0.0,     1. / 4.,  0,       0},
475:       {-2.,     -2. / 3., 2. / 3., 0},
476:       {1. / 2., 5. / 36., -2. / 9, 0}
477:     };
478:     const PetscReal b[4]  = {1. / 6., 1. / 6., 1. / 6., 1. / 2.};
479:     const PetscReal b2[4] = {1. / 8., 3. / 4., 1. / 8., 0};
480:     PetscReal       binterpt[4][3];

482:     binterpt[0][0] = 6.25;
483:     binterpt[1][0] = -30.25;
484:     binterpt[2][0] = 1.75;
485:     binterpt[3][0] = 23.25;
486:     binterpt[0][1] = -9.75;
487:     binterpt[1][1] = 58.75;
488:     binterpt[2][1] = -3.25;
489:     binterpt[3][1] = -45.75;
490:     binterpt[0][2] = 3.6666666666666666666666666666667;
491:     binterpt[1][2] = -28.333333333333333333333333333333;
492:     binterpt[2][2] = 1.6666666666666666666666666666667;
493:     binterpt[3][2] = 23.;

495:     PetscCall(TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
496:   }

498:   {
499:     const PetscReal A[4][4] = {
500:       {0,       0,       0,       0},
501:       {1. / 2., 0,       0,       0},
502:       {1. / 2., 1. / 2., 0,       0},
503:       {1. / 6., 1. / 6., 1. / 6., 0}
504:     };
505:     const PetscReal Gamma[4][4] = {
506:       {1. / 2.,  0,          0,        0},
507:       {0.0,      3. / 4.,    0,        0},
508:       {-2. / 3., -23. / 9.,  2. / 9.,  0},
509:       {1. / 18., 65. / 108., -2. / 27, 0}
510:     };
511:     const PetscReal b[4]  = {1. / 6., 1. / 6., 1. / 6., 1. / 2.};
512:     const PetscReal b2[4] = {3. / 16., 10. / 16., 3. / 16., 0};
513:     PetscReal       binterpt[4][3];

515:     binterpt[0][0] = 1.6911764705882352941176470588235;
516:     binterpt[1][0] = 3.6813725490196078431372549019608;
517:     binterpt[2][0] = 0.23039215686274509803921568627451;
518:     binterpt[3][0] = -4.6029411764705882352941176470588;
519:     binterpt[0][1] = -0.95588235294117647058823529411765;
520:     binterpt[1][1] = -6.2401960784313725490196078431373;
521:     binterpt[2][1] = -0.31862745098039215686274509803922;
522:     binterpt[3][1] = 7.5147058823529411764705882352941;
523:     binterpt[0][2] = -0.56862745098039215686274509803922;
524:     binterpt[1][2] = 2.7254901960784313725490196078431;
525:     binterpt[2][2] = 0.25490196078431372549019607843137;
526:     binterpt[3][2] = -2.4117647058823529411764705882353;

528:     PetscCall(TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
529:   }

531:   {
532:     PetscReal A[4][4], Gamma[4][4], b[4], b2[4];
533:     PetscReal binterpt[4][3];

535:     Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816;
536:     Gamma[0][1] = 0;
537:     Gamma[0][2] = 0;
538:     Gamma[0][3] = 0;
539:     Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476;
540:     Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816;
541:     Gamma[1][2] = 0;
542:     Gamma[1][3] = 0;
543:     Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903;
544:     Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131;
545:     Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816;
546:     Gamma[2][3] = 0;
547:     Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783;
548:     Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984;
549:     Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198;
550:     Gamma[3][3] = 0;

552:     A[0][0] = 0;
553:     A[0][1] = 0;
554:     A[0][2] = 0;
555:     A[0][3] = 0;
556:     A[1][0] = 0.8717330430169179988320388950590125027645343373957631;
557:     A[1][1] = 0;
558:     A[1][2] = 0;
559:     A[1][3] = 0;
560:     A[2][0] = 0.5275890119763004115618079766722914408876108660811028;
561:     A[2][1] = 0.07241098802369958843819203208518599088698057726988732;
562:     A[2][2] = 0;
563:     A[2][3] = 0;
564:     A[3][0] = 0.3990960076760701320627260685975778145384666450351314;
565:     A[3][1] = -0.4375576546135194437228463747348862825846903771419953;
566:     A[3][2] = 1.038461646937449311660120300601880176655352737312713;
567:     A[3][3] = 0;

569:     b[0] = 0.1876410243467238251612921333138006734899663569186926;
570:     b[1] = -0.5952974735769549480478230473706443582188442040780541;
571:     b[2] = 0.9717899277217721234705114616271378792182450260943198;
572:     b[3] = 0.4358665215084589994160194475295062513822671686978816;

574:     b2[0] = 0.2147402862233891404862383521089097657790734483804460;
575:     b2[1] = -0.4851622638849390928209050538171743017757490232519684;
576:     b2[2] = 0.8687250025203875511662123688667549217531982787600080;
577:     b2[3] = 0.4016969751411624011684543450940068201770721128357014;

579:     binterpt[0][0] = 2.2565812720167954547104627844105;
580:     binterpt[1][0] = 1.349166413351089573796243820819;
581:     binterpt[2][0] = -2.4695174540533503758652847586647;
582:     binterpt[3][0] = -0.13623023131453465264142184656474;
583:     binterpt[0][1] = -3.0826699111559187902922463354557;
584:     binterpt[1][1] = -2.4689115685996042534544925650515;
585:     binterpt[2][1] = 5.7428279814696677152129332773553;
586:     binterpt[3][1] = -0.19124650171414467146619437684812;
587:     binterpt[0][2] = 1.0137296634858471607430756831148;
588:     binterpt[1][2] = 0.52444768167155973161042570784064;
589:     binterpt[2][2] = -2.3015205996945452158771370439586;
590:     binterpt[3][2] = 0.76334325453713832352363565300308;

592:     PetscCall(TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
593:   }
594:   PetscCall(TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_DEFAULT, PETSC_DEFAULT, 0, -0.1282612945269037e+01));
595:   PetscCall(TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_DEFAULT, PETSC_DEFAULT, 0, 125. / 108.));
596:   PetscCall(TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.355958941201148));
597:   PetscCall(TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.093502252409163));
598:   PetscFunctionReturn(PETSC_SUCCESS);
599: }

601: /*@C
602:   TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`.

604:   Not Collective

606:   Level: advanced

608: .seealso: [](ch_ts), `TSRosWRegister()`, `TSRosWRegisterAll()`
609: @*/
610: PetscErrorCode TSRosWRegisterDestroy(void)
611: {
612:   RosWTableauLink link;

614:   PetscFunctionBegin;
615:   while ((link = RosWTableauList)) {
616:     RosWTableau t   = &link->tab;
617:     RosWTableauList = link->next;
618:     PetscCall(PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum));
619:     PetscCall(PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr));
620:     PetscCall(PetscFree2(t->bembed, t->bembedt));
621:     PetscCall(PetscFree(t->binterpt));
622:     PetscCall(PetscFree(t->name));
623:     PetscCall(PetscFree(link));
624:   }
625:   TSRosWRegisterAllCalled = PETSC_FALSE;
626:   PetscFunctionReturn(PETSC_SUCCESS);
627: }

629: /*@C
630:   TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called
631:   from `TSInitializePackage()`.

633:   Level: developer

635: .seealso: [](ch_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()`
636: @*/
637: PetscErrorCode TSRosWInitializePackage(void)
638: {
639:   PetscFunctionBegin;
640:   if (TSRosWPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
641:   TSRosWPackageInitialized = PETSC_TRUE;
642:   PetscCall(TSRosWRegisterAll());
643:   PetscCall(PetscRegisterFinalize(TSRosWFinalizePackage));
644:   PetscFunctionReturn(PETSC_SUCCESS);
645: }

647: /*@C
648:   TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is
649:   called from `PetscFinalize()`.

651:   Level: developer

653: .seealso: [](ch_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()`
654: @*/
655: PetscErrorCode TSRosWFinalizePackage(void)
656: {
657:   PetscFunctionBegin;
658:   TSRosWPackageInitialized = PETSC_FALSE;
659:   PetscCall(TSRosWRegisterDestroy());
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: /*@C
664:   TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation

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

668:   Input Parameters:
669: + name     - identifier for method
670: . order    - approximation order of method
671: . s        - number of stages, this is the dimension of the matrices below
672: . A        - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular
673: . Gamma    - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal
674: . b        - Step completion table (dimension s)
675: . bembed   - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available)
676: . pinterp  - Order of the interpolation scheme, equal to the number of columns of binterpt
677: - binterpt - Coefficients of the interpolation formula (dimension s*pinterp)

679:   Level: advanced

681:   Note:
682:   Several Rosenbrock W methods are provided, this function is only needed to create new methods.

684: .seealso: [](ch_ts), `TSROSW`
685: @*/
686: PetscErrorCode TSRosWRegister(TSRosWType name, PetscInt order, PetscInt s, const PetscReal A[], const PetscReal Gamma[], const PetscReal b[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[])
687: {
688:   RosWTableauLink link;
689:   RosWTableau     t;
690:   PetscInt        i, j, k;
691:   PetscScalar    *GammaInv;

693:   PetscFunctionBegin;
694:   PetscAssertPointer(name, 1);
695:   PetscAssertPointer(A, 4);
696:   PetscAssertPointer(Gamma, 5);
697:   PetscAssertPointer(b, 6);
698:   if (bembed) PetscAssertPointer(bembed, 7);

700:   PetscCall(TSRosWInitializePackage());
701:   PetscCall(PetscNew(&link));
702:   t = &link->tab;
703:   PetscCall(PetscStrallocpy(name, &t->name));
704:   t->order = order;
705:   t->s     = s;
706:   PetscCall(PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum));
707:   PetscCall(PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr));
708:   PetscCall(PetscArraycpy(t->A, A, s * s));
709:   PetscCall(PetscArraycpy(t->Gamma, Gamma, s * s));
710:   PetscCall(PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s));
711:   PetscCall(PetscArraycpy(t->b, b, s));
712:   if (bembed) {
713:     PetscCall(PetscMalloc2(s, &t->bembed, s, &t->bembedt));
714:     PetscCall(PetscArraycpy(t->bembed, bembed, s));
715:   }
716:   for (i = 0; i < s; i++) {
717:     t->ASum[i]     = 0;
718:     t->GammaSum[i] = 0;
719:     for (j = 0; j < s; j++) {
720:       t->ASum[i] += A[i * s + j];
721:       t->GammaSum[i] += Gamma[i * s + j];
722:     }
723:   }
724:   PetscCall(PetscMalloc1(s * s, &GammaInv)); /* Need to use Scalar for inverse, then convert back to Real */
725:   for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i];
726:   for (i = 0; i < s; i++) {
727:     if (Gamma[i * s + i] == 0.0) {
728:       GammaInv[i * s + i] = 1.0;
729:       t->GammaZeroDiag[i] = PETSC_TRUE;
730:     } else {
731:       t->GammaZeroDiag[i] = PETSC_FALSE;
732:     }
733:   }

735:   switch (s) {
736:   case 1:
737:     GammaInv[0] = 1. / GammaInv[0];
738:     break;
739:   case 2:
740:     PetscCall(PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL));
741:     break;
742:   case 3:
743:     PetscCall(PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL));
744:     break;
745:   case 4:
746:     PetscCall(PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL));
747:     break;
748:   case 5: {
749:     PetscInt  ipvt5[5];
750:     MatScalar work5[5 * 5];
751:     PetscCall(PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL));
752:     break;
753:   }
754:   case 6:
755:     PetscCall(PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL));
756:     break;
757:   case 7:
758:     PetscCall(PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL));
759:     break;
760:   default:
761:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s);
762:   }
763:   for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
764:   PetscCall(PetscFree(GammaInv));

766:   for (i = 0; i < s; i++) {
767:     for (k = 0; k < i + 1; k++) {
768:       t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]);
769:       for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]);
770:     }
771:   }

773:   for (i = 0; i < s; i++) {
774:     for (j = 0; j < s; j++) {
775:       t->At[i * s + j] = 0;
776:       for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j];
777:     }
778:     t->bt[i] = 0;
779:     for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i];
780:     if (bembed) {
781:       t->bembedt[i] = 0;
782:       for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i];
783:     }
784:   }
785:   t->ccfl = 1.0; /* Fix this */

787:   t->pinterp = pinterp;
788:   PetscCall(PetscMalloc1(s * pinterp, &t->binterpt));
789:   PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp));
790:   link->next      = RosWTableauList;
791:   RosWTableauList = link;
792:   PetscFunctionReturn(PETSC_SUCCESS);
793: }

795: /*@C
796:   TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices

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

800:   Input Parameters:
801: + name  - identifier for method
802: . gamma - leading coefficient (diagonal entry)
803: . a2    - design parameter, see Table 7.2 of Hairer&Wanner
804: . a3    - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22)
805: . b3    - design parameter, see Table 7.2 of Hairer&Wanner
806: - e4    - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer

808:   Level: developer

810:   Notes:
811:   This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2.
812:   It is used here to implement several methods from the book and can be used to experiment with new methods.
813:   It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions.

815: .seealso: [](ch_ts), `TSRosW`, `TSRosWRegister()`
816: @*/
817: PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4)
818: {
819:   /* Declare numeric constants so they can be quad precision without being truncated at double */
820:   const PetscReal one = 1, two = 2, three = 3, four = 4, five = 5, six = 6, eight = 8, twelve = 12, twenty = 20, twentyfour = 24, p32 = one / six - gamma + gamma * gamma, p42 = one / eight - gamma / three, p43 = one / twelve - gamma / three, p44 = one / twentyfour - gamma / two + three / two * gamma * gamma - gamma * gamma * gamma, p56 = one / twenty - gamma / four;
821:   PetscReal   a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp;
822:   PetscReal   A[4][4], Gamma[4][4], b[4], bm[4];
823:   PetscScalar M[3][3], rhs[3];

825:   PetscFunctionBegin;
826:   /* Step 1: choose Gamma (input) */
827:   /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */
828:   if (a3 == (PetscReal)PETSC_DEFAULT) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */
829:   a4 = a3;                                                                                       /* consequence of 7.20 */

831:   /* Solve order conditions 7.15a, 7.15c, 7.15e */
832:   M[0][0] = one;
833:   M[0][1] = one;
834:   M[0][2] = one; /* 7.15a */
835:   M[1][0] = 0.0;
836:   M[1][1] = a2 * a2;
837:   M[1][2] = a4 * a4; /* 7.15c */
838:   M[2][0] = 0.0;
839:   M[2][1] = a2 * a2 * a2;
840:   M[2][2] = a4 * a4 * a4; /* 7.15e */
841:   rhs[0]  = one - b3;
842:   rhs[1]  = one / three - a3 * a3 * b3;
843:   rhs[2]  = one / four - a3 * a3 * a3 * b3;
844:   PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL));
845:   b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
846:   b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
847:   b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);

849:   /* Step 3 */
850:   beta43       = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */
851:   beta32beta2p = p44 / (b4 * beta43);                           /* 7.15h */
852:   beta4jbetajp = (p32 - b3 * beta32beta2p) / b4;
853:   M[0][0]      = b2;
854:   M[0][1]      = b3;
855:   M[0][2]      = b4;
856:   M[1][0]      = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp;
857:   M[1][1]      = a2 * a2 * beta4jbetajp;
858:   M[1][2]      = -a2 * a2 * beta32beta2p;
859:   M[2][0]      = b4 * beta43 * a3 * a3 - p43;
860:   M[2][1]      = -b4 * beta43 * a2 * a2;
861:   M[2][2]      = 0;
862:   rhs[0]       = one / two - gamma;
863:   rhs[1]       = 0;
864:   rhs[2]       = -a2 * a2 * p32;
865:   PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL));
866:   beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
867:   beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
868:   beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);

870:   /* Step 4: back-substitute */
871:   beta32 = beta32beta2p / beta2p;
872:   beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p;

874:   /* Step 5: 7.15f and 7.20, then 7.16 */
875:   a43 = 0;
876:   a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p);
877:   a42 = a32;

879:   A[0][0]     = 0;
880:   A[0][1]     = 0;
881:   A[0][2]     = 0;
882:   A[0][3]     = 0;
883:   A[1][0]     = a2;
884:   A[1][1]     = 0;
885:   A[1][2]     = 0;
886:   A[1][3]     = 0;
887:   A[2][0]     = a3 - a32;
888:   A[2][1]     = a32;
889:   A[2][2]     = 0;
890:   A[2][3]     = 0;
891:   A[3][0]     = a4 - a43 - a42;
892:   A[3][1]     = a42;
893:   A[3][2]     = a43;
894:   A[3][3]     = 0;
895:   Gamma[0][0] = gamma;
896:   Gamma[0][1] = 0;
897:   Gamma[0][2] = 0;
898:   Gamma[0][3] = 0;
899:   Gamma[1][0] = beta2p - A[1][0];
900:   Gamma[1][1] = gamma;
901:   Gamma[1][2] = 0;
902:   Gamma[1][3] = 0;
903:   Gamma[2][0] = beta3p - beta32 - A[2][0];
904:   Gamma[2][1] = beta32 - A[2][1];
905:   Gamma[2][2] = gamma;
906:   Gamma[2][3] = 0;
907:   Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0];
908:   Gamma[3][1] = beta42 - A[3][1];
909:   Gamma[3][2] = beta43 - A[3][2];
910:   Gamma[3][3] = gamma;
911:   b[0]        = b1;
912:   b[1]        = b2;
913:   b[2]        = b3;
914:   b[3]        = b4;

916:   /* Construct embedded formula using given e4. We are solving Equation 7.18. */
917:   bm[3] = b[3] - e4 * gamma;                                              /* using definition of E4 */
918:   bm[2] = (p32 - beta4jbetajp * bm[3]) / (beta32 * beta2p);               /* fourth row of 7.18 */
919:   bm[1] = (one / two - gamma - beta3p * bm[2] - beta4p * bm[3]) / beta2p; /* second row */
920:   bm[0] = one - bm[1] - bm[2] - bm[3];                                    /* first row */

922:   {
923:     const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three;
924:     PetscCheck(PetscAbs(misfit) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "Assumptions violated, could not construct a third order embedded method");
925:   }
926:   PetscCall(TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL));
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

930: /*
931:  The step completion formula is

933:  x1 = x0 + b^T Y

935:  where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been
936:  updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write

938:  x1e = x0 + be^T Y
939:      = x1 - b^T Y + be^T Y
940:      = x1 + (be - b)^T Y

942:  so we can evaluate the method of different order even after the step has been optimistically completed.
943: */
944: static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done)
945: {
946:   TS_RosW     *ros = (TS_RosW *)ts->data;
947:   RosWTableau  tab = ros->tableau;
948:   PetscScalar *w   = ros->work;
949:   PetscInt     i;

951:   PetscFunctionBegin;
952:   if (order == tab->order) {
953:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
954:       PetscCall(VecCopy(ts->vec_sol, U));
955:       for (i = 0; i < tab->s; i++) w[i] = tab->bt[i];
956:       PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
957:     } else PetscCall(VecCopy(ts->vec_sol, U));
958:     if (done) *done = PETSC_TRUE;
959:     PetscFunctionReturn(PETSC_SUCCESS);
960:   } else if (order == tab->order - 1) {
961:     if (!tab->bembedt) goto unavailable;
962:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
963:       PetscCall(VecCopy(ts->vec_sol, U));
964:       for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i];
965:       PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
966:     } else { /* Use rollback-and-recomplete formula (bembedt - bt) */
967:       for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
968:       PetscCall(VecCopy(ts->vec_sol, U));
969:       PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
970:     }
971:     if (done) *done = PETSC_TRUE;
972:     PetscFunctionReturn(PETSC_SUCCESS);
973:   }
974: unavailable:
975:   if (done) *done = PETSC_FALSE;
976:   else
977:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Rosenbrock-W '%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,
978:             tab->order, order);
979:   PetscFunctionReturn(PETSC_SUCCESS);
980: }

982: static PetscErrorCode TSStep_RosW(TS ts)
983: {
984:   TS_RosW         *ros = (TS_RosW *)ts->data;
985:   RosWTableau      tab = ros->tableau;
986:   const PetscInt   s   = tab->s;
987:   const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv;
988:   const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
989:   const PetscBool *GammaZeroDiag     = tab->GammaZeroDiag;
990:   PetscScalar     *w                 = ros->work;
991:   Vec             *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage;
992:   SNES             snes;
993:   TSAdapt          adapt;
994:   PetscInt         i, j, its, lits;
995:   PetscInt         rejections = 0;
996:   PetscBool        stageok, accept = PETSC_TRUE;
997:   PetscReal        next_time_step = ts->time_step;
998:   PetscInt         lag;

1000:   PetscFunctionBegin;
1001:   if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev));

1003:   ros->status = TS_STEP_INCOMPLETE;
1004:   while (!ts->reason && ros->status != TS_STEP_COMPLETE) {
1005:     const PetscReal h = ts->time_step;
1006:     for (i = 0; i < s; i++) {
1007:       ros->stage_time = ts->ptime + h * ASum[i];
1008:       PetscCall(TSPreStage(ts, ros->stage_time));
1009:       if (GammaZeroDiag[i]) {
1010:         ros->stage_explicit = PETSC_TRUE;
1011:         ros->scoeff         = 1.;
1012:       } else {
1013:         ros->stage_explicit = PETSC_FALSE;
1014:         ros->scoeff         = 1. / Gamma[i * s + i];
1015:       }

1017:       PetscCall(VecCopy(ts->vec_sol, Zstage));
1018:       for (j = 0; j < i; j++) w[j] = At[i * s + j];
1019:       PetscCall(VecMAXPY(Zstage, i, w, Y));

1021:       for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j];
1022:       PetscCall(VecZeroEntries(Zdot));
1023:       PetscCall(VecMAXPY(Zdot, i, w, Y));

1025:       /* Initial guess taken from last stage */
1026:       PetscCall(VecZeroEntries(Y[i]));

1028:       if (!ros->stage_explicit) {
1029:         PetscCall(TSGetSNES(ts, &snes));
1030:         if (!ros->recompute_jacobian && !i) {
1031:           PetscCall(SNESGetLagJacobian(snes, &lag));
1032:           if (lag == 1) {                            /* use did not set a nontrivial lag, so lag over all stages */
1033:             PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */
1034:           }
1035:         }
1036:         PetscCall(SNESSolve(snes, NULL, Y[i]));
1037:         if (!ros->recompute_jacobian && i == s - 1 && lag == 1) { PetscCall(SNESSetLagJacobian(snes, lag)); /* Set lag back to 1 so we know user did not set it */ }
1038:         PetscCall(SNESGetIterationNumber(snes, &its));
1039:         PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1040:         ts->snes_its += its;
1041:         ts->ksp_its += lits;
1042:       } else {
1043:         Mat J, Jp;
1044:         PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
1045:         PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE));
1046:         PetscCall(VecScale(Y[i], -1.0));
1047:         PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/

1049:         PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
1050:         for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j];
1051:         PetscCall(VecMAXPY(Zstage, i, w, Y));

1053:         /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
1054:         PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL));
1055:         PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE));
1056:         PetscCall(MatMult(J, Zstage, Zdot));
1057:         PetscCall(VecAXPY(Y[i], -1.0, Zdot));
1058:         ts->ksp_its += 1;

1060:         PetscCall(VecScale(Y[i], h));
1061:       }
1062:       PetscCall(TSPostStage(ts, ros->stage_time, i, Y));
1063:       PetscCall(TSGetAdapt(ts, &adapt));
1064:       PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok));
1065:       if (!stageok) goto reject_step;
1066:     }

1068:     ros->status = TS_STEP_INCOMPLETE;
1069:     PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL));
1070:     ros->status = TS_STEP_PENDING;
1071:     PetscCall(TSGetAdapt(ts, &adapt));
1072:     PetscCall(TSAdaptCandidatesClear(adapt));
1073:     PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
1074:     PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
1075:     ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1076:     if (!accept) { /* Roll back the current step */
1077:       PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
1078:       ts->time_step = next_time_step;
1079:       goto reject_step;
1080:     }

1082:     ts->ptime += ts->time_step;
1083:     ts->time_step = next_time_step;
1084:     break;

1086:   reject_step:
1087:     ts->reject++;
1088:     accept = PETSC_FALSE;
1089:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1090:       ts->reason = TS_DIVERGED_STEP_REJECTED;
1091:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
1092:     }
1093:   }
1094:   PetscFunctionReturn(PETSC_SUCCESS);
1095: }

1097: static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U)
1098: {
1099:   TS_RosW         *ros = (TS_RosW *)ts->data;
1100:   PetscInt         s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j;
1101:   PetscReal        h;
1102:   PetscReal        tt, t;
1103:   PetscScalar     *bt;
1104:   const PetscReal *Bt       = ros->tableau->binterpt;
1105:   const PetscReal *GammaInv = ros->tableau->GammaInv;
1106:   PetscScalar     *w        = ros->work;
1107:   Vec             *Y        = ros->Y;

1109:   PetscFunctionBegin;
1110:   PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name);

1112:   switch (ros->status) {
1113:   case TS_STEP_INCOMPLETE:
1114:   case TS_STEP_PENDING:
1115:     h = ts->time_step;
1116:     t = (itime - ts->ptime) / h;
1117:     break;
1118:   case TS_STEP_COMPLETE:
1119:     h = ts->ptime - ts->ptime_prev;
1120:     t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1121:     break;
1122:   default:
1123:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1124:   }
1125:   PetscCall(PetscMalloc1(s, &bt));
1126:   for (i = 0; i < s; i++) bt[i] = 0;
1127:   for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1128:     for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt;
1129:   }

1131:   /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
1132:   /* U <- 0*/
1133:   PetscCall(VecZeroEntries(U));
1134:   /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
1135:   for (j = 0; j < s; j++) w[j] = 0;
1136:   for (j = 0; j < s; j++) {
1137:     for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j];
1138:   }
1139:   PetscCall(VecMAXPY(U, i, w, Y));
1140:   /* U <- y(t) + U */
1141:   PetscCall(VecAXPY(U, 1, ros->vec_sol_prev));

1143:   PetscCall(PetscFree(bt));
1144:   PetscFunctionReturn(PETSC_SUCCESS);
1145: }

1147: /*------------------------------------------------------------*/

1149: static PetscErrorCode TSRosWTableauReset(TS ts)
1150: {
1151:   TS_RosW    *ros = (TS_RosW *)ts->data;
1152:   RosWTableau tab = ros->tableau;

1154:   PetscFunctionBegin;
1155:   if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1156:   PetscCall(VecDestroyVecs(tab->s, &ros->Y));
1157:   PetscCall(PetscFree(ros->work));
1158:   PetscFunctionReturn(PETSC_SUCCESS);
1159: }

1161: static PetscErrorCode TSReset_RosW(TS ts)
1162: {
1163:   TS_RosW *ros = (TS_RosW *)ts->data;

1165:   PetscFunctionBegin;
1166:   PetscCall(TSRosWTableauReset(ts));
1167:   PetscCall(VecDestroy(&ros->Ydot));
1168:   PetscCall(VecDestroy(&ros->Ystage));
1169:   PetscCall(VecDestroy(&ros->Zdot));
1170:   PetscCall(VecDestroy(&ros->Zstage));
1171:   PetscCall(VecDestroy(&ros->vec_sol_prev));
1172:   PetscFunctionReturn(PETSC_SUCCESS);
1173: }

1175: static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1176: {
1177:   TS_RosW *rw = (TS_RosW *)ts->data;

1179:   PetscFunctionBegin;
1180:   if (Ydot) {
1181:     if (dm && dm != ts->dm) {
1182:       PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1183:     } else *Ydot = rw->Ydot;
1184:   }
1185:   if (Zdot) {
1186:     if (dm && dm != ts->dm) {
1187:       PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1188:     } else *Zdot = rw->Zdot;
1189:   }
1190:   if (Ystage) {
1191:     if (dm && dm != ts->dm) {
1192:       PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1193:     } else *Ystage = rw->Ystage;
1194:   }
1195:   if (Zstage) {
1196:     if (dm && dm != ts->dm) {
1197:       PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1198:     } else *Zstage = rw->Zstage;
1199:   }
1200:   PetscFunctionReturn(PETSC_SUCCESS);
1201: }

1203: static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1204: {
1205:   PetscFunctionBegin;
1206:   if (Ydot) {
1207:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1208:   }
1209:   if (Zdot) {
1210:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1211:   }
1212:   if (Ystage) {
1213:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1214:   }
1215:   if (Zstage) {
1216:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1217:   }
1218:   PetscFunctionReturn(PETSC_SUCCESS);
1219: }

1221: static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx)
1222: {
1223:   PetscFunctionBegin;
1224:   PetscFunctionReturn(PETSC_SUCCESS);
1225: }

1227: static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1228: {
1229:   TS  ts = (TS)ctx;
1230:   Vec Ydot, Zdot, Ystage, Zstage;
1231:   Vec Ydotc, Zdotc, Ystagec, Zstagec;

1233:   PetscFunctionBegin;
1234:   PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1235:   PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1236:   PetscCall(MatRestrict(restrct, Ydot, Ydotc));
1237:   PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc));
1238:   PetscCall(MatRestrict(restrct, Ystage, Ystagec));
1239:   PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec));
1240:   PetscCall(MatRestrict(restrct, Zdot, Zdotc));
1241:   PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc));
1242:   PetscCall(MatRestrict(restrct, Zstage, Zstagec));
1243:   PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec));
1244:   PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1245:   PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1246:   PetscFunctionReturn(PETSC_SUCCESS);
1247: }

1249: static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx)
1250: {
1251:   PetscFunctionBegin;
1252:   PetscFunctionReturn(PETSC_SUCCESS);
1253: }

1255: static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1256: {
1257:   TS  ts = (TS)ctx;
1258:   Vec Ydot, Zdot, Ystage, Zstage;
1259:   Vec Ydots, Zdots, Ystages, Zstages;

1261:   PetscFunctionBegin;
1262:   PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1263:   PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));

1265:   PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));
1266:   PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));

1268:   PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));
1269:   PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));

1271:   PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));
1272:   PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));

1274:   PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));
1275:   PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));

1277:   PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1278:   PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));
1279:   PetscFunctionReturn(PETSC_SUCCESS);
1280: }

1282: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts)
1283: {
1284:   TS_RosW  *ros = (TS_RosW *)ts->data;
1285:   Vec       Ydot, Zdot, Ystage, Zstage;
1286:   PetscReal shift = ros->scoeff / ts->time_step;
1287:   DM        dm, dmsave;

1289:   PetscFunctionBegin;
1290:   PetscCall(SNESGetDM(snes, &dm));
1291:   PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1292:   PetscCall(VecWAXPY(Ydot, shift, U, Zdot));   /* Ydot = shift*U + Zdot */
1293:   PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */
1294:   dmsave = ts->dm;
1295:   ts->dm = dm;
1296:   PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE));
1297:   ts->dm = dmsave;
1298:   PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1299:   PetscFunctionReturn(PETSC_SUCCESS);
1300: }

1302: static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts)
1303: {
1304:   TS_RosW  *ros = (TS_RosW *)ts->data;
1305:   Vec       Ydot, Zdot, Ystage, Zstage;
1306:   PetscReal shift = ros->scoeff / ts->time_step;
1307:   DM        dm, dmsave;

1309:   PetscFunctionBegin;
1310:   /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1311:   PetscCall(SNESGetDM(snes, &dm));
1312:   PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1313:   dmsave = ts->dm;
1314:   ts->dm = dm;
1315:   PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE));
1316:   ts->dm = dmsave;
1317:   PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1318:   PetscFunctionReturn(PETSC_SUCCESS);
1319: }

1321: static PetscErrorCode TSRosWTableauSetUp(TS ts)
1322: {
1323:   TS_RosW    *ros = (TS_RosW *)ts->data;
1324:   RosWTableau tab = ros->tableau;

1326:   PetscFunctionBegin;
1327:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y));
1328:   PetscCall(PetscMalloc1(tab->s, &ros->work));
1329:   PetscFunctionReturn(PETSC_SUCCESS);
1330: }

1332: static PetscErrorCode TSSetUp_RosW(TS ts)
1333: {
1334:   TS_RosW         *ros = (TS_RosW *)ts->data;
1335:   DM               dm;
1336:   SNES             snes;
1337:   TSRHSJacobianFn *rhsjacobian;

1339:   PetscFunctionBegin;
1340:   PetscCall(TSRosWTableauSetUp(ts));
1341:   PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot));
1342:   PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage));
1343:   PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot));
1344:   PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage));
1345:   PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev));
1346:   PetscCall(TSGetDM(ts, &dm));
1347:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1348:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1349:   /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1350:   PetscCall(TSGetSNES(ts, &snes));
1351:   if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1352:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1353:   if (rhsjacobian == TSComputeRHSJacobianConstant) {
1354:     Mat Amat, Pmat;

1356:     /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */
1357:     PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
1358:     if (Amat && Amat == ts->Arhs) {
1359:       if (Amat == Pmat) {
1360:         PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1361:         PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL));
1362:       } else {
1363:         PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1364:         PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
1365:         if (Pmat && Pmat == ts->Brhs) {
1366:           PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
1367:           PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
1368:           PetscCall(MatDestroy(&Pmat));
1369:         }
1370:       }
1371:       PetscCall(MatDestroy(&Amat));
1372:     }
1373:   }
1374:   PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1376: /*------------------------------------------------------------*/

1378: static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems *PetscOptionsObject)
1379: {
1380:   TS_RosW *ros = (TS_RosW *)ts->data;
1381:   SNES     snes;

1383:   PetscFunctionBegin;
1384:   PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options");
1385:   {
1386:     RosWTableauLink link;
1387:     PetscInt        count, choice;
1388:     PetscBool       flg;
1389:     const char    **namelist;

1391:     for (link = RosWTableauList, count = 0; link; link = link->next, count++);
1392:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
1393:     for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1394:     PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg));
1395:     if (flg) PetscCall(TSRosWSetType(ts, namelist[choice]));
1396:     PetscCall(PetscFree(namelist));

1398:     PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL));
1399:   }
1400:   PetscOptionsHeadEnd();
1401:   /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1402:   PetscCall(TSGetSNES(ts, &snes));
1403:   if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1404:   PetscFunctionReturn(PETSC_SUCCESS);
1405: }

1407: static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer)
1408: {
1409:   TS_RosW  *ros = (TS_RosW *)ts->data;
1410:   PetscBool iascii;

1412:   PetscFunctionBegin;
1413:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1414:   if (iascii) {
1415:     RosWTableau tab = ros->tableau;
1416:     TSRosWType  rostype;
1417:     char        buf[512];
1418:     PetscInt    i;
1419:     PetscReal   abscissa[512];
1420:     PetscCall(TSRosWGetType(ts, &rostype));
1421:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Rosenbrock-W %s\n", rostype));
1422:     PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum));
1423:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Abscissa of A       = %s\n", buf));
1424:     for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i];
1425:     PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa));
1426:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Abscissa of A+Gamma = %s\n", buf));
1427:   }
1428:   PetscFunctionReturn(PETSC_SUCCESS);
1429: }

1431: static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer)
1432: {
1433:   SNES    snes;
1434:   TSAdapt adapt;

1436:   PetscFunctionBegin;
1437:   PetscCall(TSGetAdapt(ts, &adapt));
1438:   PetscCall(TSAdaptLoad(adapt, viewer));
1439:   PetscCall(TSGetSNES(ts, &snes));
1440:   PetscCall(SNESLoad(snes, viewer));
1441:   /* function and Jacobian context for SNES when used with TS is always ts object */
1442:   PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
1443:   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
1444:   PetscFunctionReturn(PETSC_SUCCESS);
1445: }

1447: /*@C
1448:   TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme

1450:   Logically Collective

1452:   Input Parameters:
1453: + ts       - timestepping context
1454: - roswtype - type of Rosenbrock-W scheme

1456:   Level: beginner

1458: .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3`
1459: @*/
1460: PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype)
1461: {
1462:   PetscFunctionBegin;
1464:   PetscAssertPointer(roswtype, 2);
1465:   PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype));
1466:   PetscFunctionReturn(PETSC_SUCCESS);
1467: }

1469: /*@C
1470:   TSRosWGetType - Get the type of Rosenbrock-W scheme

1472:   Logically Collective

1474:   Input Parameter:
1475: . ts - timestepping context

1477:   Output Parameter:
1478: . rostype - type of Rosenbrock-W scheme

1480:   Level: intermediate

1482: .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()`
1483: @*/
1484: PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype)
1485: {
1486:   PetscFunctionBegin;
1488:   PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype));
1489:   PetscFunctionReturn(PETSC_SUCCESS);
1490: }

1492: /*@C
1493:   TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step.

1495:   Logically Collective

1497:   Input Parameters:
1498: + ts  - timestepping context
1499: - flg - `PETSC_TRUE` to recompute the Jacobian at each stage

1501:   Level: intermediate

1503: .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()`
1504: @*/
1505: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg)
1506: {
1507:   PetscFunctionBegin;
1509:   PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg));
1510:   PetscFunctionReturn(PETSC_SUCCESS);
1511: }

1513: static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype)
1514: {
1515:   TS_RosW *ros = (TS_RosW *)ts->data;

1517:   PetscFunctionBegin;
1518:   *rostype = ros->tableau->name;
1519:   PetscFunctionReturn(PETSC_SUCCESS);
1520: }

1522: static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype)
1523: {
1524:   TS_RosW        *ros = (TS_RosW *)ts->data;
1525:   PetscBool       match;
1526:   RosWTableauLink link;

1528:   PetscFunctionBegin;
1529:   if (ros->tableau) {
1530:     PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match));
1531:     if (match) PetscFunctionReturn(PETSC_SUCCESS);
1532:   }
1533:   for (link = RosWTableauList; link; link = link->next) {
1534:     PetscCall(PetscStrcmp(link->tab.name, rostype, &match));
1535:     if (match) {
1536:       if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts));
1537:       ros->tableau = &link->tab;
1538:       if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts));
1539:       ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1540:       PetscFunctionReturn(PETSC_SUCCESS);
1541:     }
1542:   }
1543:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype);
1544: }

1546: static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg)
1547: {
1548:   TS_RosW *ros = (TS_RosW *)ts->data;

1550:   PetscFunctionBegin;
1551:   ros->recompute_jacobian = flg;
1552:   PetscFunctionReturn(PETSC_SUCCESS);
1553: }

1555: static PetscErrorCode TSDestroy_RosW(TS ts)
1556: {
1557:   PetscFunctionBegin;
1558:   PetscCall(TSReset_RosW(ts));
1559:   if (ts->dm) {
1560:     PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1561:     PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1562:   }
1563:   PetscCall(PetscFree(ts->data));
1564:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL));
1565:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL));
1566:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL));
1567:   PetscFunctionReturn(PETSC_SUCCESS);
1568: }

1570: /* ------------------------------------------------------------ */
1571: /*MC
1572:       TSROSW - ODE solver using Rosenbrock-W schemes

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

1578:   Level: beginner

1580:   Notes:
1581:   This method currently only works with autonomous ODE and DAE.

1583:   Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear.

1585:   Since this uses a single linear solve per time-step if you wish to lag the jacobian or preconditioner computation you must use also -snes_lag_jacobian_persists true or -snes_lag_jacobian_preconditioner true

1587:   Developer Notes:
1588:   Rosenbrock-W methods are typically specified for autonomous ODE

1590: $  udot = f(u)

1592:   by the stage equations

1594: $  k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j

1596:   and step completion formula

1598: $  u_1 = u_0 + sum_j b_j k_j

1600:   with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u)
1601:   and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner,
1602:   we define new variables for the stage equations

1604: $  y_i = gamma_ij k_j

1606:   The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define

1608: $  A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1}

1610:   to rewrite the method as

1612: .vb
1613:   [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j
1614:   u_1 = u_0 + sum_j bt_j y_j
1615: .ve

1617:    where we have introduced the mass matrix M. Continue by defining

1619: $  ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j

1621:    or, more compactly in tensor notation

1623: $  Ydot = 1/h (Gamma^{-1} \otimes I) Y .

1625:    Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current
1626:    stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the
1627:    equation

1629: $  g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0

1631:    with initial guess y_i = 0.

1633: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`,
1634:           `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType`
1635: M*/
1636: PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1637: {
1638:   TS_RosW *ros;

1640:   PetscFunctionBegin;
1641:   PetscCall(TSRosWInitializePackage());

1643:   ts->ops->reset          = TSReset_RosW;
1644:   ts->ops->destroy        = TSDestroy_RosW;
1645:   ts->ops->view           = TSView_RosW;
1646:   ts->ops->load           = TSLoad_RosW;
1647:   ts->ops->setup          = TSSetUp_RosW;
1648:   ts->ops->step           = TSStep_RosW;
1649:   ts->ops->interpolate    = TSInterpolate_RosW;
1650:   ts->ops->evaluatestep   = TSEvaluateStep_RosW;
1651:   ts->ops->setfromoptions = TSSetFromOptions_RosW;
1652:   ts->ops->snesfunction   = SNESTSFormFunction_RosW;
1653:   ts->ops->snesjacobian   = SNESTSFormJacobian_RosW;

1655:   ts->usessnes = PETSC_TRUE;

1657:   PetscCall(PetscNew(&ros));
1658:   ts->data = (void *)ros;

1660:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW));
1661:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW));
1662:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW));

1664:   PetscCall(TSRosWSetType(ts, TSRosWDefault));
1665:   PetscFunctionReturn(PETSC_SUCCESS);
1666: }