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: [](chapter_ts), `TSROSW`
74: M*/
76: /*MC
77: TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method).
79: Only an approximate Jacobian is needed.
81: Level: intermediate
83: .seealso: [](chapter_ts), `TSROSW`
84: M*/
86: /*MC
87: TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme.
89: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P.
91: Level: intermediate
93: .seealso: [](chapter_ts), `TSROSW`
94: M*/
96: /*MC
97: TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme.
99: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M.
101: Level: intermediate
103: .seealso: [](chapter_ts), `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: Level: intermediate
115: References:
116: . * - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005.
118: .seealso: [](chapter_ts), `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: Level: intermediate
130: References:
131: . * - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005.
133: .seealso: [](chapter_ts), `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: Level: intermediate
145: References:
146: . * - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.
148: .seealso: [](chapter_ts), `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: Level: intermediate
163: References:
164: . * - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.
166: .seealso: [](chapter_ts), `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: Level: intermediate
178: References:
179: . * - Emil Constantinescu
181: .seealso: [](chapter_ts), `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: Level: intermediate
193: References:
194: . * - Emil Constantinescu
196: .seealso: [](chapter_ts), `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: Level: intermediate
208: References:
209: . * - Emil Constantinescu
211: .seealso: [](chapter_ts), `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: Level: intermediate
225: References:
226: + * - Kaps and Rentrop, Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979.
227: - * - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.
229: Hairer's code ros4.f
231: .seealso: [](chapter_ts), `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: Level: intermediate
245: References:
246: + * - Shampine, Implementation of Rosenbrock methods, 1982.
247: - * - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.
249: Hairer's code ros4.f
251: .seealso: [](chapter_ts), `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: Level: intermediate
265: References:
266: + * - van Veldhuizen, D stability and Kaps Rentrop methods, 1984.
267: - * - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.
269: Hairer's code ros4.f
271: .seealso: [](chapter_ts), `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: Level: intermediate
285: References:
286: . * - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2.
288: Hairer's code ros4.f
290: .seealso: [](chapter_ts), `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: [](chapter_ts), `TSROSW`, `TSRosWRegisterDestroy()`
301: @*/
302: PetscErrorCode TSRosWRegisterAll(void)
303: {
304: if (TSRosWRegisterAllCalled) return 0;
305: TSRosWRegisterAllCalled = PETSC_TRUE;
307: {
308: const PetscReal A = 0;
309: const PetscReal Gamma = 1;
310: const PetscReal b = 1;
311: const PetscReal binterpt = 1;
313: TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt);
314: }
316: {
317: const PetscReal A = 0;
318: const PetscReal Gamma = 0.5;
319: const PetscReal b = 1;
320: const PetscReal binterpt = 1;
322: TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt);
323: }
325: {
326: /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */
327: const PetscReal A[2][2] =
328: {
329: {0, 0},
330: {1., 0}
331: },
332: Gamma[2][2] = {{1.707106781186547524401, 0}, {-2. * 1.707106781186547524401, 1.707106781186547524401}}, b[2] = {0.5, 0.5}, b1[2] = {1.0, 0.0};
333: PetscReal binterpt[2][2];
334: binterpt[0][0] = 1.707106781186547524401 - 1.0;
335: binterpt[1][0] = 2.0 - 1.707106781186547524401;
336: binterpt[0][1] = 1.707106781186547524401 - 1.5;
337: binterpt[1][1] = 1.5 - 1.707106781186547524401;
339: TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]);
340: }
341: {
342: /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */
343: const PetscReal A[2][2] =
344: {
345: {0, 0},
346: {1., 0}
347: },
348: Gamma[2][2] = {{0.2928932188134524755992, 0}, {-2. * 0.2928932188134524755992, 0.2928932188134524755992}}, b[2] = {0.5, 0.5}, b1[2] = {1.0, 0.0};
349: PetscReal binterpt[2][2];
350: binterpt[0][0] = 0.2928932188134524755992 - 1.0;
351: binterpt[1][0] = 2.0 - 0.2928932188134524755992;
352: binterpt[0][1] = 0.2928932188134524755992 - 1.5;
353: binterpt[1][1] = 1.5 - 0.2928932188134524755992;
355: TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]);
356: }
357: {
358: /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */
359: PetscReal binterpt[3][2];
360: const PetscReal A[3][3] =
361: {
362: {0, 0, 0},
363: {1.5773502691896257e+00, 0, 0},
364: {0.5, 0, 0}
365: },
366: Gamma[3][3] = {{7.8867513459481287e-01, 0, 0}, {-1.5773502691896257e+00, 7.8867513459481287e-01, 0}, {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01}}, b[3] = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01}, b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01};
368: binterpt[0][0] = -0.8094010767585034;
369: binterpt[1][0] = -0.5;
370: binterpt[2][0] = 2.3094010767585034;
371: binterpt[0][1] = 0.9641016151377548;
372: binterpt[1][1] = 0.5;
373: binterpt[2][1] = -1.4641016151377548;
375: TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]);
376: }
377: {
378: PetscReal binterpt[4][3];
379: /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
380: const PetscReal A[4][4] =
381: {
382: {0, 0, 0, 0},
383: {8.7173304301691801e-01, 0, 0, 0},
384: {8.4457060015369423e-01, -1.1299064236484185e-01, 0, 0},
385: {0, 0, 1., 0}
386: },
387: Gamma[4][4] = {{4.3586652150845900e-01, 0, 0, 0}, {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0}, {-9.0338057013044082e-01, 5.4180672388095326e-02, 4.3586652150845900e-01, 0}, {2.4212380706095346e-01, -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01}}, b[4] = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01}, b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01};
389: binterpt[0][0] = 1.0564298455794094;
390: binterpt[1][0] = 2.296429974281067;
391: binterpt[2][0] = -1.307599564525376;
392: binterpt[3][0] = -1.045260255335102;
393: binterpt[0][1] = -1.3864882699759573;
394: binterpt[1][1] = -8.262611700275677;
395: binterpt[2][1] = 7.250979895056055;
396: binterpt[3][1] = 2.398120075195581;
397: binterpt[0][2] = 0.5721822314575016;
398: binterpt[1][2] = 4.742931142090097;
399: binterpt[2][2] = -4.398120075195578;
400: binterpt[3][2] = -0.9169932983520199;
402: TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
403: }
404: {
405: /* const PetscReal g = 0.5; Directly written in-place below */
406: const PetscReal A[4][4] =
407: {
408: {0, 0, 0, 0},
409: {0, 0, 0, 0},
410: {1., 0, 0, 0},
411: {0.75, -0.25, 0.5, 0}
412: },
413: Gamma[4][4] = {{0.5, 0, 0, 0}, {1., 0.5, 0, 0}, {-0.25, -0.25, 0.5, 0}, {1. / 12, 1. / 12, -2. / 3, 0.5}}, b[4] = {5. / 6, -1. / 6, -1. / 6, 0.5}, b2[4] = {0.75, -0.25, 0.5, 0};
415: TSRosWRegister(TSROSWRODAS3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL);
416: }
417: {
418: /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */
419: const PetscReal A[3][3] =
420: {
421: {0, 0, 0},
422: {0.43586652150845899941601945119356, 0, 0},
423: {0.43586652150845899941601945119356, 0, 0}
424: },
425: Gamma[3][3] = {{0.43586652150845899941601945119356, 0, 0}, {-0.19294655696029095575009695436041, 0.43586652150845899941601945119356, 0}, {0, 1.74927148125794685173529749738960, 0.43586652150845899941601945119356}}, b[3] = {-0.75457412385404315829818998646589, 1.94100407061964420292840123379419, -0.18642994676560104463021124732829}, b2[3] = {-1.53358745784149585370766523913002, 2.81745131148625772213931745457622, -0.28386385364476186843165221544619};
427: PetscReal binterpt[3][2];
428: binterpt[0][0] = 3.793692883777660870425141387941;
429: binterpt[1][0] = -2.918692883777660870425141387941;
430: binterpt[2][0] = 0.125;
431: binterpt[0][1] = -0.725741064379812106687651020584;
432: binterpt[1][1] = 0.559074397713145440020984353917;
433: binterpt[2][1] = 0.16666666666666666666666666666667;
435: TSRosWRegister(TSROSWSANDU3, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]);
436: }
437: {
438: /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0;
439: * Direct evaluation: s3 = 1.732050807568877293527;
440: * g = 0.7886751345948128822546;
441: * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */
442: const PetscReal A[3][3] =
443: {
444: {0, 0, 0},
445: {1, 0, 0},
446: {0.25, 0.25, 0}
447: },
448: Gamma[3][3] = {{0, 0, 0}, {(-3.0 - 1.732050807568877293527) / 6.0, 0.7886751345948128822546, 0}, {(-3.0 - 1.732050807568877293527) / 24.0, (-3.0 - 1.732050807568877293527) / 8.0, 0.7886751345948128822546}}, b[3] = {1. / 6., 1. / 6., 2. / 3.}, b2[3] = {1. / 4., 1. / 4., 1. / 2.};
449: PetscReal binterpt[3][2];
451: binterpt[0][0] = 0.089316397477040902157517886164709;
452: binterpt[1][0] = -0.91068360252295909784248211383529;
453: binterpt[2][0] = 1.8213672050459181956849642276706;
454: binterpt[0][1] = 0.077350269189625764509148780501957;
455: binterpt[1][1] = 1.077350269189625764509148780502;
456: binterpt[2][1] = -1.1547005383792515290182975610039;
458: TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]);
459: }
461: {
462: const PetscReal A[4][4] =
463: {
464: {0, 0, 0, 0},
465: {1. / 2., 0, 0, 0},
466: {1. / 2., 1. / 2., 0, 0},
467: {1. / 6., 1. / 6., 1. / 6., 0}
468: },
469: Gamma[4][4] = {{1. / 2., 0, 0, 0}, {0.0, 1. / 4., 0, 0}, {-2., -2. / 3., 2. / 3., 0}, {1. / 2., 5. / 36., -2. / 9, 0}}, b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}, b2[4] = {1. / 8., 3. / 4., 1. / 8., 0};
470: PetscReal binterpt[4][3];
472: binterpt[0][0] = 6.25;
473: binterpt[1][0] = -30.25;
474: binterpt[2][0] = 1.75;
475: binterpt[3][0] = 23.25;
476: binterpt[0][1] = -9.75;
477: binterpt[1][1] = 58.75;
478: binterpt[2][1] = -3.25;
479: binterpt[3][1] = -45.75;
480: binterpt[0][2] = 3.6666666666666666666666666666667;
481: binterpt[1][2] = -28.333333333333333333333333333333;
482: binterpt[2][2] = 1.6666666666666666666666666666667;
483: binterpt[3][2] = 23.;
485: TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
486: }
488: {
489: const PetscReal A[4][4] =
490: {
491: {0, 0, 0, 0},
492: {1. / 2., 0, 0, 0},
493: {1. / 2., 1. / 2., 0, 0},
494: {1. / 6., 1. / 6., 1. / 6., 0}
495: },
496: Gamma[4][4] = {{1. / 2., 0, 0, 0}, {0.0, 3. / 4., 0, 0}, {-2. / 3., -23. / 9., 2. / 9., 0}, {1. / 18., 65. / 108., -2. / 27, 0}}, b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}, b2[4] = {3. / 16., 10. / 16., 3. / 16., 0};
497: PetscReal binterpt[4][3];
499: binterpt[0][0] = 1.6911764705882352941176470588235;
500: binterpt[1][0] = 3.6813725490196078431372549019608;
501: binterpt[2][0] = 0.23039215686274509803921568627451;
502: binterpt[3][0] = -4.6029411764705882352941176470588;
503: binterpt[0][1] = -0.95588235294117647058823529411765;
504: binterpt[1][1] = -6.2401960784313725490196078431373;
505: binterpt[2][1] = -0.31862745098039215686274509803922;
506: binterpt[3][1] = 7.5147058823529411764705882352941;
507: binterpt[0][2] = -0.56862745098039215686274509803922;
508: binterpt[1][2] = 2.7254901960784313725490196078431;
509: binterpt[2][2] = 0.25490196078431372549019607843137;
510: binterpt[3][2] = -2.4117647058823529411764705882353;
512: TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
513: }
515: {
516: PetscReal A[4][4], Gamma[4][4], b[4], b2[4];
517: PetscReal binterpt[4][3];
519: Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816;
520: Gamma[0][1] = 0;
521: Gamma[0][2] = 0;
522: Gamma[0][3] = 0;
523: Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476;
524: Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816;
525: Gamma[1][2] = 0;
526: Gamma[1][3] = 0;
527: Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903;
528: Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131;
529: Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816;
530: Gamma[2][3] = 0;
531: Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783;
532: Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984;
533: Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198;
534: Gamma[3][3] = 0;
536: A[0][0] = 0;
537: A[0][1] = 0;
538: A[0][2] = 0;
539: A[0][3] = 0;
540: A[1][0] = 0.8717330430169179988320388950590125027645343373957631;
541: A[1][1] = 0;
542: A[1][2] = 0;
543: A[1][3] = 0;
544: A[2][0] = 0.5275890119763004115618079766722914408876108660811028;
545: A[2][1] = 0.07241098802369958843819203208518599088698057726988732;
546: A[2][2] = 0;
547: A[2][3] = 0;
548: A[3][0] = 0.3990960076760701320627260685975778145384666450351314;
549: A[3][1] = -0.4375576546135194437228463747348862825846903771419953;
550: A[3][2] = 1.038461646937449311660120300601880176655352737312713;
551: A[3][3] = 0;
553: b[0] = 0.1876410243467238251612921333138006734899663569186926;
554: b[1] = -0.5952974735769549480478230473706443582188442040780541;
555: b[2] = 0.9717899277217721234705114616271378792182450260943198;
556: b[3] = 0.4358665215084589994160194475295062513822671686978816;
558: b2[0] = 0.2147402862233891404862383521089097657790734483804460;
559: b2[1] = -0.4851622638849390928209050538171743017757490232519684;
560: b2[2] = 0.8687250025203875511662123688667549217531982787600080;
561: b2[3] = 0.4016969751411624011684543450940068201770721128357014;
563: binterpt[0][0] = 2.2565812720167954547104627844105;
564: binterpt[1][0] = 1.349166413351089573796243820819;
565: binterpt[2][0] = -2.4695174540533503758652847586647;
566: binterpt[3][0] = -0.13623023131453465264142184656474;
567: binterpt[0][1] = -3.0826699111559187902922463354557;
568: binterpt[1][1] = -2.4689115685996042534544925650515;
569: binterpt[2][1] = 5.7428279814696677152129332773553;
570: binterpt[3][1] = -0.19124650171414467146619437684812;
571: binterpt[0][2] = 1.0137296634858471607430756831148;
572: binterpt[1][2] = 0.52444768167155973161042570784064;
573: binterpt[2][2] = -2.3015205996945452158771370439586;
574: binterpt[3][2] = 0.76334325453713832352363565300308;
576: TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]);
577: }
578: TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_DEFAULT, PETSC_DEFAULT, 0, -0.1282612945269037e+01);
579: TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_DEFAULT, PETSC_DEFAULT, 0, 125. / 108.);
580: TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.355958941201148);
581: TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_DEFAULT, PETSC_DEFAULT, 0, -1.093502252409163);
582: return 0;
583: }
585: /*@C
586: TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`.
588: Not Collective
590: Level: advanced
592: .seealso: [](chapter_ts), `TSRosWRegister()`, `TSRosWRegisterAll()`
593: @*/
594: PetscErrorCode TSRosWRegisterDestroy(void)
595: {
596: RosWTableauLink link;
598: while ((link = RosWTableauList)) {
599: RosWTableau t = &link->tab;
600: RosWTableauList = link->next;
601: PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum);
602: PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr);
603: PetscFree2(t->bembed, t->bembedt);
604: PetscFree(t->binterpt);
605: PetscFree(t->name);
606: PetscFree(link);
607: }
608: TSRosWRegisterAllCalled = PETSC_FALSE;
609: return 0;
610: }
612: /*@C
613: TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called
614: from `TSInitializePackage()`.
616: Level: developer
618: .seealso: [](chapter_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()`
619: @*/
620: PetscErrorCode TSRosWInitializePackage(void)
621: {
622: if (TSRosWPackageInitialized) return 0;
623: TSRosWPackageInitialized = PETSC_TRUE;
624: TSRosWRegisterAll();
625: PetscRegisterFinalize(TSRosWFinalizePackage);
626: return 0;
627: }
629: /*@C
630: TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is
631: called from `PetscFinalize()`.
633: Level: developer
635: .seealso: [](chapter_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()`
636: @*/
637: PetscErrorCode TSRosWFinalizePackage(void)
638: {
639: TSRosWPackageInitialized = PETSC_FALSE;
640: TSRosWRegisterDestroy();
641: return 0;
642: }
644: /*@C
645: TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
647: Not Collective, but the same schemes should be registered on all processes on which they will be used
649: Input Parameters:
650: + name - identifier for method
651: . order - approximation order of method
652: . s - number of stages, this is the dimension of the matrices below
653: . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular
654: . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal
655: . b - Step completion table (dimension s)
656: . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available)
657: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt
658: - binterpt - Coefficients of the interpolation formula (dimension s*pinterp)
660: Level: advanced
662: Note:
663: Several Rosenbrock W methods are provided, this function is only needed to create new methods.
665: .seealso: [](chapter_ts), `TSROSW`
666: @*/
667: PetscErrorCode TSRosWRegister(TSRosWType name, PetscInt order, PetscInt s, const PetscReal A[], const PetscReal Gamma[], const PetscReal b[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[])
668: {
669: RosWTableauLink link;
670: RosWTableau t;
671: PetscInt i, j, k;
672: PetscScalar *GammaInv;
680: TSRosWInitializePackage();
681: PetscNew(&link);
682: t = &link->tab;
683: PetscStrallocpy(name, &t->name);
684: t->order = order;
685: t->s = s;
686: PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum);
687: PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr);
688: PetscArraycpy(t->A, A, s * s);
689: PetscArraycpy(t->Gamma, Gamma, s * s);
690: PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s);
691: PetscArraycpy(t->b, b, s);
692: if (bembed) {
693: PetscMalloc2(s, &t->bembed, s, &t->bembedt);
694: PetscArraycpy(t->bembed, bembed, s);
695: }
696: for (i = 0; i < s; i++) {
697: t->ASum[i] = 0;
698: t->GammaSum[i] = 0;
699: for (j = 0; j < s; j++) {
700: t->ASum[i] += A[i * s + j];
701: t->GammaSum[i] += Gamma[i * s + j];
702: }
703: }
704: PetscMalloc1(s * s, &GammaInv); /* Need to use Scalar for inverse, then convert back to Real */
705: for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i];
706: for (i = 0; i < s; i++) {
707: if (Gamma[i * s + i] == 0.0) {
708: GammaInv[i * s + i] = 1.0;
709: t->GammaZeroDiag[i] = PETSC_TRUE;
710: } else {
711: t->GammaZeroDiag[i] = PETSC_FALSE;
712: }
713: }
715: switch (s) {
716: case 1:
717: GammaInv[0] = 1. / GammaInv[0];
718: break;
719: case 2:
720: PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL);
721: break;
722: case 3:
723: PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL);
724: break;
725: case 4:
726: PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL);
727: break;
728: case 5: {
729: PetscInt ipvt5[5];
730: MatScalar work5[5 * 5];
731: PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL);
732: break;
733: }
734: case 6:
735: PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL);
736: break;
737: case 7:
738: PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL);
739: break;
740: default:
741: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s);
742: }
743: for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
744: PetscFree(GammaInv);
746: for (i = 0; i < s; i++) {
747: for (k = 0; k < i + 1; k++) {
748: t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]);
749: for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]);
750: }
751: }
753: for (i = 0; i < s; i++) {
754: for (j = 0; j < s; j++) {
755: t->At[i * s + j] = 0;
756: for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j];
757: }
758: t->bt[i] = 0;
759: for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i];
760: if (bembed) {
761: t->bembedt[i] = 0;
762: for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i];
763: }
764: }
765: t->ccfl = 1.0; /* Fix this */
767: t->pinterp = pinterp;
768: PetscMalloc1(s * pinterp, &t->binterpt);
769: PetscArraycpy(t->binterpt, binterpt, s * pinterp);
770: link->next = RosWTableauList;
771: RosWTableauList = link;
772: return 0;
773: }
775: /*@C
776: TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices
778: Not Collective, but the same schemes should be registered on all processes on which they will be used
780: Input Parameters:
781: + name - identifier for method
782: . gamma - leading coefficient (diagonal entry)
783: . a2 - design parameter, see Table 7.2 of Hairer&Wanner
784: . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22)
785: . b3 - design parameter, see Table 7.2 of Hairer&Wanner
786: . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner
787: - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer
789: Level: developer
791: Notes:
792: This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2.
793: It is used here to implement several methods from the book and can be used to experiment with new methods.
794: It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions.
796: .seealso: [](chapter_ts), `TSRosW`, `TSRosWRegister()`
797: @*/
798: PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4)
799: {
800: /* Declare numeric constants so they can be quad precision without being truncated at double */
801: const PetscReal one = 1, two = 2, three = 3, four = 4, five = 5, six = 6, eight = 8, twelve = 12, twenty = 20, twentyfour = 24, p32 = one / six - gamma + gamma * gamma, p42 = one / eight - gamma / three, p43 = one / twelve - gamma / three, p44 = one / twentyfour - gamma / two + three / two * gamma * gamma - gamma * gamma * gamma, p56 = one / twenty - gamma / four;
802: PetscReal a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp;
803: PetscReal A[4][4], Gamma[4][4], b[4], bm[4];
804: PetscScalar M[3][3], rhs[3];
806: /* Step 1: choose Gamma (input) */
807: /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */
808: if (a3 == PETSC_DEFAULT) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */
809: a4 = a3; /* consequence of 7.20 */
811: /* Solve order conditions 7.15a, 7.15c, 7.15e */
812: M[0][0] = one;
813: M[0][1] = one;
814: M[0][2] = one; /* 7.15a */
815: M[1][0] = 0.0;
816: M[1][1] = a2 * a2;
817: M[1][2] = a4 * a4; /* 7.15c */
818: M[2][0] = 0.0;
819: M[2][1] = a2 * a2 * a2;
820: M[2][2] = a4 * a4 * a4; /* 7.15e */
821: rhs[0] = one - b3;
822: rhs[1] = one / three - a3 * a3 * b3;
823: rhs[2] = one / four - a3 * a3 * a3 * b3;
824: PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL);
825: b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
826: b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
827: b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);
829: /* Step 3 */
830: beta43 = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */
831: beta32beta2p = p44 / (b4 * beta43); /* 7.15h */
832: beta4jbetajp = (p32 - b3 * beta32beta2p) / b4;
833: M[0][0] = b2;
834: M[0][1] = b3;
835: M[0][2] = b4;
836: M[1][0] = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp;
837: M[1][1] = a2 * a2 * beta4jbetajp;
838: M[1][2] = -a2 * a2 * beta32beta2p;
839: M[2][0] = b4 * beta43 * a3 * a3 - p43;
840: M[2][1] = -b4 * beta43 * a2 * a2;
841: M[2][2] = 0;
842: rhs[0] = one / two - gamma;
843: rhs[1] = 0;
844: rhs[2] = -a2 * a2 * p32;
845: PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL);
846: beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
847: beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
848: beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);
850: /* Step 4: back-substitute */
851: beta32 = beta32beta2p / beta2p;
852: beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p;
854: /* Step 5: 7.15f and 7.20, then 7.16 */
855: a43 = 0;
856: a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p);
857: a42 = a32;
859: A[0][0] = 0;
860: A[0][1] = 0;
861: A[0][2] = 0;
862: A[0][3] = 0;
863: A[1][0] = a2;
864: A[1][1] = 0;
865: A[1][2] = 0;
866: A[1][3] = 0;
867: A[2][0] = a3 - a32;
868: A[2][1] = a32;
869: A[2][2] = 0;
870: A[2][3] = 0;
871: A[3][0] = a4 - a43 - a42;
872: A[3][1] = a42;
873: A[3][2] = a43;
874: A[3][3] = 0;
875: Gamma[0][0] = gamma;
876: Gamma[0][1] = 0;
877: Gamma[0][2] = 0;
878: Gamma[0][3] = 0;
879: Gamma[1][0] = beta2p - A[1][0];
880: Gamma[1][1] = gamma;
881: Gamma[1][2] = 0;
882: Gamma[1][3] = 0;
883: Gamma[2][0] = beta3p - beta32 - A[2][0];
884: Gamma[2][1] = beta32 - A[2][1];
885: Gamma[2][2] = gamma;
886: Gamma[2][3] = 0;
887: Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0];
888: Gamma[3][1] = beta42 - A[3][1];
889: Gamma[3][2] = beta43 - A[3][2];
890: Gamma[3][3] = gamma;
891: b[0] = b1;
892: b[1] = b2;
893: b[2] = b3;
894: b[3] = b4;
896: /* Construct embedded formula using given e4. We are solving Equation 7.18. */
897: bm[3] = b[3] - e4 * gamma; /* using definition of E4 */
898: bm[2] = (p32 - beta4jbetajp * bm[3]) / (beta32 * beta2p); /* fourth row of 7.18 */
899: bm[1] = (one / two - gamma - beta3p * bm[2] - beta4p * bm[3]) / beta2p; /* second row */
900: bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */
902: {
903: const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three;
905: }
906: TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL);
907: return 0;
908: }
910: /*
911: The step completion formula is
913: x1 = x0 + b^T Y
915: where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been
916: updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write
918: x1e = x0 + be^T Y
919: = x1 - b^T Y + be^T Y
920: = x1 + (be - b)^T Y
922: so we can evaluate the method of different order even after the step has been optimistically completed.
923: */
924: static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done)
925: {
926: TS_RosW *ros = (TS_RosW *)ts->data;
927: RosWTableau tab = ros->tableau;
928: PetscScalar *w = ros->work;
929: PetscInt i;
931: if (order == tab->order) {
932: if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
933: VecCopy(ts->vec_sol, U);
934: for (i = 0; i < tab->s; i++) w[i] = tab->bt[i];
935: VecMAXPY(U, tab->s, w, ros->Y);
936: } else VecCopy(ts->vec_sol, U);
937: if (done) *done = PETSC_TRUE;
938: return 0;
939: } else if (order == tab->order - 1) {
940: if (!tab->bembedt) goto unavailable;
941: if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
942: VecCopy(ts->vec_sol, U);
943: for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i];
944: VecMAXPY(U, tab->s, w, ros->Y);
945: } else { /* Use rollback-and-recomplete formula (bembedt - bt) */
946: for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
947: VecCopy(ts->vec_sol, U);
948: VecMAXPY(U, tab->s, w, ros->Y);
949: }
950: if (done) *done = PETSC_TRUE;
951: return 0;
952: }
953: unavailable:
954: if (done) *done = PETSC_FALSE;
955: else
956: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Rosenbrock-W '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", tab->name,
957: tab->order, order);
958: return 0;
959: }
961: static PetscErrorCode TSRollBack_RosW(TS ts)
962: {
963: TS_RosW *ros = (TS_RosW *)ts->data;
965: VecCopy(ros->vec_sol_prev, ts->vec_sol);
966: return 0;
967: }
969: static PetscErrorCode TSStep_RosW(TS ts)
970: {
971: TS_RosW *ros = (TS_RosW *)ts->data;
972: RosWTableau tab = ros->tableau;
973: const PetscInt s = tab->s;
974: const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv;
975: const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
976: const PetscBool *GammaZeroDiag = tab->GammaZeroDiag;
977: PetscScalar *w = ros->work;
978: Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage;
979: SNES snes;
980: TSAdapt adapt;
981: PetscInt i, j, its, lits;
982: PetscInt rejections = 0;
983: PetscBool stageok, accept = PETSC_TRUE;
984: PetscReal next_time_step = ts->time_step;
985: PetscInt lag;
987: if (!ts->steprollback) VecCopy(ts->vec_sol, ros->vec_sol_prev);
989: ros->status = TS_STEP_INCOMPLETE;
990: while (!ts->reason && ros->status != TS_STEP_COMPLETE) {
991: const PetscReal h = ts->time_step;
992: for (i = 0; i < s; i++) {
993: ros->stage_time = ts->ptime + h * ASum[i];
994: TSPreStage(ts, ros->stage_time);
995: if (GammaZeroDiag[i]) {
996: ros->stage_explicit = PETSC_TRUE;
997: ros->scoeff = 1.;
998: } else {
999: ros->stage_explicit = PETSC_FALSE;
1000: ros->scoeff = 1. / Gamma[i * s + i];
1001: }
1003: VecCopy(ts->vec_sol, Zstage);
1004: for (j = 0; j < i; j++) w[j] = At[i * s + j];
1005: VecMAXPY(Zstage, i, w, Y);
1007: for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j];
1008: VecZeroEntries(Zdot);
1009: VecMAXPY(Zdot, i, w, Y);
1011: /* Initial guess taken from last stage */
1012: VecZeroEntries(Y[i]);
1014: if (!ros->stage_explicit) {
1015: TSGetSNES(ts, &snes);
1016: if (!ros->recompute_jacobian && !i) {
1017: SNESGetLagJacobian(snes, &lag);
1018: if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */
1019: SNESSetLagJacobian(snes, -2); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */
1020: }
1021: }
1022: SNESSolve(snes, NULL, Y[i]);
1023: if (!ros->recompute_jacobian && i == s - 1 && lag == 1) { SNESSetLagJacobian(snes, lag); /* Set lag back to 1 so we know user did not set it */ }
1024: SNESGetIterationNumber(snes, &its);
1025: SNESGetLinearSolveIterations(snes, &lits);
1026: ts->snes_its += its;
1027: ts->ksp_its += lits;
1028: } else {
1029: Mat J, Jp;
1030: VecZeroEntries(Ydot); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
1031: TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE);
1032: VecScale(Y[i], -1.0);
1033: VecAXPY(Y[i], -1.0, Zdot); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/
1035: VecZeroEntries(Zstage); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
1036: for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j];
1037: VecMAXPY(Zstage, i, w, Y);
1039: /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
1040: TSGetIJacobian(ts, &J, &Jp, NULL, NULL);
1041: TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE);
1042: MatMult(J, Zstage, Zdot);
1043: VecAXPY(Y[i], -1.0, Zdot);
1044: ts->ksp_its += 1;
1046: VecScale(Y[i], h);
1047: }
1048: TSPostStage(ts, ros->stage_time, i, Y);
1049: TSGetAdapt(ts, &adapt);
1050: TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok);
1051: if (!stageok) goto reject_step;
1052: }
1054: ros->status = TS_STEP_INCOMPLETE;
1055: TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL);
1056: ros->status = TS_STEP_PENDING;
1057: TSGetAdapt(ts, &adapt);
1058: TSAdaptCandidatesClear(adapt);
1059: TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE);
1060: TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
1061: ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1062: if (!accept) { /* Roll back the current step */
1063: TSRollBack_RosW(ts);
1064: ts->time_step = next_time_step;
1065: goto reject_step;
1066: }
1068: ts->ptime += ts->time_step;
1069: ts->time_step = next_time_step;
1070: break;
1072: reject_step:
1073: ts->reject++;
1074: accept = PETSC_FALSE;
1075: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1076: ts->reason = TS_DIVERGED_STEP_REJECTED;
1077: PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
1078: }
1079: }
1080: return 0;
1081: }
1083: static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U)
1084: {
1085: TS_RosW *ros = (TS_RosW *)ts->data;
1086: PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j;
1087: PetscReal h;
1088: PetscReal tt, t;
1089: PetscScalar *bt;
1090: const PetscReal *Bt = ros->tableau->binterpt;
1091: const PetscReal *GammaInv = ros->tableau->GammaInv;
1092: PetscScalar *w = ros->work;
1093: Vec *Y = ros->Y;
1097: switch (ros->status) {
1098: case TS_STEP_INCOMPLETE:
1099: case TS_STEP_PENDING:
1100: h = ts->time_step;
1101: t = (itime - ts->ptime) / h;
1102: break;
1103: case TS_STEP_COMPLETE:
1104: h = ts->ptime - ts->ptime_prev;
1105: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1106: break;
1107: default:
1108: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1109: }
1110: PetscMalloc1(s, &bt);
1111: for (i = 0; i < s; i++) bt[i] = 0;
1112: for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1113: for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt;
1114: }
1116: /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
1117: /* U <- 0*/
1118: VecZeroEntries(U);
1119: /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
1120: for (j = 0; j < s; j++) w[j] = 0;
1121: for (j = 0; j < s; j++) {
1122: for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j];
1123: }
1124: VecMAXPY(U, i, w, Y);
1125: /* U <- y(t) + U */
1126: VecAXPY(U, 1, ros->vec_sol_prev);
1128: PetscFree(bt);
1129: return 0;
1130: }
1132: /*------------------------------------------------------------*/
1134: static PetscErrorCode TSRosWTableauReset(TS ts)
1135: {
1136: TS_RosW *ros = (TS_RosW *)ts->data;
1137: RosWTableau tab = ros->tableau;
1139: if (!tab) return 0;
1140: VecDestroyVecs(tab->s, &ros->Y);
1141: PetscFree(ros->work);
1142: return 0;
1143: }
1145: static PetscErrorCode TSReset_RosW(TS ts)
1146: {
1147: TS_RosW *ros = (TS_RosW *)ts->data;
1149: TSRosWTableauReset(ts);
1150: VecDestroy(&ros->Ydot);
1151: VecDestroy(&ros->Ystage);
1152: VecDestroy(&ros->Zdot);
1153: VecDestroy(&ros->Zstage);
1154: VecDestroy(&ros->vec_sol_prev);
1155: return 0;
1156: }
1158: static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1159: {
1160: TS_RosW *rw = (TS_RosW *)ts->data;
1162: if (Ydot) {
1163: if (dm && dm != ts->dm) {
1164: DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot);
1165: } else *Ydot = rw->Ydot;
1166: }
1167: if (Zdot) {
1168: if (dm && dm != ts->dm) {
1169: DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot);
1170: } else *Zdot = rw->Zdot;
1171: }
1172: if (Ystage) {
1173: if (dm && dm != ts->dm) {
1174: DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage);
1175: } else *Ystage = rw->Ystage;
1176: }
1177: if (Zstage) {
1178: if (dm && dm != ts->dm) {
1179: DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage);
1180: } else *Zstage = rw->Zstage;
1181: }
1182: return 0;
1183: }
1185: static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1186: {
1187: if (Ydot) {
1188: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot);
1189: }
1190: if (Zdot) {
1191: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot);
1192: }
1193: if (Ystage) {
1194: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage);
1195: }
1196: if (Zstage) {
1197: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage);
1198: }
1199: return 0;
1200: }
1202: static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx)
1203: {
1204: return 0;
1205: }
1207: static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1208: {
1209: TS ts = (TS)ctx;
1210: Vec Ydot, Zdot, Ystage, Zstage;
1211: Vec Ydotc, Zdotc, Ystagec, Zstagec;
1213: TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage);
1214: TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec);
1215: MatRestrict(restrct, Ydot, Ydotc);
1216: VecPointwiseMult(Ydotc, rscale, Ydotc);
1217: MatRestrict(restrct, Ystage, Ystagec);
1218: VecPointwiseMult(Ystagec, rscale, Ystagec);
1219: MatRestrict(restrct, Zdot, Zdotc);
1220: VecPointwiseMult(Zdotc, rscale, Zdotc);
1221: MatRestrict(restrct, Zstage, Zstagec);
1222: VecPointwiseMult(Zstagec, rscale, Zstagec);
1223: TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage);
1224: TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec);
1225: return 0;
1226: }
1228: static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx)
1229: {
1230: return 0;
1231: }
1233: static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1234: {
1235: TS ts = (TS)ctx;
1236: Vec Ydot, Zdot, Ystage, Zstage;
1237: Vec Ydots, Zdots, Ystages, Zstages;
1239: TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage);
1240: TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages);
1242: VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD);
1243: VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD);
1245: VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD);
1246: VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD);
1248: VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD);
1249: VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD);
1251: VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD);
1252: VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD);
1254: TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage);
1255: TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages);
1256: return 0;
1257: }
1259: /*
1260: This defines the nonlinear equation that is to be solved with SNES
1261: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
1262: */
1263: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts)
1264: {
1265: TS_RosW *ros = (TS_RosW *)ts->data;
1266: Vec Ydot, Zdot, Ystage, Zstage;
1267: PetscReal shift = ros->scoeff / ts->time_step;
1268: DM dm, dmsave;
1270: SNESGetDM(snes, &dm);
1271: TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1272: VecWAXPY(Ydot, shift, U, Zdot); /* Ydot = shift*U + Zdot */
1273: VecWAXPY(Ystage, 1.0, U, Zstage); /* Ystage = U + Zstage */
1274: dmsave = ts->dm;
1275: ts->dm = dm;
1276: TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE);
1277: ts->dm = dmsave;
1278: TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1279: return 0;
1280: }
1282: static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts)
1283: {
1284: TS_RosW *ros = (TS_RosW *)ts->data;
1285: Vec Ydot, Zdot, Ystage, Zstage;
1286: PetscReal shift = ros->scoeff / ts->time_step;
1287: DM dm, dmsave;
1289: /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1290: SNESGetDM(snes, &dm);
1291: TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1292: dmsave = ts->dm;
1293: ts->dm = dm;
1294: TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE);
1295: ts->dm = dmsave;
1296: TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage);
1297: return 0;
1298: }
1300: static PetscErrorCode TSRosWTableauSetUp(TS ts)
1301: {
1302: TS_RosW *ros = (TS_RosW *)ts->data;
1303: RosWTableau tab = ros->tableau;
1305: VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y);
1306: PetscMalloc1(tab->s, &ros->work);
1307: return 0;
1308: }
1310: static PetscErrorCode TSSetUp_RosW(TS ts)
1311: {
1312: TS_RosW *ros = (TS_RosW *)ts->data;
1313: DM dm;
1314: SNES snes;
1315: TSRHSJacobian rhsjacobian;
1317: TSRosWTableauSetUp(ts);
1318: VecDuplicate(ts->vec_sol, &ros->Ydot);
1319: VecDuplicate(ts->vec_sol, &ros->Ystage);
1320: VecDuplicate(ts->vec_sol, &ros->Zdot);
1321: VecDuplicate(ts->vec_sol, &ros->Zstage);
1322: VecDuplicate(ts->vec_sol, &ros->vec_sol_prev);
1323: TSGetDM(ts, &dm);
1324: DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts);
1325: DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts);
1326: /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1327: TSGetSNES(ts, &snes);
1328: if (!((PetscObject)snes)->type_name) SNESSetType(snes, SNESKSPONLY);
1329: DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);
1330: if (rhsjacobian == TSComputeRHSJacobianConstant) {
1331: Mat Amat, Pmat;
1333: /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */
1334: SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL);
1335: if (Amat && Amat == ts->Arhs) {
1336: if (Amat == Pmat) {
1337: MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
1338: SNESSetJacobian(snes, Amat, Amat, NULL, NULL);
1339: } else {
1340: MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
1341: SNESSetJacobian(snes, Amat, NULL, NULL, NULL);
1342: if (Pmat && Pmat == ts->Brhs) {
1343: MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat);
1344: SNESSetJacobian(snes, NULL, Pmat, NULL, NULL);
1345: MatDestroy(&Pmat);
1346: }
1347: }
1348: MatDestroy(&Amat);
1349: }
1350: }
1351: return 0;
1352: }
1353: /*------------------------------------------------------------*/
1355: static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems *PetscOptionsObject)
1356: {
1357: TS_RosW *ros = (TS_RosW *)ts->data;
1358: SNES snes;
1360: PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options");
1361: {
1362: RosWTableauLink link;
1363: PetscInt count, choice;
1364: PetscBool flg;
1365: const char **namelist;
1367: for (link = RosWTableauList, count = 0; link; link = link->next, count++)
1368: ;
1369: PetscMalloc1(count, (char ***)&namelist);
1370: for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1371: PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg);
1372: if (flg) TSRosWSetType(ts, namelist[choice]);
1373: PetscFree(namelist);
1375: PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL);
1376: }
1377: PetscOptionsHeadEnd();
1378: /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1379: TSGetSNES(ts, &snes);
1380: if (!((PetscObject)snes)->type_name) SNESSetType(snes, SNESKSPONLY);
1381: return 0;
1382: }
1384: static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer)
1385: {
1386: TS_RosW *ros = (TS_RosW *)ts->data;
1387: PetscBool iascii;
1389: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1390: if (iascii) {
1391: RosWTableau tab = ros->tableau;
1392: TSRosWType rostype;
1393: char buf[512];
1394: PetscInt i;
1395: PetscReal abscissa[512];
1396: TSRosWGetType(ts, &rostype);
1397: PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype);
1398: PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum);
1399: PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf);
1400: for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i];
1401: PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa);
1402: PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf);
1403: }
1404: return 0;
1405: }
1407: static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer)
1408: {
1409: SNES snes;
1410: TSAdapt adapt;
1412: TSGetAdapt(ts, &adapt);
1413: TSAdaptLoad(adapt, viewer);
1414: TSGetSNES(ts, &snes);
1415: SNESLoad(snes, viewer);
1416: /* function and Jacobian context for SNES when used with TS is always ts object */
1417: SNESSetFunction(snes, NULL, NULL, ts);
1418: SNESSetJacobian(snes, NULL, NULL, NULL, ts);
1419: return 0;
1420: }
1422: /*@C
1423: TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme
1425: Logically collective
1427: Input Parameters:
1428: + ts - timestepping context
1429: - roswtype - type of Rosenbrock-W scheme
1431: Level: beginner
1433: .seealso: [](chapter_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3`
1434: @*/
1435: PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype)
1436: {
1439: PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype));
1440: return 0;
1441: }
1443: /*@C
1444: TSRosWGetType - Get the type of Rosenbrock-W scheme
1446: Logically collective
1448: Input Parameter:
1449: . ts - timestepping context
1451: Output Parameter:
1452: . rostype - type of Rosenbrock-W scheme
1454: Level: intermediate
1456: .seealso: [](chapter_ts), `TSRosWType`, `TSRosWSetType()`
1457: @*/
1458: PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype)
1459: {
1461: PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype));
1462: return 0;
1463: }
1465: /*@C
1466: TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step.
1468: Logically collective
1470: Input Parameters:
1471: + ts - timestepping context
1472: - flg - `PETSC_TRUE` to recompute the Jacobian at each stage
1474: Level: intermediate
1476: .seealso: [](chapter_ts), `TSRosWType`, `TSRosWGetType()`
1477: @*/
1478: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg)
1479: {
1481: PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg));
1482: return 0;
1483: }
1485: static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype)
1486: {
1487: TS_RosW *ros = (TS_RosW *)ts->data;
1489: *rostype = ros->tableau->name;
1490: return 0;
1491: }
1493: static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype)
1494: {
1495: TS_RosW *ros = (TS_RosW *)ts->data;
1496: PetscBool match;
1497: RosWTableauLink link;
1499: if (ros->tableau) {
1500: PetscStrcmp(ros->tableau->name, rostype, &match);
1501: if (match) return 0;
1502: }
1503: for (link = RosWTableauList; link; link = link->next) {
1504: PetscStrcmp(link->tab.name, rostype, &match);
1505: if (match) {
1506: if (ts->setupcalled) TSRosWTableauReset(ts);
1507: ros->tableau = &link->tab;
1508: if (ts->setupcalled) TSRosWTableauSetUp(ts);
1509: ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1510: return 0;
1511: }
1512: }
1513: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype);
1514: }
1516: static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg)
1517: {
1518: TS_RosW *ros = (TS_RosW *)ts->data;
1520: ros->recompute_jacobian = flg;
1521: return 0;
1522: }
1524: static PetscErrorCode TSDestroy_RosW(TS ts)
1525: {
1526: TSReset_RosW(ts);
1527: if (ts->dm) {
1528: DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts);
1529: DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts);
1530: }
1531: PetscFree(ts->data);
1532: PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL);
1533: PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL);
1534: PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL);
1535: return 0;
1536: }
1538: /* ------------------------------------------------------------ */
1539: /*MC
1540: TSROSW - ODE solver using Rosenbrock-W schemes
1542: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1543: nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1544: of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.
1546: Level: beginner
1548: Notes:
1549: This method currently only works with autonomous ODE and DAE.
1551: Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear.
1553: 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
1555: Developer Notes:
1556: Rosenbrock-W methods are typically specified for autonomous ODE
1558: $ udot = f(u)
1560: by the stage equations
1562: $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j
1564: and step completion formula
1566: $ u_1 = u_0 + sum_j b_j k_j
1568: with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u)
1569: and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner,
1570: we define new variables for the stage equations
1572: $ y_i = gamma_ij k_j
1574: The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define
1576: $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1}
1578: to rewrite the method as
1580: $ [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j
1581: $ u_1 = u_0 + sum_j bt_j y_j
1583: where we have introduced the mass matrix M. Continue by defining
1585: $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j
1587: or, more compactly in tensor notation
1589: $ Ydot = 1/h (Gamma^{-1} \otimes I) Y .
1591: Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current
1592: stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the
1593: equation
1595: $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0
1597: with initial guess y_i = 0.
1599: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`,
1600: `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType`
1601: M*/
1602: PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1603: {
1604: TS_RosW *ros;
1606: TSRosWInitializePackage();
1608: ts->ops->reset = TSReset_RosW;
1609: ts->ops->destroy = TSDestroy_RosW;
1610: ts->ops->view = TSView_RosW;
1611: ts->ops->load = TSLoad_RosW;
1612: ts->ops->setup = TSSetUp_RosW;
1613: ts->ops->step = TSStep_RosW;
1614: ts->ops->interpolate = TSInterpolate_RosW;
1615: ts->ops->evaluatestep = TSEvaluateStep_RosW;
1616: ts->ops->rollback = TSRollBack_RosW;
1617: ts->ops->setfromoptions = TSSetFromOptions_RosW;
1618: ts->ops->snesfunction = SNESTSFormFunction_RosW;
1619: ts->ops->snesjacobian = SNESTSFormJacobian_RosW;
1621: ts->usessnes = PETSC_TRUE;
1623: PetscNew(&ros);
1624: ts->data = (void *)ros;
1626: PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW);
1627: PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW);
1628: PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW);
1630: TSRosWSetType(ts, TSRosWDefault);
1631: return 0;
1632: }