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