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: .  1. -   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: .  1. -   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: .  1. -   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: .  1. -   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: +   1. -  Kaps and Rentrop, Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979.
225: -   2. -  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: +   1. -  Shampine, Implementation of Rosenbrock methods, 1982.
245: -   2. -  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: +   1. -  van Veldhuizen, D stability and Kaps Rentrop methods, 1984.
265: -   2. -  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: .  1. -   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: {

307:   if (TSRosWRegisterAllCalled) return(0);
308:   TSRosWRegisterAllCalled = PETSC_TRUE;

310:   {
311:     const PetscReal A = 0;
312:     const PetscReal Gamma = 1;
313:     const PetscReal b = 1;
314:     const PetscReal binterpt=1;

316:     TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,NULL,1,&binterpt);
317:   }

319:   {
320:     const PetscReal A = 0;
321:     const PetscReal Gamma = 0.5;
322:     const PetscReal b = 1;
323:     const PetscReal binterpt=1;

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

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

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

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

371:       binterpt[0][0] = -0.8094010767585034;
372:       binterpt[1][0] = -0.5;
373:       binterpt[2][0] = 2.3094010767585034;
374:       binterpt[0][1] = 0.9641016151377548;
375:       binterpt[1][1] = 0.5;
376:       binterpt[2][1] = -1.4641016151377548;

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

395:     binterpt[0][0]=1.0564298455794094;
396:     binterpt[1][0]=2.296429974281067;
397:     binterpt[2][0]=-1.307599564525376;
398:     binterpt[3][0]=-1.045260255335102;
399:     binterpt[0][1]=-1.3864882699759573;
400:     binterpt[1][1]=-8.262611700275677;
401:     binterpt[2][1]=7.250979895056055;
402:     binterpt[3][1]=2.398120075195581;
403:     binterpt[0][2]=0.5721822314575016;
404:     binterpt[1][2]=4.742931142090097;
405:     binterpt[2][2]=-4.398120075195578;
406:     binterpt[3][2]=-0.9169932983520199;

408:     TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
409:   }
410:   {
411:     /* const PetscReal g = 0.5;       Directly written in-place below */
412:     const PetscReal
413:       A[4][4] = {{0,0,0,0},
414:                  {0,0,0,0},
415:                  {1.,0,0,0},
416:                  {0.75,-0.25,0.5,0}},
417:       Gamma[4][4] = {{0.5,0,0,0},
418:                      {1.,0.5,0,0},
419:                      {-0.25,-0.25,0.5,0},
420:                      {1./12,1./12,-2./3,0.5}},
421:       b[4]  = {5./6,-1./6,-1./6,0.5},
422:       b2[4] = {0.75,-0.25,0.5,0};

424:     TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,NULL);
425:   }
426:   {
427:     /*const PetscReal g = 0.43586652150845899941601945119356;       Directly written in-place below */
428:     const PetscReal
429:       A[3][3] = {{0,0,0},
430:                  {0.43586652150845899941601945119356,0,0},
431:                  {0.43586652150845899941601945119356,0,0}},
432:       Gamma[3][3] = {{0.43586652150845899941601945119356,0,0},
433:                      {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0},
434:                      {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}},
435:       b[3]  = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829},
436:       b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619};

438:     PetscReal binterpt[3][2];
439:     binterpt[0][0] = 3.793692883777660870425141387941;
440:     binterpt[1][0] = -2.918692883777660870425141387941;
441:     binterpt[2][0] = 0.125;
442:     binterpt[0][1] = -0.725741064379812106687651020584;
443:     binterpt[1][1] = 0.559074397713145440020984353917;
444:     binterpt[2][1] = 0.16666666666666666666666666666667;

446:     TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);
447:   }
448:   {
449:     /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0;
450:      * Direct evaluation: s3 = 1.732050807568877293527;
451:      *                     g = 0.7886751345948128822546;
452:      * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */
453:     const PetscReal
454:       A[3][3] = {{0,0,0},
455:                  {1,0,0},
456:                  {0.25,0.25,0}},
457:       Gamma[3][3] = {{0,0,0},
458:                      {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0},
459:                      {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}},
460:       b[3]  = {1./6.,1./6.,2./3.},
461:       b2[3] = {1./4.,1./4.,1./2.};
462:     PetscReal binterpt[3][2];

464:     binterpt[0][0]=0.089316397477040902157517886164709;
465:     binterpt[1][0]=-0.91068360252295909784248211383529;
466:     binterpt[2][0]=1.8213672050459181956849642276706;
467:     binterpt[0][1]=0.077350269189625764509148780501957;
468:     binterpt[1][1]=1.077350269189625764509148780502;
469:     binterpt[2][1]=-1.1547005383792515290182975610039;

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

474:   {
475:     const PetscReal
476:       A[4][4] = {{0,0,0,0},
477:                  {1./2.,0,0,0},
478:                  {1./2.,1./2.,0,0},
479:                  {1./6.,1./6.,1./6.,0}},
480:       Gamma[4][4] = {{1./2.,0,0,0},
481:                      {0.0,1./4.,0,0},
482:                      {-2.,-2./3.,2./3.,0},
483:                      {1./2.,5./36.,-2./9,0}},
484:       b[4]  = {1./6.,1./6.,1./6.,1./2.},
485:       b2[4] = {1./8.,3./4.,1./8.,0};
486:     PetscReal binterpt[4][3];

488:     binterpt[0][0]=6.25;
489:     binterpt[1][0]=-30.25;
490:     binterpt[2][0]=1.75;
491:     binterpt[3][0]=23.25;
492:     binterpt[0][1]=-9.75;
493:     binterpt[1][1]=58.75;
494:     binterpt[2][1]=-3.25;
495:     binterpt[3][1]=-45.75;
496:     binterpt[0][2]=3.6666666666666666666666666666667;
497:     binterpt[1][2]=-28.333333333333333333333333333333;
498:     binterpt[2][2]=1.6666666666666666666666666666667;
499:     binterpt[3][2]=23.;

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

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

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

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

534:   {
535:     PetscReal A[4][4],Gamma[4][4],b[4],b2[4];
536:     PetscReal binterpt[4][3];

538:     Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816;
539:     Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0;
540:     Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476;
541:     Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816;
542:     Gamma[1][2]=0; 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; A[0][1]=0; A[0][2]=0; A[0][3]=0;
553:     A[1][0]=0.8717330430169179988320388950590125027645343373957631;
554:     A[1][1]=0; A[1][2]=0; A[1][3]=0;
555:     A[2][0]=0.5275890119763004115618079766722914408876108660811028;
556:     A[2][1]=0.07241098802369958843819203208518599088698057726988732;
557:     A[2][2]=0; A[2][3]=0;
558:     A[3][0]=0.3990960076760701320627260685975778145384666450351314;
559:     A[3][1]=-0.4375576546135194437228463747348862825846903771419953;
560:     A[3][2]=1.038461646937449311660120300601880176655352737312713;
561:     A[3][3]=0;

563:     b[0]=0.1876410243467238251612921333138006734899663569186926;
564:     b[1]=-0.5952974735769549480478230473706443582188442040780541;
565:     b[2]=0.9717899277217721234705114616271378792182450260943198;
566:     b[3]=0.4358665215084589994160194475295062513822671686978816;

568:     b2[0]=0.2147402862233891404862383521089097657790734483804460;
569:     b2[1]=-0.4851622638849390928209050538171743017757490232519684;
570:     b2[2]=0.8687250025203875511662123688667549217531982787600080;
571:     b2[3]=0.4016969751411624011684543450940068201770721128357014;

573:     binterpt[0][0]=2.2565812720167954547104627844105;
574:     binterpt[1][0]=1.349166413351089573796243820819;
575:     binterpt[2][0]=-2.4695174540533503758652847586647;
576:     binterpt[3][0]=-0.13623023131453465264142184656474;
577:     binterpt[0][1]=-3.0826699111559187902922463354557;
578:     binterpt[1][1]=-2.4689115685996042534544925650515;
579:     binterpt[2][1]=5.7428279814696677152129332773553;
580:     binterpt[3][1]=-0.19124650171414467146619437684812;
581:     binterpt[0][2]=1.0137296634858471607430756831148;
582:     binterpt[1][2]=0.52444768167155973161042570784064;
583:     binterpt[2][2]=-2.3015205996945452158771370439586;
584:     binterpt[3][2]=0.76334325453713832352363565300308;

586:     TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
587:   }
588:   TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);
589:   TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);
590:   TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);
591:   TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);
592:   return(0);
593: }

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

598:    Not Collective

600:    Level: advanced

602: .seealso: TSRosWRegister(), TSRosWRegisterAll()
603: @*/
604: PetscErrorCode TSRosWRegisterDestroy(void)
605: {
606:   PetscErrorCode  ierr;
607:   RosWTableauLink link;

610:   while ((link = RosWTableauList)) {
611:     RosWTableau t = &link->tab;
612:     RosWTableauList = link->next;
613:     PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);
614:     PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);
615:     PetscFree2(t->bembed,t->bembedt);
616:     PetscFree(t->binterpt);
617:     PetscFree(t->name);
618:     PetscFree(link);
619:   }
620:   TSRosWRegisterAllCalled = PETSC_FALSE;
621:   return(0);
622: }

624: /*@C
625:   TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called
626:   from TSInitializePackage().

628:   Level: developer

630: .seealso: PetscInitialize()
631: @*/
632: PetscErrorCode TSRosWInitializePackage(void)
633: {

637:   if (TSRosWPackageInitialized) return(0);
638:   TSRosWPackageInitialized = PETSC_TRUE;
639:   TSRosWRegisterAll();
640:   PetscRegisterFinalize(TSRosWFinalizePackage);
641:   return(0);
642: }

644: /*@C
645:   TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is
646:   called from PetscFinalize().

648:   Level: developer

650: .seealso: PetscFinalize()
651: @*/
652: PetscErrorCode TSRosWFinalizePackage(void)
653: {

657:   TSRosWPackageInitialized = PETSC_FALSE;
658:   TSRosWRegisterDestroy();
659:   return(0);
660: }

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

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

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

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

681:    Level: advanced

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


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

736:   switch (s) {
737:   case 1: GammaInv[0] = 1./GammaInv[0]; break;
738:   case 2: PetscKernel_A_gets_inverse_A_2(GammaInv,0,PETSC_FALSE,NULL); break;
739:   case 3: PetscKernel_A_gets_inverse_A_3(GammaInv,0,PETSC_FALSE,NULL); break;
740:   case 4: PetscKernel_A_gets_inverse_A_4(GammaInv,0,PETSC_FALSE,NULL); break;
741:   case 5: {
742:     PetscInt  ipvt5[5];
743:     MatScalar work5[5*5];
744:     PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0,PETSC_FALSE,NULL); break;
745:   }
746:   case 6: PetscKernel_A_gets_inverse_A_6(GammaInv,0,PETSC_FALSE,NULL); break;
747:   case 7: PetscKernel_A_gets_inverse_A_7(GammaInv,0,PETSC_FALSE,NULL); break;
748:   default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s);
749:   }
750:   for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
751:   PetscFree(GammaInv);

753:   for (i=0; i<s; i++) {
754:     for (k=0; k<i+1; k++) {
755:       t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]);
756:       for (j=k+1; j<i+1; j++) {
757:         t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]);
758:       }
759:     }
760:   }

762:   for (i=0; i<s; i++) {
763:     for (j=0; j<s; j++) {
764:       t->At[i*s+j] = 0;
765:       for (k=0; k<s; k++) {
766:         t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j];
767:       }
768:     }
769:     t->bt[i] = 0;
770:     for (j=0; j<s; j++) {
771:       t->bt[i] += t->b[j] * t->GammaInv[j*s+i];
772:     }
773:     if (bembed) {
774:       t->bembedt[i] = 0;
775:       for (j=0; j<s; j++) {
776:         t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i];
777:       }
778:     }
779:   }
780:   t->ccfl = 1.0;                /* Fix this */

782:   t->pinterp = pinterp;
783:   PetscMalloc1(s*pinterp,&t->binterpt);
784:   PetscArraycpy(t->binterpt,binterpt,s*pinterp);
785:   link->next = RosWTableauList;
786:   RosWTableauList = link;
787:   return(0);
788: }

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

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

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

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

809:    Level: developer

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

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

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

845:   /* Step 3 */
846:   beta43       = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */
847:   beta32beta2p =  p44 / (b4*beta43);                    /* 7.15h */
848:   beta4jbetajp = (p32 - b3*beta32beta2p) / b4;
849:   M[0][0]      = b2;                                    M[0][1] = b3;                 M[0][2] = b4;
850:   M[1][0]      = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p;
851:   M[2][0]      = b4*beta43*a3*a3-p43;                   M[2][1] = -b4*beta43*a2*a2;   M[2][2] = 0;
852:   rhs[0]       = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32;
853:   PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);
854:   beta2p       = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]);
855:   beta3p       = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]);
856:   beta4p       = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]);

858:   /* Step 4: back-substitute */
859:   beta32 = beta32beta2p / beta2p;
860:   beta42 = (beta4jbetajp - beta43*beta3p) / beta2p;

862:   /* Step 5: 7.15f and 7.20, then 7.16 */
863:   a43 = 0;
864:   a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p);
865:   a42 = a32;

867:   A[0][0]     = 0;          A[0][1] = 0;   A[0][2] = 0;   A[0][3] = 0;
868:   A[1][0]     = a2;         A[1][1] = 0;   A[1][2] = 0;   A[1][3] = 0;
869:   A[2][0]     = a3-a32;     A[2][1] = a32; A[2][2] = 0;   A[2][3] = 0;
870:   A[3][0]     = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0;
871:   Gamma[0][0] = gamma;                        Gamma[0][1] = 0;              Gamma[0][2] = 0;              Gamma[0][3] = 0;
872:   Gamma[1][0] = beta2p-A[1][0];               Gamma[1][1] = gamma;          Gamma[1][2] = 0;              Gamma[1][3] = 0;
873:   Gamma[2][0] = beta3p-beta32-A[2][0];        Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma;          Gamma[2][3] = 0;
874:   Gamma[3][0] = beta4p-beta42-beta43-A[3][0]; Gamma[3][1] = beta42-A[3][1]; Gamma[3][2] = beta43-A[3][2]; Gamma[3][3] = gamma;
875:   b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4;

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

883:   {
884:     const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three;
885:     if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method");
886:   }
887:   TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);
888:   return(0);
889: }

891: /*
892:  The step completion formula is

894:  x1 = x0 + b^T Y

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

899:  x1e = x0 + be^T Y
900:      = x1 - b^T Y + be^T Y
901:      = x1 + (be - b)^T Y

903:  so we can evaluate the method of different order even after the step has been optimistically completed.
904: */
905: static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done)
906: {
907:   TS_RosW        *ros = (TS_RosW*)ts->data;
908:   RosWTableau    tab  = ros->tableau;
909:   PetscScalar    *w   = ros->work;
910:   PetscInt       i;

914:   if (order == tab->order) {
915:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
916:       VecCopy(ts->vec_sol,U);
917:       for (i=0; i<tab->s; i++) w[i] = tab->bt[i];
918:       VecMAXPY(U,tab->s,w,ros->Y);
919:     } else {VecCopy(ts->vec_sol,U);}
920:     if (done) *done = PETSC_TRUE;
921:     return(0);
922:   } else if (order == tab->order-1) {
923:     if (!tab->bembedt) goto unavailable;
924:     if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
925:       VecCopy(ts->vec_sol,U);
926:       for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i];
927:       VecMAXPY(U,tab->s,w,ros->Y);
928:     } else {                    /* Use rollback-and-recomplete formula (bembedt - bt) */
929:       for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
930:       VecCopy(ts->vec_sol,U);
931:       VecMAXPY(U,tab->s,w,ros->Y);
932:     }
933:     if (done) *done = PETSC_TRUE;
934:     return(0);
935:   }
936:   unavailable:
937:   if (done) *done = PETSC_FALSE;
938:   else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order);
939:   return(0);
940: }

942: static PetscErrorCode TSRollBack_RosW(TS ts)
943: {
944:   TS_RosW        *ros = (TS_RosW*)ts->data;

948:   VecCopy(ros->vec_sol_prev,ts->vec_sol);
949:   return(0);
950: }

952: static PetscErrorCode TSStep_RosW(TS ts)
953: {
954:   TS_RosW         *ros = (TS_RosW*)ts->data;
955:   RosWTableau     tab  = ros->tableau;
956:   const PetscInt  s    = tab->s;
957:   const PetscReal *At  = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv;
958:   const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
959:   const PetscBool *GammaZeroDiag = tab->GammaZeroDiag;
960:   PetscScalar     *w   = ros->work;
961:   Vec             *Y   = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage;
962:   SNES            snes;
963:   TSAdapt         adapt;
964:   PetscInt        i,j,its,lits;
965:   PetscInt        rejections = 0;
966:   PetscBool       stageok,accept = PETSC_TRUE;
967:   PetscReal       next_time_step = ts->time_step;
968:   PetscErrorCode  ierr;
969:   PetscInt        lag;

972:   if (!ts->steprollback) {
973:     VecCopy(ts->vec_sol,ros->vec_sol_prev);
974:   }

976:   ros->status = TS_STEP_INCOMPLETE;
977:   while (!ts->reason && ros->status != TS_STEP_COMPLETE) {
978:     const PetscReal h = ts->time_step;
979:     for (i=0; i<s; i++) {
980:       ros->stage_time = ts->ptime + h*ASum[i];
981:       TSPreStage(ts,ros->stage_time);
982:       if (GammaZeroDiag[i]) {
983:         ros->stage_explicit = PETSC_TRUE;
984:         ros->scoeff         = 1.;
985:       } else {
986:         ros->stage_explicit = PETSC_FALSE;
987:         ros->scoeff         = 1./Gamma[i*s+i];
988:       }

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

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

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

1001:       if (!ros->stage_explicit) {
1002:         TSGetSNES(ts,&snes);
1003:         if (!ros->recompute_jacobian && !i) {
1004:           SNESGetLagJacobian(snes,&lag);
1005:           if (lag == 1) {  /* use did not set a nontrival lag, so lag over all stages */
1006:             SNESSetLagJacobian(snes,-2); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */
1007:           }
1008:         }
1009:         SNESSolve(snes,NULL,Y[i]);
1010:         if (!ros->recompute_jacobian && i == s-1 && lag == 1) {
1011:           SNESSetLagJacobian(snes,lag); /* Set lag back to 1 so we know user did not set it */
1012:         }
1013:         SNESGetIterationNumber(snes,&its);
1014:         SNESGetLinearSolveIterations(snes,&lits);
1015:         ts->snes_its += its; ts->ksp_its += lits;
1016:       } else {
1017:         Mat J,Jp;
1018:         VecZeroEntries(Ydot); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
1019:         TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);
1020:         VecScale(Y[i],-1.0);
1021:         VecAXPY(Y[i],-1.0,Zdot); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/

1023:         VecZeroEntries(Zstage); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
1024:         for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j];
1025:         VecMAXPY(Zstage,i,w,Y);

1027:         /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
1028:         TSGetIJacobian(ts,&J,&Jp,NULL,NULL);
1029:         TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,J,Jp,PETSC_FALSE);
1030:         MatMult(J,Zstage,Zdot);
1031:         VecAXPY(Y[i],-1.0,Zdot);
1032:         ts->ksp_its += 1;

1034:         VecScale(Y[i],h);
1035:       }
1036:       TSPostStage(ts,ros->stage_time,i,Y);
1037:       TSGetAdapt(ts,&adapt);
1038:       TSAdaptCheckStage(adapt,ts,ros->stage_time,Y[i],&stageok);
1039:       if (!stageok) goto reject_step;
1040:     }

1042:     ros->status = TS_STEP_INCOMPLETE;
1043:     TSEvaluateStep_RosW(ts,tab->order,ts->vec_sol,NULL);
1044:     ros->status = TS_STEP_PENDING;
1045:     TSGetAdapt(ts,&adapt);
1046:     TSAdaptCandidatesClear(adapt);
1047:     TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
1048:     TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
1049:     ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1050:     if (!accept) { /* Roll back the current step */
1051:       TSRollBack_RosW(ts);
1052:       ts->time_step = next_time_step;
1053:       goto reject_step;
1054:     }

1056:     ts->ptime += ts->time_step;
1057:     ts->time_step = next_time_step;
1058:     break;

1060:   reject_step:
1061:     ts->reject++; accept = PETSC_FALSE;
1062:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1063:       ts->reason = TS_DIVERGED_STEP_REJECTED;
1064:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
1065:     }
1066:   }
1067:   return(0);
1068: }

1070: static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U)
1071: {
1072:   TS_RosW         *ros = (TS_RosW*)ts->data;
1073:   PetscInt        s    = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j;
1074:   PetscReal       h;
1075:   PetscReal       tt,t;
1076:   PetscScalar     *bt;
1077:   const PetscReal *Bt = ros->tableau->binterpt;
1078:   PetscErrorCode  ierr;
1079:   const PetscReal *GammaInv = ros->tableau->GammaInv;
1080:   PetscScalar     *w        = ros->work;
1081:   Vec             *Y        = ros->Y;

1084:   if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name);

1086:   switch (ros->status) {
1087:   case TS_STEP_INCOMPLETE:
1088:   case TS_STEP_PENDING:
1089:     h = ts->time_step;
1090:     t = (itime - ts->ptime)/h;
1091:     break;
1092:   case TS_STEP_COMPLETE:
1093:     h = ts->ptime - ts->ptime_prev;
1094:     t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
1095:     break;
1096:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
1097:   }
1098:   PetscMalloc1(s,&bt);
1099:   for (i=0; i<s; i++) bt[i] = 0;
1100:   for (j=0,tt=t; j<pinterp; j++,tt*=t) {
1101:     for (i=0; i<s; i++) {
1102:       bt[i] += Bt[i*pinterp+j] * tt;
1103:     }
1104:   }

1106:   /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
1107:   /* U <- 0*/
1108:   VecZeroEntries(U);
1109:   /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
1110:   for (j=0; j<s; j++) w[j] = 0;
1111:   for (j=0; j<s; j++) {
1112:     for (i=j; i<s; i++) {
1113:       w[j] +=  bt[i]*GammaInv[i*s+j];
1114:     }
1115:   }
1116:   VecMAXPY(U,i,w,Y);
1117:   /* U <- y(t) + U */
1118:   VecAXPY(U,1,ros->vec_sol_prev);

1120:   PetscFree(bt);
1121:   return(0);
1122: }

1124: /*------------------------------------------------------------*/

1126: static PetscErrorCode TSRosWTableauReset(TS ts)
1127: {
1128:   TS_RosW        *ros = (TS_RosW*)ts->data;
1129:   RosWTableau    tab  = ros->tableau;

1133:   if (!tab) return(0);
1134:   VecDestroyVecs(tab->s,&ros->Y);
1135:   PetscFree(ros->work);
1136:   return(0);
1137: }

1139: static PetscErrorCode TSReset_RosW(TS ts)
1140: {
1141:   TS_RosW        *ros = (TS_RosW*)ts->data;

1145:   TSRosWTableauReset(ts);
1146:   VecDestroy(&ros->Ydot);
1147:   VecDestroy(&ros->Ystage);
1148:   VecDestroy(&ros->Zdot);
1149:   VecDestroy(&ros->Zstage);
1150:   VecDestroy(&ros->vec_sol_prev);
1151:   return(0);
1152: }

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

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

1183: static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage)
1184: {

1188:   if (Ydot) {
1189:     if (dm && dm != ts->dm) {
1190:       DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);
1191:     }
1192:   }
1193:   if (Zdot) {
1194:     if (dm && dm != ts->dm) {
1195:       DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);
1196:     }
1197:   }
1198:   if (Ystage) {
1199:     if (dm && dm != ts->dm) {
1200:       DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);
1201:     }
1202:   }
1203:   if (Zstage) {
1204:     if (dm && dm != ts->dm) {
1205:       DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);
1206:     }
1207:   }
1208:   return(0);
1209: }

1211: static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx)
1212: {
1214:   return(0);
1215: }

1217: static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
1218: {
1219:   TS             ts = (TS)ctx;
1221:   Vec            Ydot,Zdot,Ystage,Zstage;
1222:   Vec            Ydotc,Zdotc,Ystagec,Zstagec;

1225:   TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);
1226:   TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);
1227:   MatRestrict(restrct,Ydot,Ydotc);
1228:   VecPointwiseMult(Ydotc,rscale,Ydotc);
1229:   MatRestrict(restrct,Ystage,Ystagec);
1230:   VecPointwiseMult(Ystagec,rscale,Ystagec);
1231:   MatRestrict(restrct,Zdot,Zdotc);
1232:   VecPointwiseMult(Zdotc,rscale,Zdotc);
1233:   MatRestrict(restrct,Zstage,Zstagec);
1234:   VecPointwiseMult(Zstagec,rscale,Zstagec);
1235:   TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);
1236:   TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);
1237:   return(0);
1238: }

1240: static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx)
1241: {
1243:   return(0);
1244: }

1246: static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
1247: {
1248:   TS             ts = (TS)ctx;
1250:   Vec            Ydot,Zdot,Ystage,Zstage;
1251:   Vec            Ydots,Zdots,Ystages,Zstages;

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

1257:   VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);
1258:   VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);

1260:   VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);
1261:   VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);

1263:   VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);
1264:   VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);

1266:   VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);
1267:   VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);

1269:   TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);
1270:   TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);
1271:   return(0);
1272: }

1274: /*
1275:   This defines the nonlinear equation that is to be solved with SNES
1276:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
1277: */
1278: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts)
1279: {
1280:   TS_RosW        *ros = (TS_RosW*)ts->data;
1282:   Vec            Ydot,Zdot,Ystage,Zstage;
1283:   PetscReal      shift = ros->scoeff / ts->time_step;
1284:   DM             dm,dmsave;

1287:   SNESGetDM(snes,&dm);
1288:   TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);
1289:   VecWAXPY(Ydot,shift,U,Zdot);    /* Ydot = shift*U + Zdot */
1290:   VecWAXPY(Ystage,1.0,U,Zstage);  /* Ystage = U + Zstage */
1291:   dmsave = ts->dm;
1292:   ts->dm = dm;
1293:   TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);
1294:   ts->dm = dmsave;
1295:   TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);
1296:   return(0);
1297: }

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

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

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

1326:   VecDuplicateVecs(ts->vec_sol,tab->s,&ros->Y);
1327:   PetscMalloc1(tab->s,&ros->work);
1328:   return(0);
1329: }

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

1340:   TSRosWTableauSetUp(ts);
1341:   VecDuplicate(ts->vec_sol,&ros->Ydot);
1342:   VecDuplicate(ts->vec_sol,&ros->Ystage);
1343:   VecDuplicate(ts->vec_sol,&ros->Zdot);
1344:   VecDuplicate(ts->vec_sol,&ros->Zstage);
1345:   VecDuplicate(ts->vec_sol,&ros->vec_sol_prev);
1346:   TSGetDM(ts,&dm);
1347:   DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);
1348:   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:   TSGetSNES(ts,&snes);
1351:   if (!((PetscObject)snes)->type_name) {
1352:     SNESSetType(snes,SNESKSPONLY);
1353:   }
1354:   DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);
1355:   if (rhsjacobian == TSComputeRHSJacobianConstant) {
1356:     Mat Amat,Pmat;

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

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

1387:   PetscOptionsHead(PetscOptionsObject,"RosW ODE solver options");
1388:   {
1389:     RosWTableauLink link;
1390:     PetscInt        count,choice;
1391:     PetscBool       flg;
1392:     const char      **namelist;

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

1401:     PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);
1402:   }
1403:   PetscOptionsTail();
1404:   /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1405:   TSGetSNES(ts,&snes);
1406:   if (!((PetscObject)snes)->type_name) {
1407:     SNESSetType(snes,SNESKSPONLY);
1408:   }
1409:   return(0);
1410: }

1412: static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer)
1413: {
1414:   TS_RosW        *ros = (TS_RosW*)ts->data;
1415:   PetscBool      iascii;

1419:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1420:   if (iascii) {
1421:     RosWTableau tab  = ros->tableau;
1422:     TSRosWType  rostype;
1423:     char        buf[512];
1424:     PetscInt    i;
1425:     PetscReal   abscissa[512];
1426:     TSRosWGetType(ts,&rostype);
1427:     PetscViewerASCIIPrintf(viewer,"  Rosenbrock-W %s\n",rostype);
1428:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);
1429:     PetscViewerASCIIPrintf(viewer,"  Abscissa of A       = %s\n",buf);
1430:     for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i];
1431:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);
1432:     PetscViewerASCIIPrintf(viewer,"  Abscissa of A+Gamma = %s\n",buf);
1433:   }
1434:   return(0);
1435: }

1437: static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer)
1438: {
1440:   SNES           snes;
1441:   TSAdapt        adapt;

1444:   TSGetAdapt(ts,&adapt);
1445:   TSAdaptLoad(adapt,viewer);
1446:   TSGetSNES(ts,&snes);
1447:   SNESLoad(snes,viewer);
1448:   /* function and Jacobian context for SNES when used with TS is always ts object */
1449:   SNESSetFunction(snes,NULL,NULL,ts);
1450:   SNESSetJacobian(snes,NULL,NULL,NULL,ts);
1451:   return(0);
1452: }

1454: /*@C
1455:   TSRosWSetType - Set the type of Rosenbrock-W scheme

1457:   Logically collective

1459:   Input Parameters:
1460: +  ts - timestepping context
1461: -  roswtype - type of Rosenbrock-W scheme

1463:   Level: beginner

1465: .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3
1466: @*/
1467: PetscErrorCode TSRosWSetType(TS ts,TSRosWType roswtype)
1468: {

1474:   PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,roswtype));
1475:   return(0);
1476: }

1478: /*@C
1479:   TSRosWGetType - Get the type of Rosenbrock-W scheme

1481:   Logically collective

1483:   Input Parameter:
1484: .  ts - timestepping context

1486:   Output Parameter:
1487: .  rostype - type of Rosenbrock-W scheme

1489:   Level: intermediate

1491: .seealso: TSRosWGetType()
1492: @*/
1493: PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype)
1494: {

1499:   PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));
1500:   return(0);
1501: }

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

1506:   Logically collective

1508:   Input Parameters:
1509: +  ts - timestepping context
1510: -  flg - PETSC_TRUE to recompute the Jacobian at each stage

1512:   Level: intermediate

1514: .seealso: TSRosWGetType()
1515: @*/
1516: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg)
1517: {

1522:   PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));
1523:   return(0);
1524: }

1526: static PetscErrorCode  TSRosWGetType_RosW(TS ts,TSRosWType *rostype)
1527: {
1528:   TS_RosW        *ros = (TS_RosW*)ts->data;

1531:   *rostype = ros->tableau->name;
1532:   return(0);
1533: }

1535: static PetscErrorCode  TSRosWSetType_RosW(TS ts,TSRosWType rostype)
1536: {
1537:   TS_RosW         *ros = (TS_RosW*)ts->data;
1538:   PetscErrorCode  ierr;
1539:   PetscBool       match;
1540:   RosWTableauLink link;

1543:   if (ros->tableau) {
1544:     PetscStrcmp(ros->tableau->name,rostype,&match);
1545:     if (match) return(0);
1546:   }
1547:   for (link = RosWTableauList; link; link=link->next) {
1548:     PetscStrcmp(link->tab.name,rostype,&match);
1549:     if (match) {
1550:       if (ts->setupcalled) {TSRosWTableauReset(ts);}
1551:       ros->tableau = &link->tab;
1552:       if (ts->setupcalled) {TSRosWTableauSetUp(ts);}
1553:       ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1554:       return(0);
1555:     }
1556:   }
1557:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype);
1558: }

1560: static PetscErrorCode  TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg)
1561: {
1562:   TS_RosW *ros = (TS_RosW*)ts->data;

1565:   ros->recompute_jacobian = flg;
1566:   return(0);
1567: }

1569: static PetscErrorCode TSDestroy_RosW(TS ts)
1570: {

1574:   TSReset_RosW(ts);
1575:   if (ts->dm) {
1576:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);
1577:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);
1578:   }
1579:   PetscFree(ts->data);
1580:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);
1581:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);
1582:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);
1583:   return(0);
1584: }

1586: /* ------------------------------------------------------------ */
1587: /*MC
1588:       TSROSW - ODE solver using Rosenbrock-W schemes

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

1594:   Notes:
1595:   This method currently only works with autonomous ODE and DAE.

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

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

1601:   Developer Notes:
1602:   Rosenbrock-W methods are typically specified for autonomous ODE

1604: $  udot = f(u)

1606:   by the stage equations

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

1610:   and step completion formula

1612: $  u_1 = u_0 + sum_j b_j k_j

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

1618: $  y_i = gamma_ij k_j

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

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

1624:   to rewrite the method as

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

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

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

1633:    or, more compactly in tensor notation

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

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

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

1643:    with initial guess y_i = 0.

1645:   Level: beginner

1647: .seealso:  TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSWTHETA1, TSROSWTHETA2, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3,
1648:            TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L
1649: M*/
1650: PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1651: {
1652:   TS_RosW        *ros;

1656:   TSRosWInitializePackage();

1658:   ts->ops->reset          = TSReset_RosW;
1659:   ts->ops->destroy        = TSDestroy_RosW;
1660:   ts->ops->view           = TSView_RosW;
1661:   ts->ops->load           = TSLoad_RosW;
1662:   ts->ops->setup          = TSSetUp_RosW;
1663:   ts->ops->step           = TSStep_RosW;
1664:   ts->ops->interpolate    = TSInterpolate_RosW;
1665:   ts->ops->evaluatestep   = TSEvaluateStep_RosW;
1666:   ts->ops->rollback       = TSRollBack_RosW;
1667:   ts->ops->setfromoptions = TSSetFromOptions_RosW;
1668:   ts->ops->snesfunction   = SNESTSFormFunction_RosW;
1669:   ts->ops->snesjacobian   = SNESTSFormJacobian_RosW;

1671:   ts->usessnes = PETSC_TRUE;

1673:   PetscNewLog(ts,&ros);
1674:   ts->data = (void*)ros;

1676:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);
1677:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);
1678:   PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);

1680:   TSRosWSetType(ts,TSRosWDefault);
1681:   return(0);
1682: }