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: `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: `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: `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: `TSROSW`
104: M*/

106: /*MC
107:      TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1.

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:      References:
114: .  * - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005.

116:      Level: intermediate

118: .seealso: `TSROSW`
119: M*/

121: /*MC
122:      TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1.

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

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

128:      References:
129: .  * - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005.

131:      Level: intermediate

133: .seealso: `TSROSW`
134: M*/

136: /*MC
137:      TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme

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

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

143:      References:
144: .  * - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.

146:      Level: intermediate

148: .seealso: `TSROSW`, `TSROSWSANDU3`
149: M*/

151: /*MC
152:      TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme

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

156:      The third order method is L-stable, but not stiffly accurate.
157:      The second order embedded method is strongly A-stable with R(infty) = 0.5.
158:      The internal stages are L-stable.
159:      This method is called ROS3 in the paper.

161:      References:
162: .  * - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.

164:      Level: intermediate

166: .seealso: `TSROSW`, `TSROSWRODAS3`
167: M*/

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

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

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

176:      References:
177: . * - Emil Constantinescu

179:      Level: intermediate

181: .seealso: `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP`
182: M*/

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

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

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

191:      References:
192: . * - Emil Constantinescu

194:      Level: intermediate

196: .seealso: `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP`
197: M*/

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

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

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

206:      References:
207: . * - Emil Constantinescu

209:      Level: intermediate

211: .seealso: `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP`
212: M*/

214: /*MC
215:      TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop

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

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

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

223:      References:
224: +   * -  Kaps and Rentrop, Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979.
225: -   * -  Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.

227:      Hairer's code ros4.f

229:      Level: intermediate

231: .seealso: `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`
232: M*/

234: /*MC
235:      TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine

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

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

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

243:      References:
244: +   * -  Shampine, Implementation of Rosenbrock methods, 1982.
245: -   * -  Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.

247:      Hairer's code ros4.f

249:      Level: intermediate

251: .seealso: `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L`
252: M*/

254: /*MC
255:      TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen

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

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

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

263:      References:
264: +   * -  van Veldhuizen, D stability and Kaps Rentrop methods, 1984.
265: -   * -  Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.

267:      Hairer's code ros4.f

269:      Level: intermediate

271: .seealso: `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
272: M*/

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

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

279:      A-stable and L-stable

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

283:      References:
284: .  * -   Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.

286:      Hairer's code ros4.f

288:      Level: intermediate

290: .seealso: `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
291: M*/

293: /*@C
294:   TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in TSRosW

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

298:   Level: advanced

300: .seealso: `TSRosWRegisterDestroy()`
301: @*/
302: PetscErrorCode TSRosWRegisterAll(void)
303: {
304:   if (TSRosWRegisterAllCalled) return 0;
305:   TSRosWRegisterAllCalled = PETSC_TRUE;

307:   {
308:     const PetscReal A        = 0;
309:     const PetscReal Gamma    = 1;
310:     const PetscReal b        = 1;
311:     const PetscReal binterpt = 1;

313:     TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt);
314:   }

316:   {
317:     const PetscReal A        = 0;
318:     const PetscReal Gamma    = 0.5;
319:     const PetscReal b        = 1;
320:     const PetscReal binterpt = 1;

322:     TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt);
323:   }

325:   {
326:     /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0);   Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */
327:     const PetscReal A[2][2] =
328:       {
329:         {0,  0},
330:         {1., 0}
331:     },
332:                     Gamma[2][2] = {{1.707106781186547524401, 0}, {-2. * 1.707106781186547524401, 1.707106781186547524401}}, b[2] = {0.5, 0.5}, b1[2] = {1.0, 0.0};
333:     PetscReal binterpt[2][2];
334:     binterpt[0][0] = 1.707106781186547524401 - 1.0;
335:     binterpt[1][0] = 2.0 - 1.707106781186547524401;
336:     binterpt[0][1] = 1.707106781186547524401 - 1.5;
337:     binterpt[1][1] = 1.5 - 1.707106781186547524401;

339:     TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]);
340:   }
341:   {
342:     /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0);   Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */
343:     const PetscReal A[2][2] =
344:       {
345:         {0,  0},
346:         {1., 0}
347:     },
348:                     Gamma[2][2] = {{0.2928932188134524755992, 0}, {-2. * 0.2928932188134524755992, 0.2928932188134524755992}}, b[2] = {0.5, 0.5}, b1[2] = {1.0, 0.0};
349:     PetscReal binterpt[2][2];
350:     binterpt[0][0] = 0.2928932188134524755992 - 1.0;
351:     binterpt[1][0] = 2.0 - 0.2928932188134524755992;
352:     binterpt[0][1] = 0.2928932188134524755992 - 1.5;
353:     binterpt[1][1] = 1.5 - 0.2928932188134524755992;

355:     TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]);
356:   }
357:   {
358:     /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */
359:     PetscReal       binterpt[3][2];
360:     const PetscReal A[3][3] =
361:       {
362:         {0,                      0, 0},
363:         {1.5773502691896257e+00, 0, 0},
364:         {0.5,                    0, 0}
365:     },
366:                     Gamma[3][3] = {{7.8867513459481287e-01, 0, 0}, {-1.5773502691896257e+00, 7.8867513459481287e-01, 0}, {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01}}, b[3] = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01}, b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01};

368:     binterpt[0][0] = -0.8094010767585034;
369:     binterpt[1][0] = -0.5;
370:     binterpt[2][0] = 2.3094010767585034;
371:     binterpt[0][1] = 0.9641016151377548;
372:     binterpt[1][1] = 0.5;
373:     binterpt[2][1] = -1.4641016151377548;

375:     TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]);
376:   }
377:   {
378:     PetscReal binterpt[4][3];
379:     /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
380:     const PetscReal A[4][4] =
381:       {
382:         {0,                      0,                       0,  0},
383:         {8.7173304301691801e-01, 0,                       0,  0},
384:         {8.4457060015369423e-01, -1.1299064236484185e-01, 0,  0},
385:         {0,                      0,                       1., 0}
386:     },
387:                     Gamma[4][4] = {{4.3586652150845900e-01, 0, 0, 0}, {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0}, {-9.0338057013044082e-01, 5.4180672388095326e-02, 4.3586652150845900e-01, 0}, {2.4212380706095346e-01, -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01}}, b[4] = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01}, b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01};

389:     binterpt[0][0] = 1.0564298455794094;
390:     binterpt[1][0] = 2.296429974281067;
391:     binterpt[2][0] = -1.307599564525376;
392:     binterpt[3][0] = -1.045260255335102;
393:     binterpt[0][1] = -1.3864882699759573;
394:     binterpt[1][1] = -8.262611700275677;
395:     binterpt[2][1] = 7.250979895056055;
396:     binterpt[3][1] = 2.398120075195581;
397:     binterpt[0][2] = 0.5721822314575016;
398:     binterpt[1][2] = 4.742931142090097;
399:     binterpt[2][2] = -4.398120075195578;
400:     binterpt[3][2] = -0.9169932983520199;

402:     TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
403:   }
404:   {
405:     /* const PetscReal g = 0.5;       Directly written in-place below */
406:     const PetscReal A[4][4] =
407:       {
408:         {0,    0,     0,   0},
409:         {0,    0,     0,   0},
410:         {1.,   0,     0,   0},
411:         {0.75, -0.25, 0.5, 0}
412:     },
413:                     Gamma[4][4] = {{0.5, 0, 0, 0}, {1., 0.5, 0, 0}, {-0.25, -0.25, 0.5, 0}, {1. / 12, 1. / 12, -2. / 3, 0.5}}, b[4] = {5. / 6, -1. / 6, -1. / 6, 0.5}, b2[4] = {0.75, -0.25, 0.5, 0};

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

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

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

451:     binterpt[0][0] = 0.089316397477040902157517886164709;
452:     binterpt[1][0] = -0.91068360252295909784248211383529;
453:     binterpt[2][0] = 1.8213672050459181956849642276706;
454:     binterpt[0][1] = 0.077350269189625764509148780501957;
455:     binterpt[1][1] = 1.077350269189625764509148780502;
456:     binterpt[2][1] = -1.1547005383792515290182975610039;

458:     TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]);
459:   }

461:   {
462:     const PetscReal A[4][4] =
463:       {
464:         {0,       0,       0,       0},
465:         {1. / 2., 0,       0,       0},
466:         {1. / 2., 1. / 2., 0,       0},
467:         {1. / 6., 1. / 6., 1. / 6., 0}
468:     },
469:                     Gamma[4][4] = {{1. / 2., 0, 0, 0}, {0.0, 1. / 4., 0, 0}, {-2., -2. / 3., 2. / 3., 0}, {1. / 2., 5. / 36., -2. / 9, 0}}, b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}, b2[4] = {1. / 8., 3. / 4., 1. / 8., 0};
470:     PetscReal binterpt[4][3];

472:     binterpt[0][0] = 6.25;
473:     binterpt[1][0] = -30.25;
474:     binterpt[2][0] = 1.75;
475:     binterpt[3][0] = 23.25;
476:     binterpt[0][1] = -9.75;
477:     binterpt[1][1] = 58.75;
478:     binterpt[2][1] = -3.25;
479:     binterpt[3][1] = -45.75;
480:     binterpt[0][2] = 3.6666666666666666666666666666667;
481:     binterpt[1][2] = -28.333333333333333333333333333333;
482:     binterpt[2][2] = 1.6666666666666666666666666666667;
483:     binterpt[3][2] = 23.;

485:     TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
486:   }

488:   {
489:     const PetscReal A[4][4] =
490:       {
491:         {0,       0,       0,       0},
492:         {1. / 2., 0,       0,       0},
493:         {1. / 2., 1. / 2., 0,       0},
494:         {1. / 6., 1. / 6., 1. / 6., 0}
495:     },
496:                     Gamma[4][4] = {{1. / 2., 0, 0, 0}, {0.0, 3. / 4., 0, 0}, {-2. / 3., -23. / 9., 2. / 9., 0}, {1. / 18., 65. / 108., -2. / 27, 0}}, b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}, b2[4] = {3. / 16., 10. / 16., 3. / 16., 0};
497:     PetscReal binterpt[4][3];

499:     binterpt[0][0] = 1.6911764705882352941176470588235;
500:     binterpt[1][0] = 3.6813725490196078431372549019608;
501:     binterpt[2][0] = 0.23039215686274509803921568627451;
502:     binterpt[3][0] = -4.6029411764705882352941176470588;
503:     binterpt[0][1] = -0.95588235294117647058823529411765;
504:     binterpt[1][1] = -6.2401960784313725490196078431373;
505:     binterpt[2][1] = -0.31862745098039215686274509803922;
506:     binterpt[3][1] = 7.5147058823529411764705882352941;
507:     binterpt[0][2] = -0.56862745098039215686274509803922;
508:     binterpt[1][2] = 2.7254901960784313725490196078431;
509:     binterpt[2][2] = 0.25490196078431372549019607843137;
510:     binterpt[3][2] = -2.4117647058823529411764705882353;

512:     TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
513:   }

515:   {
516:     PetscReal A[4][4], Gamma[4][4], b[4], b2[4];
517:     PetscReal binterpt[4][3];

519:     Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816;
520:     Gamma[0][1] = 0;
521:     Gamma[0][2] = 0;
522:     Gamma[0][3] = 0;
523:     Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476;
524:     Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816;
525:     Gamma[1][2] = 0;
526:     Gamma[1][3] = 0;
527:     Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903;
528:     Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131;
529:     Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816;
530:     Gamma[2][3] = 0;
531:     Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783;
532:     Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984;
533:     Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198;
534:     Gamma[3][3] = 0;

536:     A[0][0] = 0;
537:     A[0][1] = 0;
538:     A[0][2] = 0;
539:     A[0][3] = 0;
540:     A[1][0] = 0.8717330430169179988320388950590125027645343373957631;
541:     A[1][1] = 0;
542:     A[1][2] = 0;
543:     A[1][3] = 0;
544:     A[2][0] = 0.5275890119763004115618079766722914408876108660811028;
545:     A[2][1] = 0.07241098802369958843819203208518599088698057726988732;
546:     A[2][2] = 0;
547:     A[2][3] = 0;
548:     A[3][0] = 0.3990960076760701320627260685975778145384666450351314;
549:     A[3][1] = -0.4375576546135194437228463747348862825846903771419953;
550:     A[3][2] = 1.038461646937449311660120300601880176655352737312713;
551:     A[3][3] = 0;

553:     b[0] = 0.1876410243467238251612921333138006734899663569186926;
554:     b[1] = -0.5952974735769549480478230473706443582188442040780541;
555:     b[2] = 0.9717899277217721234705114616271378792182450260943198;
556:     b[3] = 0.4358665215084589994160194475295062513822671686978816;

558:     b2[0] = 0.2147402862233891404862383521089097657790734483804460;
559:     b2[1] = -0.4851622638849390928209050538171743017757490232519684;
560:     b2[2] = 0.8687250025203875511662123688667549217531982787600080;
561:     b2[3] = 0.4016969751411624011684543450940068201770721128357014;

563:     binterpt[0][0] = 2.2565812720167954547104627844105;
564:     binterpt[1][0] = 1.349166413351089573796243820819;
565:     binterpt[2][0] = -2.4695174540533503758652847586647;
566:     binterpt[3][0] = -0.13623023131453465264142184656474;
567:     binterpt[0][1] = -3.0826699111559187902922463354557;
568:     binterpt[1][1] = -2.4689115685996042534544925650515;
569:     binterpt[2][1] = 5.7428279814696677152129332773553;
570:     binterpt[3][1] = -0.19124650171414467146619437684812;
571:     binterpt[0][2] = 1.0137296634858471607430756831148;
572:     binterpt[1][2] = 0.52444768167155973161042570784064;
573:     binterpt[2][2] = -2.3015205996945452158771370439586;
574:     binterpt[3][2] = 0.76334325453713832352363565300308;

576:     TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
577:   }
578:   TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_DEFAULT, PETSC_DEFAULT, 0, -0.1282612945269037e+01);
579:   TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_DEFAULT, PETSC_DEFAULT, 0, 125. / 108.);
580:   TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.355958941201148);
581:   TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.093502252409163);
582:   return 0;
583: }

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

588:    Not Collective

590:    Level: advanced

592: .seealso: `TSRosWRegister()`, `TSRosWRegisterAll()`
593: @*/
594: PetscErrorCode TSRosWRegisterDestroy(void)
595: {
596:   RosWTableauLink link;

598:   while ((link = RosWTableauList)) {
599:     RosWTableau t   = &link->tab;
600:     RosWTableauList = link->next;
601:     PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum);
602:     PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr);
603:     PetscFree2(t->bembed, t->bembedt);
604:     PetscFree(t->binterpt);
605:     PetscFree(t->name);
606:     PetscFree(link);
607:   }
608:   TSRosWRegisterAllCalled = PETSC_FALSE;
609:   return 0;
610: }

612: /*@C
613:   TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called
614:   from TSInitializePackage().

616:   Level: developer

618: .seealso: `PetscInitialize()`
619: @*/
620: PetscErrorCode TSRosWInitializePackage(void)
621: {
622:   if (TSRosWPackageInitialized) return 0;
623:   TSRosWPackageInitialized = PETSC_TRUE;
624:   TSRosWRegisterAll();
625:   PetscRegisterFinalize(TSRosWFinalizePackage);
626:   return 0;
627: }

629: /*@C
630:   TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is
631:   called from PetscFinalize().

633:   Level: developer

635: .seealso: `PetscFinalize()`
636: @*/
637: PetscErrorCode TSRosWFinalizePackage(void)
638: {
639:   TSRosWPackageInitialized = PETSC_FALSE;
640:   TSRosWRegisterDestroy();
641:   return 0;
642: }

644: /*@C
645:    TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation

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

649:    Input Parameters:
650: +  name - identifier for method
651: .  order - approximation order of method
652: .  s - number of stages, this is the dimension of the matrices below
653: .  A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular
654: .  Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal
655: .  b - Step completion table (dimension s)
656: .  bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available)
657: .  pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt
658: -  binterpt - Coefficients of the interpolation formula (dimension s*pinterp)

660:    Notes:
661:    Several Rosenbrock W methods are provided, this function is only needed to create new methods.

663:    Level: advanced

665: .seealso: `TSRosW`
666: @*/
667: 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[])
668: {
669:   RosWTableauLink link;
670:   RosWTableau     t;
671:   PetscInt        i, j, k;
672:   PetscScalar    *GammaInv;


680:   TSRosWInitializePackage();
681:   PetscNew(&link);
682:   t = &link->tab;
683:   PetscStrallocpy(name, &t->name);
684:   t->order = order;
685:   t->s     = s;
686:   PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum);
687:   PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr);
688:   PetscArraycpy(t->A, A, s * s);
689:   PetscArraycpy(t->Gamma, Gamma, s * s);
690:   PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s);
691:   PetscArraycpy(t->b, b, s);
692:   if (bembed) {
693:     PetscMalloc2(s, &t->bembed, s, &t->bembedt);
694:     PetscArraycpy(t->bembed, bembed, s);
695:   }
696:   for (i = 0; i < s; i++) {
697:     t->ASum[i]     = 0;
698:     t->GammaSum[i] = 0;
699:     for (j = 0; j < s; j++) {
700:       t->ASum[i] += A[i * s + j];
701:       t->GammaSum[i] += Gamma[i * s + j];
702:     }
703:   }
704:   PetscMalloc1(s * s, &GammaInv); /* Need to use Scalar for inverse, then convert back to Real */
705:   for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i];
706:   for (i = 0; i < s; i++) {
707:     if (Gamma[i * s + i] == 0.0) {
708:       GammaInv[i * s + i] = 1.0;
709:       t->GammaZeroDiag[i] = PETSC_TRUE;
710:     } else {
711:       t->GammaZeroDiag[i] = PETSC_FALSE;
712:     }
713:   }

715:   switch (s) {
716:   case 1:
717:     GammaInv[0] = 1. / GammaInv[0];
718:     break;
719:   case 2:
720:     PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL);
721:     break;
722:   case 3:
723:     PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL);
724:     break;
725:   case 4:
726:     PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL);
727:     break;
728:   case 5: {
729:     PetscInt  ipvt5[5];
730:     MatScalar work5[5 * 5];
731:     PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL);
732:     break;
733:   }
734:   case 6:
735:     PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL);
736:     break;
737:   case 7:
738:     PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL);
739:     break;
740:   default:
741:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s);
742:   }
743:   for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
744:   PetscFree(GammaInv);

746:   for (i = 0; i < s; i++) {
747:     for (k = 0; k < i + 1; k++) {
748:       t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]);
749:       for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]);
750:     }
751:   }

753:   for (i = 0; i < s; i++) {
754:     for (j = 0; j < s; j++) {
755:       t->At[i * s + j] = 0;
756:       for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j];
757:     }
758:     t->bt[i] = 0;
759:     for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i];
760:     if (bembed) {
761:       t->bembedt[i] = 0;
762:       for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i];
763:     }
764:   }
765:   t->ccfl = 1.0; /* Fix this */

767:   t->pinterp = pinterp;
768:   PetscMalloc1(s * pinterp, &t->binterpt);
769:   PetscArraycpy(t->binterpt, binterpt, s * pinterp);
770:   link->next      = RosWTableauList;
771:   RosWTableauList = link;
772:   return 0;
773: }

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

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

780:    Input Parameters:
781: +  name - identifier for method
782: .  gamma - leading coefficient (diagonal entry)
783: .  a2 - design parameter, see Table 7.2 of Hairer&Wanner
784: .  a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22)
785: .  b3 - design parameter, see Table 7.2 of Hairer&Wanner
786: .  beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner
787: -  e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer

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

794:    Level: developer

796: .seealso: `TSRosW`, `TSRosWRegister()`
797: @*/
798: PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4)
799: {
800:   /* Declare numeric constants so they can be quad precision without being truncated at double */
801:   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;
802:   PetscReal   a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp;
803:   PetscReal   A[4][4], Gamma[4][4], b[4], bm[4];
804:   PetscScalar M[3][3], rhs[3];

806:   /* Step 1: choose Gamma (input) */
807:   /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */
808:   if (a3 == PETSC_DEFAULT) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */
809:   a4 = a3;                                                                            /* consequence of 7.20 */

811:   /* Solve order conditions 7.15a, 7.15c, 7.15e */
812:   M[0][0] = one;
813:   M[0][1] = one;
814:   M[0][2] = one; /* 7.15a */
815:   M[1][0] = 0.0;
816:   M[1][1] = a2 * a2;
817:   M[1][2] = a4 * a4; /* 7.15c */
818:   M[2][0] = 0.0;
819:   M[2][1] = a2 * a2 * a2;
820:   M[2][2] = a4 * a4 * a4; /* 7.15e */
821:   rhs[0]  = one - b3;
822:   rhs[1]  = one / three - a3 * a3 * b3;
823:   rhs[2]  = one / four - a3 * a3 * a3 * b3;
824:   PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL);
825:   b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
826:   b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
827:   b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);

829:   /* Step 3 */
830:   beta43       = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */
831:   beta32beta2p = p44 / (b4 * beta43);                           /* 7.15h */
832:   beta4jbetajp = (p32 - b3 * beta32beta2p) / b4;
833:   M[0][0]      = b2;
834:   M[0][1]      = b3;
835:   M[0][2]      = b4;
836:   M[1][0]      = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp;
837:   M[1][1]      = a2 * a2 * beta4jbetajp;
838:   M[1][2]      = -a2 * a2 * beta32beta2p;
839:   M[2][0]      = b4 * beta43 * a3 * a3 - p43;
840:   M[2][1]      = -b4 * beta43 * a2 * a2;
841:   M[2][2]      = 0;
842:   rhs[0]       = one / two - gamma;
843:   rhs[1]       = 0;
844:   rhs[2]       = -a2 * a2 * p32;
845:   PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL);
846:   beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
847:   beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
848:   beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);

850:   /* Step 4: back-substitute */
851:   beta32 = beta32beta2p / beta2p;
852:   beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p;

854:   /* Step 5: 7.15f and 7.20, then 7.16 */
855:   a43 = 0;
856:   a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p);
857:   a42 = a32;

859:   A[0][0]     = 0;
860:   A[0][1]     = 0;
861:   A[0][2]     = 0;
862:   A[0][3]     = 0;
863:   A[1][0]     = a2;
864:   A[1][1]     = 0;
865:   A[1][2]     = 0;
866:   A[1][3]     = 0;
867:   A[2][0]     = a3 - a32;
868:   A[2][1]     = a32;
869:   A[2][2]     = 0;
870:   A[2][3]     = 0;
871:   A[3][0]     = a4 - a43 - a42;
872:   A[3][1]     = a42;
873:   A[3][2]     = a43;
874:   A[3][3]     = 0;
875:   Gamma[0][0] = gamma;
876:   Gamma[0][1] = 0;
877:   Gamma[0][2] = 0;
878:   Gamma[0][3] = 0;
879:   Gamma[1][0] = beta2p - A[1][0];
880:   Gamma[1][1] = gamma;
881:   Gamma[1][2] = 0;
882:   Gamma[1][3] = 0;
883:   Gamma[2][0] = beta3p - beta32 - A[2][0];
884:   Gamma[2][1] = beta32 - A[2][1];
885:   Gamma[2][2] = gamma;
886:   Gamma[2][3] = 0;
887:   Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0];
888:   Gamma[3][1] = beta42 - A[3][1];
889:   Gamma[3][2] = beta43 - A[3][2];
890:   Gamma[3][3] = gamma;
891:   b[0]        = b1;
892:   b[1]        = b2;
893:   b[2]        = b3;
894:   b[3]        = b4;

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

902:   {
903:     const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three;
905:   }
906:   TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL);
907:   return 0;
908: }

910: /*
911:  The step completion formula is

913:  x1 = x0 + b^T Y

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

918:  x1e = x0 + be^T Y
919:      = x1 - b^T Y + be^T Y
920:      = x1 + (be - b)^T Y

922:  so we can evaluate the method of different order even after the step has been optimistically completed.
923: */
924: static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done)
925: {
926:   TS_RosW     *ros = (TS_RosW *)ts->data;
927:   RosWTableau  tab = ros->tableau;
928:   PetscScalar *w   = ros->work;
929:   PetscInt     i;

931:   if (order == tab->order) {
932:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
933:       VecCopy(ts->vec_sol, U);
934:       for (i = 0; i < tab->s; i++) w[i] = tab->bt[i];
935:       VecMAXPY(U, tab->s, w, ros->Y);
936:     } else VecCopy(ts->vec_sol, U);
937:     if (done) *done = PETSC_TRUE;
938:     return 0;
939:   } else if (order == tab->order - 1) {
940:     if (!tab->bembedt) goto unavailable;
941:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
942:       VecCopy(ts->vec_sol, U);
943:       for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i];
944:       VecMAXPY(U, tab->s, w, ros->Y);
945:     } else { /* Use rollback-and-recomplete formula (bembedt - bt) */
946:       for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
947:       VecCopy(ts->vec_sol, U);
948:       VecMAXPY(U, tab->s, w, ros->Y);
949:     }
950:     if (done) *done = PETSC_TRUE;
951:     return 0;
952:   }
953: unavailable:
954:   if (done) *done = PETSC_FALSE;
955:   else
956:     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,
957:             tab->order, order);
958:   return 0;
959: }

961: static PetscErrorCode TSRollBack_RosW(TS ts)
962: {
963:   TS_RosW *ros = (TS_RosW *)ts->data;

965:   VecCopy(ros->vec_sol_prev, ts->vec_sol);
966:   return 0;
967: }

969: static PetscErrorCode TSStep_RosW(TS ts)
970: {
971:   TS_RosW         *ros = (TS_RosW *)ts->data;
972:   RosWTableau      tab = ros->tableau;
973:   const PetscInt   s   = tab->s;
974:   const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv;
975:   const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
976:   const PetscBool *GammaZeroDiag     = tab->GammaZeroDiag;
977:   PetscScalar     *w                 = ros->work;
978:   Vec             *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage;
979:   SNES             snes;
980:   TSAdapt          adapt;
981:   PetscInt         i, j, its, lits;
982:   PetscInt         rejections = 0;
983:   PetscBool        stageok, accept = PETSC_TRUE;
984:   PetscReal        next_time_step = ts->time_step;
985:   PetscInt         lag;

987:   if (!ts->steprollback) VecCopy(ts->vec_sol, ros->vec_sol_prev);

989:   ros->status = TS_STEP_INCOMPLETE;
990:   while (!ts->reason && ros->status != TS_STEP_COMPLETE) {
991:     const PetscReal h = ts->time_step;
992:     for (i = 0; i < s; i++) {
993:       ros->stage_time = ts->ptime + h * ASum[i];
994:       TSPreStage(ts, ros->stage_time);
995:       if (GammaZeroDiag[i]) {
996:         ros->stage_explicit = PETSC_TRUE;
997:         ros->scoeff         = 1.;
998:       } else {
999:         ros->stage_explicit = PETSC_FALSE;
1000:         ros->scoeff         = 1. / Gamma[i * s + i];
1001:       }

1003:       VecCopy(ts->vec_sol, Zstage);
1004:       for (j = 0; j < i; j++) w[j] = At[i * s + j];
1005:       VecMAXPY(Zstage, i, w, Y);

1007:       for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j];
1008:       VecZeroEntries(Zdot);
1009:       VecMAXPY(Zdot, i, w, Y);

1011:       /* Initial guess taken from last stage */
1012:       VecZeroEntries(Y[i]);

1014:       if (!ros->stage_explicit) {
1015:         TSGetSNES(ts, &snes);
1016:         if (!ros->recompute_jacobian && !i) {
1017:           SNESGetLagJacobian(snes, &lag);
1018:           if (lag == 1) {                            /* use did not set a nontrivial lag, so lag over all stages */
1019:             SNESSetLagJacobian(snes, -2); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */
1020:           }
1021:         }
1022:         SNESSolve(snes, NULL, Y[i]);
1023:         if (!ros->recompute_jacobian && i == s - 1 && lag == 1) { SNESSetLagJacobian(snes, lag); /* Set lag back to 1 so we know user did not set it */ }
1024:         SNESGetIterationNumber(snes, &its);
1025:         SNESGetLinearSolveIterations(snes, &lits);
1026:         ts->snes_its += its;
1027:         ts->ksp_its += lits;
1028:       } else {
1029:         Mat J, Jp;
1030:         VecZeroEntries(Ydot); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
1031:         TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE);
1032:         VecScale(Y[i], -1.0);
1033:         VecAXPY(Y[i], -1.0, Zdot); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/

1035:         VecZeroEntries(Zstage); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
1036:         for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j];
1037:         VecMAXPY(Zstage, i, w, Y);

1039:         /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
1040:         TSGetIJacobian(ts, &J, &Jp, NULL, NULL);
1041:         TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE);
1042:         MatMult(J, Zstage, Zdot);
1043:         VecAXPY(Y[i], -1.0, Zdot);
1044:         ts->ksp_its += 1;

1046:         VecScale(Y[i], h);
1047:       }
1048:       TSPostStage(ts, ros->stage_time, i, Y);
1049:       TSGetAdapt(ts, &adapt);
1050:       TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok);
1051:       if (!stageok) goto reject_step;
1052:     }

1054:     ros->status = TS_STEP_INCOMPLETE;
1055:     TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL);
1056:     ros->status = TS_STEP_PENDING;
1057:     TSGetAdapt(ts, &adapt);
1058:     TSAdaptCandidatesClear(adapt);
1059:     TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE);
1060:     TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
1061:     ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1062:     if (!accept) { /* Roll back the current step */
1063:       TSRollBack_RosW(ts);
1064:       ts->time_step = next_time_step;
1065:       goto reject_step;
1066:     }

1068:     ts->ptime += ts->time_step;
1069:     ts->time_step = next_time_step;
1070:     break;

1072:   reject_step:
1073:     ts->reject++;
1074:     accept = PETSC_FALSE;
1075:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1076:       ts->reason = TS_DIVERGED_STEP_REJECTED;
1077:       PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
1078:     }
1079:   }
1080:   return 0;
1081: }

1083: static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U)
1084: {
1085:   TS_RosW         *ros = (TS_RosW *)ts->data;
1086:   PetscInt         s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j;
1087:   PetscReal        h;
1088:   PetscReal        tt, t;
1089:   PetscScalar     *bt;
1090:   const PetscReal *Bt       = ros->tableau->binterpt;
1091:   const PetscReal *GammaInv = ros->tableau->GammaInv;
1092:   PetscScalar     *w        = ros->work;
1093:   Vec             *Y        = ros->Y;


1097:   switch (ros->status) {
1098:   case TS_STEP_INCOMPLETE:
1099:   case TS_STEP_PENDING:
1100:     h = ts->time_step;
1101:     t = (itime - ts->ptime) / h;
1102:     break;
1103:   case TS_STEP_COMPLETE:
1104:     h = ts->ptime - ts->ptime_prev;
1105:     t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1106:     break;
1107:   default:
1108:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1109:   }
1110:   PetscMalloc1(s, &bt);
1111:   for (i = 0; i < s; i++) bt[i] = 0;
1112:   for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1113:     for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt;
1114:   }

1116:   /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
1117:   /* U <- 0*/
1118:   VecZeroEntries(U);
1119:   /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
1120:   for (j = 0; j < s; j++) w[j] = 0;
1121:   for (j = 0; j < s; j++) {
1122:     for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j];
1123:   }
1124:   VecMAXPY(U, i, w, Y);
1125:   /* U <- y(t) + U */
1126:   VecAXPY(U, 1, ros->vec_sol_prev);

1128:   PetscFree(bt);
1129:   return 0;
1130: }

1132: /*------------------------------------------------------------*/

1134: static PetscErrorCode TSRosWTableauReset(TS ts)
1135: {
1136:   TS_RosW    *ros = (TS_RosW *)ts->data;
1137:   RosWTableau tab = ros->tableau;

1139:   if (!tab) return 0;
1140:   VecDestroyVecs(tab->s, &ros->Y);
1141:   PetscFree(ros->work);
1142:   return 0;
1143: }

1145: static PetscErrorCode TSReset_RosW(TS ts)
1146: {
1147:   TS_RosW *ros = (TS_RosW *)ts->data;

1149:   TSRosWTableauReset(ts);
1150:   VecDestroy(&ros->Ydot);
1151:   VecDestroy(&ros->Ystage);
1152:   VecDestroy(&ros->Zdot);
1153:   VecDestroy(&ros->Zstage);
1154:   VecDestroy(&ros->vec_sol_prev);
1155:   return 0;
1156: }

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

1162:   if (Ydot) {
1163:     if (dm && dm != ts->dm) {
1164:       DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot);
1165:     } else *Ydot = rw->Ydot;
1166:   }
1167:   if (Zdot) {
1168:     if (dm && dm != ts->dm) {
1169:       DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot);
1170:     } else *Zdot = rw->Zdot;
1171:   }
1172:   if (Ystage) {
1173:     if (dm && dm != ts->dm) {
1174:       DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage);
1175:     } else *Ystage = rw->Ystage;
1176:   }
1177:   if (Zstage) {
1178:     if (dm && dm != ts->dm) {
1179:       DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage);
1180:     } else *Zstage = rw->Zstage;
1181:   }
1182:   return 0;
1183: }

1185: static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1186: {
1187:   if (Ydot) {
1188:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot);
1189:   }
1190:   if (Zdot) {
1191:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot);
1192:   }
1193:   if (Ystage) {
1194:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage);
1195:   }
1196:   if (Zstage) {
1197:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage);
1198:   }
1199:   return 0;
1200: }

1202: static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx)
1203: {
1204:   return 0;
1205: }

1207: static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1208: {
1209:   TS  ts = (TS)ctx;
1210:   Vec Ydot, Zdot, Ystage, Zstage;
1211:   Vec Ydotc, Zdotc, Ystagec, Zstagec;

1213:   TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage);
1214:   TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec);
1215:   MatRestrict(restrct, Ydot, Ydotc);
1216:   VecPointwiseMult(Ydotc, rscale, Ydotc);
1217:   MatRestrict(restrct, Ystage, Ystagec);
1218:   VecPointwiseMult(Ystagec, rscale, Ystagec);
1219:   MatRestrict(restrct, Zdot, Zdotc);
1220:   VecPointwiseMult(Zdotc, rscale, Zdotc);
1221:   MatRestrict(restrct, Zstage, Zstagec);
1222:   VecPointwiseMult(Zstagec, rscale, Zstagec);
1223:   TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage);
1224:   TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec);
1225:   return 0;
1226: }

1228: static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx)
1229: {
1230:   return 0;
1231: }

1233: static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1234: {
1235:   TS  ts = (TS)ctx;
1236:   Vec Ydot, Zdot, Ystage, Zstage;
1237:   Vec Ydots, Zdots, Ystages, Zstages;

1239:   TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage);
1240:   TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages);

1242:   VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD);
1243:   VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD);

1245:   VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD);
1246:   VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD);

1248:   VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD);
1249:   VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD);

1251:   VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD);
1252:   VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD);

1254:   TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage);
1255:   TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages);
1256:   return 0;
1257: }

1259: /*
1260:   This defines the nonlinear equation that is to be solved with SNES
1261:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
1262: */
1263: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts)
1264: {
1265:   TS_RosW  *ros = (TS_RosW *)ts->data;
1266:   Vec       Ydot, Zdot, Ystage, Zstage;
1267:   PetscReal shift = ros->scoeff / ts->time_step;
1268:   DM        dm, dmsave;

1270:   SNESGetDM(snes, &dm);
1271:   TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1272:   VecWAXPY(Ydot, shift, U, Zdot);   /* Ydot = shift*U + Zdot */
1273:   VecWAXPY(Ystage, 1.0, U, Zstage); /* Ystage = U + Zstage */
1274:   dmsave = ts->dm;
1275:   ts->dm = dm;
1276:   TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE);
1277:   ts->dm = dmsave;
1278:   TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1279:   return 0;
1280: }

1282: static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, 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:   /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1290:   SNESGetDM(snes, &dm);
1291:   TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1292:   dmsave = ts->dm;
1293:   ts->dm = dm;
1294:   TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE);
1295:   ts->dm = dmsave;
1296:   TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1297:   return 0;
1298: }

1300: static PetscErrorCode TSRosWTableauSetUp(TS ts)
1301: {
1302:   TS_RosW    *ros = (TS_RosW *)ts->data;
1303:   RosWTableau tab = ros->tableau;

1305:   VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y);
1306:   PetscMalloc1(tab->s, &ros->work);
1307:   return 0;
1308: }

1310: static PetscErrorCode TSSetUp_RosW(TS ts)
1311: {
1312:   TS_RosW      *ros = (TS_RosW *)ts->data;
1313:   DM            dm;
1314:   SNES          snes;
1315:   TSRHSJacobian rhsjacobian;

1317:   TSRosWTableauSetUp(ts);
1318:   VecDuplicate(ts->vec_sol, &ros->Ydot);
1319:   VecDuplicate(ts->vec_sol, &ros->Ystage);
1320:   VecDuplicate(ts->vec_sol, &ros->Zdot);
1321:   VecDuplicate(ts->vec_sol, &ros->Zstage);
1322:   VecDuplicate(ts->vec_sol, &ros->vec_sol_prev);
1323:   TSGetDM(ts, &dm);
1324:   DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts);
1325:   DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts);
1326:   /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1327:   TSGetSNES(ts, &snes);
1328:   if (!((PetscObject)snes)->type_name) SNESSetType(snes, SNESKSPONLY);
1329:   DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);
1330:   if (rhsjacobian == TSComputeRHSJacobianConstant) {
1331:     Mat Amat, Pmat;

1333:     /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */
1334:     SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL);
1335:     if (Amat && Amat == ts->Arhs) {
1336:       if (Amat == Pmat) {
1337:         MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
1338:         SNESSetJacobian(snes, Amat, Amat, NULL, NULL);
1339:       } else {
1340:         MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
1341:         SNESSetJacobian(snes, Amat, NULL, NULL, NULL);
1342:         if (Pmat && Pmat == ts->Brhs) {
1343:           MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat);
1344:           SNESSetJacobian(snes, NULL, Pmat, NULL, NULL);
1345:           MatDestroy(&Pmat);
1346:         }
1347:       }
1348:       MatDestroy(&Amat);
1349:     }
1350:   }
1351:   return 0;
1352: }
1353: /*------------------------------------------------------------*/

1355: static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems *PetscOptionsObject)
1356: {
1357:   TS_RosW *ros = (TS_RosW *)ts->data;
1358:   SNES     snes;

1360:   PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options");
1361:   {
1362:     RosWTableauLink link;
1363:     PetscInt        count, choice;
1364:     PetscBool       flg;
1365:     const char    **namelist;

1367:     for (link = RosWTableauList, count = 0; link; link = link->next, count++)
1368:       ;
1369:     PetscMalloc1(count, (char ***)&namelist);
1370:     for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1371:     PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg);
1372:     if (flg) TSRosWSetType(ts, namelist[choice]);
1373:     PetscFree(namelist);

1375:     PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL);
1376:   }
1377:   PetscOptionsHeadEnd();
1378:   /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1379:   TSGetSNES(ts, &snes);
1380:   if (!((PetscObject)snes)->type_name) SNESSetType(snes, SNESKSPONLY);
1381:   return 0;
1382: }

1384: static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer)
1385: {
1386:   TS_RosW  *ros = (TS_RosW *)ts->data;
1387:   PetscBool iascii;

1389:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1390:   if (iascii) {
1391:     RosWTableau tab = ros->tableau;
1392:     TSRosWType  rostype;
1393:     char        buf[512];
1394:     PetscInt    i;
1395:     PetscReal   abscissa[512];
1396:     TSRosWGetType(ts, &rostype);
1397:     PetscViewerASCIIPrintf(viewer, "  Rosenbrock-W %s\n", rostype);
1398:     PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum);
1399:     PetscViewerASCIIPrintf(viewer, "  Abscissa of A       = %s\n", buf);
1400:     for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i];
1401:     PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa);
1402:     PetscViewerASCIIPrintf(viewer, "  Abscissa of A+Gamma = %s\n", buf);
1403:   }
1404:   return 0;
1405: }

1407: static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer)
1408: {
1409:   SNES    snes;
1410:   TSAdapt adapt;

1412:   TSGetAdapt(ts, &adapt);
1413:   TSAdaptLoad(adapt, viewer);
1414:   TSGetSNES(ts, &snes);
1415:   SNESLoad(snes, viewer);
1416:   /* function and Jacobian context for SNES when used with TS is always ts object */
1417:   SNESSetFunction(snes, NULL, NULL, ts);
1418:   SNESSetJacobian(snes, NULL, NULL, NULL, ts);
1419:   return 0;
1420: }

1422: /*@C
1423:   TSRosWSetType - Set the type of Rosenbrock-W scheme

1425:   Logically collective

1427:   Input Parameters:
1428: +  ts - timestepping context
1429: -  roswtype - type of Rosenbrock-W scheme

1431:   Level: beginner

1433: .seealso: `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3`
1434: @*/
1435: PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype)
1436: {
1439:   PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype));
1440:   return 0;
1441: }

1443: /*@C
1444:   TSRosWGetType - Get the type of Rosenbrock-W scheme

1446:   Logically collective

1448:   Input Parameter:
1449: .  ts - timestepping context

1451:   Output Parameter:
1452: .  rostype - type of Rosenbrock-W scheme

1454:   Level: intermediate

1456: .seealso: `TSRosWGetType()`
1457: @*/
1458: PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype)
1459: {
1461:   PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype));
1462:   return 0;
1463: }

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

1468:   Logically collective

1470:   Input Parameters:
1471: +  ts - timestepping context
1472: -  flg - PETSC_TRUE to recompute the Jacobian at each stage

1474:   Level: intermediate

1476: .seealso: `TSRosWGetType()`
1477: @*/
1478: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg)
1479: {
1481:   PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg));
1482:   return 0;
1483: }

1485: static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype)
1486: {
1487:   TS_RosW *ros = (TS_RosW *)ts->data;

1489:   *rostype = ros->tableau->name;
1490:   return 0;
1491: }

1493: static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype)
1494: {
1495:   TS_RosW        *ros = (TS_RosW *)ts->data;
1496:   PetscBool       match;
1497:   RosWTableauLink link;

1499:   if (ros->tableau) {
1500:     PetscStrcmp(ros->tableau->name, rostype, &match);
1501:     if (match) return 0;
1502:   }
1503:   for (link = RosWTableauList; link; link = link->next) {
1504:     PetscStrcmp(link->tab.name, rostype, &match);
1505:     if (match) {
1506:       if (ts->setupcalled) TSRosWTableauReset(ts);
1507:       ros->tableau = &link->tab;
1508:       if (ts->setupcalled) TSRosWTableauSetUp(ts);
1509:       ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1510:       return 0;
1511:     }
1512:   }
1513:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype);
1514: }

1516: static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg)
1517: {
1518:   TS_RosW *ros = (TS_RosW *)ts->data;

1520:   ros->recompute_jacobian = flg;
1521:   return 0;
1522: }

1524: static PetscErrorCode TSDestroy_RosW(TS ts)
1525: {
1526:   TSReset_RosW(ts);
1527:   if (ts->dm) {
1528:     DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts);
1529:     DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts);
1530:   }
1531:   PetscFree(ts->data);
1532:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL);
1533:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL);
1534:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL);
1535:   return 0;
1536: }

1538: /* ------------------------------------------------------------ */
1539: /*MC
1540:       TSROSW - ODE solver using Rosenbrock-W schemes

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

1546:   Notes:
1547:   This method currently only works with autonomous ODE and DAE.

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

1551:   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

1553:   Developer Notes:
1554:   Rosenbrock-W methods are typically specified for autonomous ODE

1556: $  udot = f(u)

1558:   by the stage equations

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

1562:   and step completion formula

1564: $  u_1 = u_0 + sum_j b_j k_j

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

1570: $  y_i = gamma_ij k_j

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

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

1576:   to rewrite the method as

1578: $  [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j
1579: $  u_1 = u_0 + sum_j bt_j y_j

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

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

1585:    or, more compactly in tensor notation

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

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

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

1595:    with initial guess y_i = 0.

1597:   Level: beginner

1599: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`,
1600:           `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`
1601: M*/
1602: PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1603: {
1604:   TS_RosW *ros;

1606:   TSRosWInitializePackage();

1608:   ts->ops->reset          = TSReset_RosW;
1609:   ts->ops->destroy        = TSDestroy_RosW;
1610:   ts->ops->view           = TSView_RosW;
1611:   ts->ops->load           = TSLoad_RosW;
1612:   ts->ops->setup          = TSSetUp_RosW;
1613:   ts->ops->step           = TSStep_RosW;
1614:   ts->ops->interpolate    = TSInterpolate_RosW;
1615:   ts->ops->evaluatestep   = TSEvaluateStep_RosW;
1616:   ts->ops->rollback       = TSRollBack_RosW;
1617:   ts->ops->setfromoptions = TSSetFromOptions_RosW;
1618:   ts->ops->snesfunction   = SNESTSFormFunction_RosW;
1619:   ts->ops->snesjacobian   = SNESTSFormJacobian_RosW;

1621:   ts->usessnes = PETSC_TRUE;

1623:   PetscNew(&ros);
1624:   ts->data = (void *)ros;

1626:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW);
1627:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW);
1628:   PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW);

1630:   TSRosWSetType(ts, TSRosWDefault);
1631:   return 0;
1632: }