Actual source code: rosw.c
1: /*
2: Code for timestepping with Rosenbrock W methods
4: Notes:
5: The general system is written as
7: F(t,U,Udot) = G(t,U)
9: where F represents the stiff part of the physics and G represents the non-stiff part.
10: This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian.
12: */
13: #include <petsc/private/tsimpl.h>
14: #include <petscdm.h>
16: #include <petsc/private/kernels/blockinvert.h>
18: static TSRosWType TSRosWDefault = TSROSWRA34PW2;
19: static PetscBool TSRosWRegisterAllCalled;
20: static PetscBool TSRosWPackageInitialized;
22: typedef struct _RosWTableau *RosWTableau;
23: struct _RosWTableau {
24: char *name;
25: PetscInt order; /* Classical approximation order of the method */
26: PetscInt s; /* Number of stages */
27: PetscInt pinterp; /* Interpolation order */
28: PetscReal *A; /* Propagation table, strictly lower triangular */
29: PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */
30: PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */
31: PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/
32: PetscReal *b; /* Step completion table */
33: PetscReal *bembed; /* Step completion table for embedded method of order one less */
34: PetscReal *ASum; /* Row sum of A */
35: PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */
36: PetscReal *At; /* Propagation table in transformed variables */
37: PetscReal *bt; /* Step completion table in transformed variables */
38: PetscReal *bembedt; /* Step completion table of order one less in transformed variables */
39: PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */
40: PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
41: PetscReal *binterpt; /* Dense output formula */
42: };
43: typedef struct _RosWTableauLink *RosWTableauLink;
44: struct _RosWTableauLink {
45: struct _RosWTableau tab;
46: RosWTableauLink next;
47: };
48: static RosWTableauLink RosWTableauList;
50: typedef struct {
51: RosWTableau tableau;
52: Vec *Y; /* States computed during the step, used to complete the step */
53: Vec Ydot; /* Work vector holding Ydot during residual evaluation */
54: Vec Ystage; /* Work vector for the state value at each stage */
55: Vec Zdot; /* Ydot = Zdot + shift*Y */
56: Vec Zstage; /* Y = Zstage + Y */
57: Vec vec_sol_prev; /* Solution from the previous step (used for interpolation and rollback)*/
58: PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */
59: PetscReal scoeff; /* shift = scoeff/dt */
60: PetscReal stage_time;
61: PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */
62: PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
63: TSStepStatus status;
64: } TS_RosW;
66: /*MC
67: TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method).
69: Only an approximate Jacobian is needed.
71: Level: intermediate
73: .seealso: [](ch_ts), `TSROSW`
74: M*/
76: /*MC
77: TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method).
79: Only an approximate Jacobian is needed.
81: Level: intermediate
83: .seealso: [](ch_ts), `TSROSW`
84: M*/
86: /*MC
87: TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme.
89: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P.
91: Level: intermediate
93: .seealso: [](ch_ts), `TSROSW`
94: M*/
96: /*MC
97: TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme.
99: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M.
101: Level: intermediate
103: .seealso: [](ch_ts), `TSROSW`
104: M*/
106: /*MC
107: TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`
109: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
111: This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73.
113: Level: intermediate
115: .seealso: [](ch_ts), `TSROSW`
116: M*/
118: /*MC
119: TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`.
121: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
123: This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48.
125: Level: intermediate
127: .seealso: [](ch_ts), `TSROSW`
128: M*/
130: /*MC
131: TSROSWR34PRW - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`.
133: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
135: This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25.
136: This method is B_{PR} consistent of order 3.
137: This method is spelled "ROS34PRw" in the paper, an improvement to an earlier "ROS34PRW" method from the same author with B_{PR} order 2.
139: Level: intermediate
141: .seealso: [](ch_ts), `TSROSW`
142: M*/
144: /*MC
145: TSROSWR3PRL2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`.
147: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
149: This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25.
150: This method is B_{PR} consistent of order 3.
152: Level: intermediate
154: .seealso: [](ch_ts), `TSROSW`
155: M*/
157: /*MC
158: TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`
160: By default, the Jacobian is only recomputed once per step.
162: Both the third order and embedded second order methods are stiffly accurate and L-stable.
164: Level: intermediate
166: .seealso: [](ch_ts), `TSROSW`, `TSROSWSANDU3`
167: M*/
169: /*MC
170: TSROSWRODASPR - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved`
172: By default, the Jacobian is only recomputed once per step.
174: Both the fourth order and embedded third order methods are stiffly accurate and L-stable.
175: The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems.
177: Level: intermediate
179: .seealso: [](ch_ts), `TSROSW`, `TSROSWR34PRW`, `TSROSWR3PRL2`
180: M*/
182: /*MC
183: TSROSWRODASPR2 - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved`
185: By default, the Jacobian is only recomputed once per step.
187: Both the fourth order and embedded third order methods are stiffly accurate and L-stable.
188: The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems.
189: This method is similar to `TSROSWRODASPR`, but satisfies one extra B_{PR} order condition.
191: Developer Note:
192: In numerical experiments with ts/tutorials/ex22.c, I (Jed) find this to produce surprisingly poor results.
193: Although the coefficients pass basic smoke tests, I'm not confident it was tabulated correctly in the paper.
194: It would be informative if someone could reproduce tests from the paper and/or reach out to the author to understand why it fails on this test problem.
195: If the method is implemented correctly, doing so might shed light on an additional analysis lens (or further conditions) for robustness on such problems.
197: Level: intermediate
199: .seealso: [](ch_ts), `TSROSW`, `TSROSWRODASPR`
200: M*/
202: /*MC
203: TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`
205: By default, the Jacobian is only recomputed once per step.
207: The third order method is L-stable, but not stiffly accurate.
208: The second order embedded method is strongly A-stable with R(infty) = 0.5.
209: The internal stages are L-stable.
210: This method is called ROS3 in {cite}`sandu_1997`.
212: Level: intermediate
214: .seealso: [](ch_ts), `TSROSW`, `TSROSWRODAS3`
215: M*/
217: /*MC
218: TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages
220: By default, the Jacobian is only recomputed once per step.
222: A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3)
224: Level: intermediate
226: .seealso: [](ch_ts), `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP`
227: M*/
229: /*MC
230: TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
232: By default, the Jacobian is only recomputed once per step.
234: L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
236: Level: intermediate
238: .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP`
239: M*/
241: /*MC
242: TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
244: By default, the Jacobian is only recomputed once per step.
246: L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
248: Level: intermediate
250: .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP`
251: M*/
253: /*MC
254: TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop {cite}`kaps1979generalized`
256: By default, the Jacobian is only recomputed once per step.
258: A(89.3 degrees)-stable, |R(infty)| = 0.454.
260: This method does not provide a dense output formula.
262: Level: intermediate
264: Note:
265: See Section 4 Table 7.2 in {cite}`wanner1996solving`
267: .seealso: [](ch_ts), `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`
268: M*/
270: /*MC
271: TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine {cite}`shampine1982implementation`
273: By default, the Jacobian is only recomputed once per step.
275: A-stable, |R(infty)| = 1/3.
277: This method does not provide a dense output formula.
279: Level: intermediate
281: Note:
282: See Section 4 Table 7.2 in {cite}`wanner1996solving`
284: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L`
285: M*/
287: /*MC
288: TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen {cite}`veldhuizen1984d`
290: By default, the Jacobian is only recomputed once per step.
292: A(89.5 degrees)-stable, |R(infty)| = 0.24.
294: This method does not provide a dense output formula.
296: Level: intermediate
298: Note:
299: See Section 4 Table 7.2 in {cite}`wanner1996solving`
301: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
302: M*/
304: /*MC
305: TSROSW4L - four stage, fourth order Rosenbrock (not W) method
307: By default, the Jacobian is only recomputed once per step.
309: A-stable and L-stable
311: This method does not provide a dense output formula.
313: Level: intermediate
315: Note:
316: See Section 4 Table 7.2 in {cite}`wanner1996solving`
318: .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
319: M*/
321: /*@C
322: TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in `TSROSW`
324: Not Collective, but should be called by all MPI processes which will need the schemes to be registered
326: Level: advanced
328: .seealso: [](ch_ts), `TSROSW`, `TSRosWRegisterDestroy()`
329: @*/
330: PetscErrorCode TSRosWRegisterAll(void)
331: {
332: PetscFunctionBegin;
333: if (TSRosWRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
334: TSRosWRegisterAllCalled = PETSC_TRUE;
336: {
337: const PetscReal A = 0;
338: const PetscReal Gamma = 1;
339: const PetscReal b = 1;
340: const PetscReal binterpt = 1;
342: PetscCall(TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
343: }
345: {
346: const PetscReal A = 0;
347: const PetscReal Gamma = 0.5;
348: const PetscReal b = 1;
349: const PetscReal binterpt = 1;
351: PetscCall(TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
352: }
354: {
355: /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */
356: const PetscReal A[2][2] = {
357: {0, 0},
358: {1., 0}
359: };
360: const PetscReal Gamma[2][2] = {
361: {1.707106781186547524401, 0 },
362: {-2. * 1.707106781186547524401, 1.707106781186547524401}
363: };
364: const PetscReal b[2] = {0.5, 0.5};
365: const PetscReal b1[2] = {1.0, 0.0};
366: PetscReal binterpt[2][2];
367: binterpt[0][0] = 1.707106781186547524401 - 1.0;
368: binterpt[1][0] = 2.0 - 1.707106781186547524401;
369: binterpt[0][1] = 1.707106781186547524401 - 1.5;
370: binterpt[1][1] = 1.5 - 1.707106781186547524401;
372: PetscCall(TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]));
373: }
374: {
375: /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */
376: const PetscReal A[2][2] = {
377: {0, 0},
378: {1., 0}
379: };
380: const PetscReal Gamma[2][2] = {
381: {0.2928932188134524755992, 0 },
382: {-2. * 0.2928932188134524755992, 0.2928932188134524755992}
383: };
384: const PetscReal b[2] = {0.5, 0.5};
385: const PetscReal b1[2] = {1.0, 0.0};
386: PetscReal binterpt[2][2];
387: binterpt[0][0] = 0.2928932188134524755992 - 1.0;
388: binterpt[1][0] = 2.0 - 0.2928932188134524755992;
389: binterpt[0][1] = 0.2928932188134524755992 - 1.5;
390: binterpt[1][1] = 1.5 - 0.2928932188134524755992;
392: PetscCall(TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0]));
393: }
394: {
395: /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */
396: PetscReal binterpt[3][2];
397: const PetscReal A[3][3] = {
398: {0, 0, 0},
399: {1.5773502691896257e+00, 0, 0},
400: {0.5, 0, 0}
401: };
402: const PetscReal Gamma[3][3] = {
403: {7.8867513459481287e-01, 0, 0 },
404: {-1.5773502691896257e+00, 7.8867513459481287e-01, 0 },
405: {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01}
406: };
407: const PetscReal b[3] = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01};
408: const PetscReal b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01};
410: binterpt[0][0] = -0.8094010767585034;
411: binterpt[1][0] = -0.5;
412: binterpt[2][0] = 2.3094010767585034;
413: binterpt[0][1] = 0.9641016151377548;
414: binterpt[1][1] = 0.5;
415: binterpt[2][1] = -1.4641016151377548;
417: PetscCall(TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
418: }
419: {
420: PetscReal binterpt[4][3];
421: /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
422: const PetscReal A[4][4] = {
423: {0, 0, 0, 0},
424: {8.7173304301691801e-01, 0, 0, 0},
425: {8.4457060015369423e-01, -1.1299064236484185e-01, 0, 0},
426: {0, 0, 1., 0}
427: };
428: const PetscReal Gamma[4][4] = {
429: {4.3586652150845900e-01, 0, 0, 0 },
430: {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 },
431: {-9.0338057013044082e-01, 5.4180672388095326e-02, 4.3586652150845900e-01, 0 },
432: {2.4212380706095346e-01, -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01}
433: };
434: const PetscReal b[4] = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01};
435: const PetscReal b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01};
437: binterpt[0][0] = 1.0564298455794094;
438: binterpt[1][0] = 2.296429974281067;
439: binterpt[2][0] = -1.307599564525376;
440: binterpt[3][0] = -1.045260255335102;
441: binterpt[0][1] = -1.3864882699759573;
442: binterpt[1][1] = -8.262611700275677;
443: binterpt[2][1] = 7.250979895056055;
444: binterpt[3][1] = 2.398120075195581;
445: binterpt[0][2] = 0.5721822314575016;
446: binterpt[1][2] = 4.742931142090097;
447: binterpt[2][2] = -4.398120075195578;
448: binterpt[3][2] = -0.9169932983520199;
450: PetscCall(TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
451: }
452: {
453: /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
454: const PetscReal A[4][4] = {
455: {0, 0, 0, 0},
456: {8.7173304301691801e-01, 0, 0, 0},
457: {1.4722022879435914e+00, -3.1840250568090289e-01, 0, 0},
458: {8.1505192016694938e-01, 5.0000000000000000e-01, -3.1505192016694938e-01, 0}
459: };
460: const PetscReal Gamma[4][4] = {
461: {4.3586652150845900e-01, 0, 0, 0 },
462: {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 },
463: {-1.2855347382089872e+00, 5.0507005541550687e-01, 4.3586652150845900e-01, 0 },
464: {-4.8201449182864348e-01, 2.1793326075422950e-01, -1.7178529043404503e-01, 4.3586652150845900e-01}
465: };
466: const PetscReal b[4] = {3.3303742833830591e-01, 7.1793326075422947e-01, -4.8683721060099439e-01, 4.3586652150845900e-01};
467: const PetscReal b2[4] = {2.5000000000000000e-01, 7.4276119608319180e-01, -3.1472922970066219e-01, 3.2196803361747034e-01};
469: PetscCall(TSRosWRegister(TSROSWR34PRW, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
470: }
471: {
472: /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */
473: const PetscReal A[4][4] = {
474: {0, 0, 0, 0},
475: {1.3075995645253771e+00, 0, 0, 0},
476: {0.5, 0.5, 0, 0},
477: {0.5, 0.5, 0, 0}
478: };
479: const PetscReal Gamma[4][4] = {
480: {4.3586652150845900e-01, 0, 0, 0 },
481: {-1.3075995645253771e+00, 4.3586652150845900e-01, 0, 0 },
482: {-7.0988575860972170e-01, -5.5996735960277766e-01, 4.3586652150845900e-01, 0 },
483: {-1.5550856807552085e-01, -9.5388516575112225e-01, 6.7352721231818413e-01, 4.3586652150845900e-01}
484: };
485: const PetscReal b[4] = {3.4449143192447917e-01, -4.5388516575112231e-01, 6.7352721231818413e-01, 4.3586652150845900e-01};
486: const PetscReal b2[4] = {5.0000000000000000e-01, -2.5738812086522078e-01, 4.3542008724775044e-01, 3.2196803361747034e-01};
488: PetscCall(TSRosWRegister(TSROSWR3PRL2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
489: }
490: {
491: /* const PetscReal g = 0.5; Directly written in-place below */
492: const PetscReal A[4][4] = {
493: {0, 0, 0, 0},
494: {0, 0, 0, 0},
495: {1., 0, 0, 0},
496: {0.75, -0.25, 0.5, 0}
497: };
498: const PetscReal Gamma[4][4] = {
499: {0.5, 0, 0, 0 },
500: {1., 0.5, 0, 0 },
501: {-0.25, -0.25, 0.5, 0 },
502: {1. / 12, 1. / 12, -2. / 3, 0.5}
503: };
504: const PetscReal b[4] = {5. / 6, -1. / 6, -1. / 6, 0.5};
505: const PetscReal b2[4] = {0.75, -0.25, 0.5, 0};
507: PetscCall(TSRosWRegister(TSROSWRODAS3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
508: }
509: {
510: /* const PetscReal g = 0.25; Directly written in-place below */
511: const PetscReal A[6][6] = {
512: {0, 0, 0, 0, 0, 0},
513: {0.75, 0, 0, 0, 0, 0},
514: {7.5162877593868457e-02, 2.4837122406131545e-02, 0, 0, 0, 0},
515: {1.6532708886396510e+00, 2.1545706385445562e-01, -1.3157488872766792e+00, 0, 0, 0},
516: {1.9385003738039885e+01, 1.2007117225835324e+00, -1.9337924059522791e+01, -2.4779140110062559e-01, 0, 0},
517: {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0}
518: };
519: const PetscReal Gamma[6][6] = {
520: {0.25, 0, 0, 0, 0, 0 },
521: {-7.5000000000000000e-01, 0.25, 0, 0, 0, 0 },
522: {-8.8644359075349941e-02, -2.8688974257983398e-02, 0.25, 0, 0, 0 },
523: {-4.8470034585330284e+00, -3.1583244269672095e-01, 4.9536568360123221e+00, 0.25, 0, 0 },
524: {-2.6769456904577400e+01, -1.5066459128852787e+00, 2.7200131480460591e+01, 8.2597133700208525e-01, 0.25, 0 },
525: {6.5876206496361416e+00, 3.6807059172993878e-01, -6.7423520694658121e+00, -1.0619631475741095e-01, -3.5714285714285715e-01, 0.25}
526: };
527: const PetscReal b[6] = {-7.9683251690137014e-01, 6.2136401428192344e-02, 1.1198553514719862e+00, 4.7198362114404874e-01, -1.0714285714285714e-01, 0.25};
528: const PetscReal b2[6] = {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0.0};
530: PetscCall(TSRosWRegister(TSROSWRODASPR, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
531: }
532: {
533: /* const PetscReal g = 0.3125; Directly written in-place below */
534: const PetscReal A[6][6] = {
535: {0, 0, 0, 0, 0, 0},
536: {9.3750000000000000e-01, 0, 0, 0, 0, 0},
537: {-4.7145892646261345e-02, 5.4531286650471122e-01, 0, 0, 0, 0},
538: {4.6915543899742240e-01, 4.4490537602383673e-01, -2.2498239334061121e-01, 0, 0, 0},
539: {1.0950372887345903e+00, 6.3223023457294381e-01, -8.9232966090485821e-01, 1.6506213759732410e-01, 0, 0},
540: {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 3.1250000000000000e-01, 0}
541: };
542: const PetscReal Gamma[6][6] = {
543: {0.3125, 0, 0, 0, 0, 0 },
544: {-9.3750000000000000e-01, 0.3125, 0, 0, 0, 0 },
545: {-9.7580572085994507e-02, -5.8666328499964138e-01, 0.3125, 0, 0, 0 },
546: {-4.9407065013256957e-01, -5.6819726428975503e-01, 5.0318949274167679e-01, 0.3125, 0, 0 },
547: {-1.2725031394709183e+00, -1.2146444240989676e+00, 1.5741357867872399e+00, 6.0051177678264578e-01, 0.3125, 0 },
548: {6.9690744901421153e-01, 6.2237005730756434e-01, -1.1553701989197045e+00, 1.8350029013386296e-01, -6.5990759753593431e-01, 0.3125}
549: };
550: const PetscReal b[6] = {5.1944159827788361e-01, 3.9955867781540699e-02, -4.7356407303732290e-01, 9.4907420451383284e-01, -3.4740759753593431e-01, 0.3125};
551: const PetscReal b2[6] = {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 0.3125, 0};
553: PetscCall(TSRosWRegister(TSROSWRODASPR2, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL));
554: }
555: {
556: /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */
557: const PetscReal A[3][3] = {
558: {0, 0, 0},
559: {0.43586652150845899941601945119356, 0, 0},
560: {0.43586652150845899941601945119356, 0, 0}
561: };
562: const PetscReal Gamma[3][3] = {
563: {0.43586652150845899941601945119356, 0, 0 },
564: {-0.19294655696029095575009695436041, 0.43586652150845899941601945119356, 0 },
565: {0, 1.74927148125794685173529749738960, 0.43586652150845899941601945119356}
566: };
567: const PetscReal b[3] = {-0.75457412385404315829818998646589, 1.94100407061964420292840123379419, -0.18642994676560104463021124732829};
568: const PetscReal b2[3] = {-1.53358745784149585370766523913002, 2.81745131148625772213931745457622, -0.28386385364476186843165221544619};
570: PetscReal binterpt[3][2];
571: binterpt[0][0] = 3.793692883777660870425141387941;
572: binterpt[1][0] = -2.918692883777660870425141387941;
573: binterpt[2][0] = 0.125;
574: binterpt[0][1] = -0.725741064379812106687651020584;
575: binterpt[1][1] = 0.559074397713145440020984353917;
576: binterpt[2][1] = 0.16666666666666666666666666666667;
578: PetscCall(TSRosWRegister(TSROSWSANDU3, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
579: }
580: {
581: /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0;
582: * Direct evaluation: s3 = 1.732050807568877293527;
583: * g = 0.7886751345948128822546;
584: * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */
585: const PetscReal A[3][3] = {
586: {0, 0, 0},
587: {1, 0, 0},
588: {0.25, 0.25, 0}
589: };
590: const PetscReal Gamma[3][3] = {
591: {0, 0, 0 },
592: {(-3.0 - 1.732050807568877293527) / 6.0, 0.7886751345948128822546, 0 },
593: {(-3.0 - 1.732050807568877293527) / 24.0, (-3.0 - 1.732050807568877293527) / 8.0, 0.7886751345948128822546}
594: };
595: const PetscReal b[3] = {1. / 6., 1. / 6., 2. / 3.};
596: const PetscReal b2[3] = {1. / 4., 1. / 4., 1. / 2.};
597: PetscReal binterpt[3][2];
599: binterpt[0][0] = 0.089316397477040902157517886164709;
600: binterpt[1][0] = -0.91068360252295909784248211383529;
601: binterpt[2][0] = 1.8213672050459181956849642276706;
602: binterpt[0][1] = 0.077350269189625764509148780501957;
603: binterpt[1][1] = 1.077350269189625764509148780502;
604: binterpt[2][1] = -1.1547005383792515290182975610039;
606: PetscCall(TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
607: }
609: {
610: const PetscReal A[4][4] = {
611: {0, 0, 0, 0},
612: {1. / 2., 0, 0, 0},
613: {1. / 2., 1. / 2., 0, 0},
614: {1. / 6., 1. / 6., 1. / 6., 0}
615: };
616: const PetscReal Gamma[4][4] = {
617: {1. / 2., 0, 0, 0},
618: {0.0, 1. / 4., 0, 0},
619: {-2., -2. / 3., 2. / 3., 0},
620: {1. / 2., 5. / 36., -2. / 9, 0}
621: };
622: const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.};
623: const PetscReal b2[4] = {1. / 8., 3. / 4., 1. / 8., 0};
624: PetscReal binterpt[4][3];
626: binterpt[0][0] = 6.25;
627: binterpt[1][0] = -30.25;
628: binterpt[2][0] = 1.75;
629: binterpt[3][0] = 23.25;
630: binterpt[0][1] = -9.75;
631: binterpt[1][1] = 58.75;
632: binterpt[2][1] = -3.25;
633: binterpt[3][1] = -45.75;
634: binterpt[0][2] = 3.6666666666666666666666666666667;
635: binterpt[1][2] = -28.333333333333333333333333333333;
636: binterpt[2][2] = 1.6666666666666666666666666666667;
637: binterpt[3][2] = 23.;
639: PetscCall(TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
640: }
642: {
643: const PetscReal A[4][4] = {
644: {0, 0, 0, 0},
645: {1. / 2., 0, 0, 0},
646: {1. / 2., 1. / 2., 0, 0},
647: {1. / 6., 1. / 6., 1. / 6., 0}
648: };
649: const PetscReal Gamma[4][4] = {
650: {1. / 2., 0, 0, 0},
651: {0.0, 3. / 4., 0, 0},
652: {-2. / 3., -23. / 9., 2. / 9., 0},
653: {1. / 18., 65. / 108., -2. / 27, 0}
654: };
655: const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.};
656: const PetscReal b2[4] = {3. / 16., 10. / 16., 3. / 16., 0};
657: PetscReal binterpt[4][3];
659: binterpt[0][0] = 1.6911764705882352941176470588235;
660: binterpt[1][0] = 3.6813725490196078431372549019608;
661: binterpt[2][0] = 0.23039215686274509803921568627451;
662: binterpt[3][0] = -4.6029411764705882352941176470588;
663: binterpt[0][1] = -0.95588235294117647058823529411765;
664: binterpt[1][1] = -6.2401960784313725490196078431373;
665: binterpt[2][1] = -0.31862745098039215686274509803922;
666: binterpt[3][1] = 7.5147058823529411764705882352941;
667: binterpt[0][2] = -0.56862745098039215686274509803922;
668: binterpt[1][2] = 2.7254901960784313725490196078431;
669: binterpt[2][2] = 0.25490196078431372549019607843137;
670: binterpt[3][2] = -2.4117647058823529411764705882353;
672: PetscCall(TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
673: }
675: {
676: PetscReal A[4][4], Gamma[4][4], b[4], b2[4];
677: PetscReal binterpt[4][3];
679: Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816;
680: Gamma[0][1] = 0;
681: Gamma[0][2] = 0;
682: Gamma[0][3] = 0;
683: Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476;
684: Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816;
685: Gamma[1][2] = 0;
686: Gamma[1][3] = 0;
687: Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903;
688: Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131;
689: Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816;
690: Gamma[2][3] = 0;
691: Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783;
692: Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984;
693: Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198;
694: Gamma[3][3] = 0;
696: A[0][0] = 0;
697: A[0][1] = 0;
698: A[0][2] = 0;
699: A[0][3] = 0;
700: A[1][0] = 0.8717330430169179988320388950590125027645343373957631;
701: A[1][1] = 0;
702: A[1][2] = 0;
703: A[1][3] = 0;
704: A[2][0] = 0.5275890119763004115618079766722914408876108660811028;
705: A[2][1] = 0.07241098802369958843819203208518599088698057726988732;
706: A[2][2] = 0;
707: A[2][3] = 0;
708: A[3][0] = 0.3990960076760701320627260685975778145384666450351314;
709: A[3][1] = -0.4375576546135194437228463747348862825846903771419953;
710: A[3][2] = 1.038461646937449311660120300601880176655352737312713;
711: A[3][3] = 0;
713: b[0] = 0.1876410243467238251612921333138006734899663569186926;
714: b[1] = -0.5952974735769549480478230473706443582188442040780541;
715: b[2] = 0.9717899277217721234705114616271378792182450260943198;
716: b[3] = 0.4358665215084589994160194475295062513822671686978816;
718: b2[0] = 0.2147402862233891404862383521089097657790734483804460;
719: b2[1] = -0.4851622638849390928209050538171743017757490232519684;
720: b2[2] = 0.8687250025203875511662123688667549217531982787600080;
721: b2[3] = 0.4016969751411624011684543450940068201770721128357014;
723: binterpt[0][0] = 2.2565812720167954547104627844105;
724: binterpt[1][0] = 1.349166413351089573796243820819;
725: binterpt[2][0] = -2.4695174540533503758652847586647;
726: binterpt[3][0] = -0.13623023131453465264142184656474;
727: binterpt[0][1] = -3.0826699111559187902922463354557;
728: binterpt[1][1] = -2.4689115685996042534544925650515;
729: binterpt[2][1] = 5.7428279814696677152129332773553;
730: binterpt[3][1] = -0.19124650171414467146619437684812;
731: binterpt[0][2] = 1.0137296634858471607430756831148;
732: binterpt[1][2] = 0.52444768167155973161042570784064;
733: binterpt[2][2] = -2.3015205996945452158771370439586;
734: binterpt[3][2] = 0.76334325453713832352363565300308;
736: PetscCall(TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
737: }
738: PetscCall(TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_CURRENT, PETSC_CURRENT, 0, -0.1282612945269037e+01));
739: PetscCall(TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_CURRENT, PETSC_CURRENT, 0, 125. / 108.));
740: PetscCall(TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_CURRENT, PETSC_CURRENT, 0, -1.355958941201148));
741: PetscCall(TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_CURRENT, PETSC_CURRENT, 0, -1.093502252409163));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: /*@C
746: TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`.
748: Not Collective
750: Level: advanced
752: .seealso: [](ch_ts), `TSRosWRegister()`, `TSRosWRegisterAll()`
753: @*/
754: PetscErrorCode TSRosWRegisterDestroy(void)
755: {
756: RosWTableauLink link;
758: PetscFunctionBegin;
759: while ((link = RosWTableauList)) {
760: RosWTableau t = &link->tab;
761: RosWTableauList = link->next;
762: PetscCall(PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum));
763: PetscCall(PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr));
764: PetscCall(PetscFree2(t->bembed, t->bembedt));
765: PetscCall(PetscFree(t->binterpt));
766: PetscCall(PetscFree(t->name));
767: PetscCall(PetscFree(link));
768: }
769: TSRosWRegisterAllCalled = PETSC_FALSE;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@C
774: TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called
775: from `TSInitializePackage()`.
777: Level: developer
779: .seealso: [](ch_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()`
780: @*/
781: PetscErrorCode TSRosWInitializePackage(void)
782: {
783: PetscFunctionBegin;
784: if (TSRosWPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
785: TSRosWPackageInitialized = PETSC_TRUE;
786: PetscCall(TSRosWRegisterAll());
787: PetscCall(PetscRegisterFinalize(TSRosWFinalizePackage));
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /*@C
792: TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is
793: called from `PetscFinalize()`.
795: Level: developer
797: .seealso: [](ch_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()`
798: @*/
799: PetscErrorCode TSRosWFinalizePackage(void)
800: {
801: PetscFunctionBegin;
802: TSRosWPackageInitialized = PETSC_FALSE;
803: PetscCall(TSRosWRegisterDestroy());
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*@C
808: TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
810: Not Collective, but the same schemes should be registered on all processes on which they will be used
812: Input Parameters:
813: + name - identifier for method
814: . order - approximation order of method
815: . s - number of stages, this is the dimension of the matrices below
816: . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular
817: . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal
818: . b - Step completion table (dimension s)
819: . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available)
820: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt
821: - binterpt - Coefficients of the interpolation formula (dimension s*pinterp)
823: Level: advanced
825: Note:
826: Several Rosenbrock W methods are provided, this function is only needed to create new methods.
828: .seealso: [](ch_ts), `TSROSW`
829: @*/
830: 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[])
831: {
832: RosWTableauLink link;
833: RosWTableau t;
834: PetscInt i, j, k;
835: PetscScalar *GammaInv;
837: PetscFunctionBegin;
838: PetscAssertPointer(name, 1);
839: PetscAssertPointer(A, 4);
840: PetscAssertPointer(Gamma, 5);
841: PetscAssertPointer(b, 6);
842: if (bembed) PetscAssertPointer(bembed, 7);
844: PetscCall(TSRosWInitializePackage());
845: PetscCall(PetscNew(&link));
846: t = &link->tab;
847: PetscCall(PetscStrallocpy(name, &t->name));
848: t->order = order;
849: t->s = s;
850: PetscCall(PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum));
851: PetscCall(PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr));
852: PetscCall(PetscArraycpy(t->A, A, s * s));
853: PetscCall(PetscArraycpy(t->Gamma, Gamma, s * s));
854: PetscCall(PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s));
855: PetscCall(PetscArraycpy(t->b, b, s));
856: if (bembed) {
857: PetscCall(PetscMalloc2(s, &t->bembed, s, &t->bembedt));
858: PetscCall(PetscArraycpy(t->bembed, bembed, s));
859: }
860: for (i = 0; i < s; i++) {
861: t->ASum[i] = 0;
862: t->GammaSum[i] = 0;
863: for (j = 0; j < s; j++) {
864: t->ASum[i] += A[i * s + j];
865: t->GammaSum[i] += Gamma[i * s + j];
866: }
867: }
868: PetscCall(PetscMalloc1(s * s, &GammaInv)); /* Need to use Scalar for inverse, then convert back to Real */
869: for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i];
870: for (i = 0; i < s; i++) {
871: if (Gamma[i * s + i] == 0.0) {
872: GammaInv[i * s + i] = 1.0;
873: t->GammaZeroDiag[i] = PETSC_TRUE;
874: } else {
875: t->GammaZeroDiag[i] = PETSC_FALSE;
876: }
877: }
879: switch (s) {
880: case 1:
881: GammaInv[0] = 1. / GammaInv[0];
882: break;
883: case 2:
884: PetscCall(PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL));
885: break;
886: case 3:
887: PetscCall(PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL));
888: break;
889: case 4:
890: PetscCall(PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL));
891: break;
892: case 5: {
893: PetscInt ipvt5[5];
894: MatScalar work5[5 * 5];
895: PetscCall(PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL));
896: break;
897: }
898: case 6:
899: PetscCall(PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL));
900: break;
901: case 7:
902: PetscCall(PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL));
903: break;
904: default:
905: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s);
906: }
907: for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
908: PetscCall(PetscFree(GammaInv));
910: for (i = 0; i < s; i++) {
911: for (k = 0; k < i + 1; k++) {
912: t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]);
913: for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]);
914: }
915: }
917: for (i = 0; i < s; i++) {
918: for (j = 0; j < s; j++) {
919: t->At[i * s + j] = 0;
920: for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j];
921: }
922: t->bt[i] = 0;
923: for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i];
924: if (bembed) {
925: t->bembedt[i] = 0;
926: for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i];
927: }
928: }
929: t->ccfl = 1.0; /* Fix this */
931: t->pinterp = pinterp;
932: PetscCall(PetscMalloc1(s * pinterp, &t->binterpt));
933: PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp));
934: link->next = RosWTableauList;
935: RosWTableauList = link;
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: /*@C
940: TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices
942: Not Collective, but the same schemes should be registered on all processes on which they will be used
944: Input Parameters:
945: + name - identifier for method
946: . gamma - leading coefficient (diagonal entry)
947: . a2 - design parameter, see Table 7.2 of {cite}`wanner1996solving`
948: . a3 - design parameter or `PETSC_DETERMINE` to satisfy one of the order five conditions (Eq 7.22)
949: . b3 - design parameter, see Table 7.2 of {cite}`wanner1996solving`
950: - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer
952: Level: developer
954: Notes:
955: This routine encodes the design of fourth order Rosenbrock methods as described in {cite}`wanner1996solving`
956: It is used here to implement several methods from the book and can be used to experiment with new methods.
957: It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions.
959: .seealso: [](ch_ts), `TSROSW`, `TSRosWRegister()`
960: @*/
961: PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4)
962: {
963: /* Declare numeric constants so they can be quad precision without being truncated at double */
964: 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;
965: PetscReal a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp;
966: PetscReal A[4][4], Gamma[4][4], b[4], bm[4];
967: PetscScalar M[3][3], rhs[3];
969: PetscFunctionBegin;
970: /* Step 1: choose Gamma (input) */
971: /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */
972: if (a3 == (PetscReal)PETSC_DEFAULT || a3 == (PetscReal)PETSC_DETERMINE) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */
973: a4 = a3; /* consequence of 7.20 */
975: /* Solve order conditions 7.15a, 7.15c, 7.15e */
976: M[0][0] = one;
977: M[0][1] = one;
978: M[0][2] = one; /* 7.15a */
979: M[1][0] = 0.0;
980: M[1][1] = a2 * a2;
981: M[1][2] = a4 * a4; /* 7.15c */
982: M[2][0] = 0.0;
983: M[2][1] = a2 * a2 * a2;
984: M[2][2] = a4 * a4 * a4; /* 7.15e */
985: rhs[0] = one - b3;
986: rhs[1] = one / three - a3 * a3 * b3;
987: rhs[2] = one / four - a3 * a3 * a3 * b3;
988: PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL));
989: b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
990: b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
991: b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);
993: /* Step 3 */
994: beta43 = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */
995: beta32beta2p = p44 / (b4 * beta43); /* 7.15h */
996: beta4jbetajp = (p32 - b3 * beta32beta2p) / b4;
997: M[0][0] = b2;
998: M[0][1] = b3;
999: M[0][2] = b4;
1000: M[1][0] = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp;
1001: M[1][1] = a2 * a2 * beta4jbetajp;
1002: M[1][2] = -a2 * a2 * beta32beta2p;
1003: M[2][0] = b4 * beta43 * a3 * a3 - p43;
1004: M[2][1] = -b4 * beta43 * a2 * a2;
1005: M[2][2] = 0;
1006: rhs[0] = one / two - gamma;
1007: rhs[1] = 0;
1008: rhs[2] = -a2 * a2 * p32;
1009: PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL));
1010: beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]);
1011: beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]);
1012: beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]);
1014: /* Step 4: back-substitute */
1015: beta32 = beta32beta2p / beta2p;
1016: beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p;
1018: /* Step 5: 7.15f and 7.20, then 7.16 */
1019: a43 = 0;
1020: a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p);
1021: a42 = a32;
1023: A[0][0] = 0;
1024: A[0][1] = 0;
1025: A[0][2] = 0;
1026: A[0][3] = 0;
1027: A[1][0] = a2;
1028: A[1][1] = 0;
1029: A[1][2] = 0;
1030: A[1][3] = 0;
1031: A[2][0] = a3 - a32;
1032: A[2][1] = a32;
1033: A[2][2] = 0;
1034: A[2][3] = 0;
1035: A[3][0] = a4 - a43 - a42;
1036: A[3][1] = a42;
1037: A[3][2] = a43;
1038: A[3][3] = 0;
1039: Gamma[0][0] = gamma;
1040: Gamma[0][1] = 0;
1041: Gamma[0][2] = 0;
1042: Gamma[0][3] = 0;
1043: Gamma[1][0] = beta2p - A[1][0];
1044: Gamma[1][1] = gamma;
1045: Gamma[1][2] = 0;
1046: Gamma[1][3] = 0;
1047: Gamma[2][0] = beta3p - beta32 - A[2][0];
1048: Gamma[2][1] = beta32 - A[2][1];
1049: Gamma[2][2] = gamma;
1050: Gamma[2][3] = 0;
1051: Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0];
1052: Gamma[3][1] = beta42 - A[3][1];
1053: Gamma[3][2] = beta43 - A[3][2];
1054: Gamma[3][3] = gamma;
1055: b[0] = b1;
1056: b[1] = b2;
1057: b[2] = b3;
1058: b[3] = b4;
1060: /* Construct embedded formula using given e4. We are solving Equation 7.18. */
1061: bm[3] = b[3] - e4 * gamma; /* using definition of E4 */
1062: bm[2] = (p32 - beta4jbetajp * bm[3]) / (beta32 * beta2p); /* fourth row of 7.18 */
1063: bm[1] = (one / two - gamma - beta3p * bm[2] - beta4p * bm[3]) / beta2p; /* second row */
1064: bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */
1066: {
1067: const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three;
1068: PetscCheck(PetscAbs(misfit) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "Assumptions violated, could not construct a third order embedded method");
1069: }
1070: PetscCall(TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*
1075: The step completion formula is
1077: x1 = x0 + b^T Y
1079: where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been
1080: updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write
1082: x1e = x0 + be^T Y
1083: = x1 - b^T Y + be^T Y
1084: = x1 + (be - b)^T Y
1086: so we can evaluate the method of different order even after the step has been optimistically completed.
1087: */
1088: static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done)
1089: {
1090: TS_RosW *ros = (TS_RosW *)ts->data;
1091: RosWTableau tab = ros->tableau;
1092: PetscScalar *w = ros->work;
1093: PetscInt i;
1095: PetscFunctionBegin;
1096: if (order == tab->order) {
1097: if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
1098: PetscCall(VecCopy(ts->vec_sol, U));
1099: for (i = 0; i < tab->s; i++) w[i] = tab->bt[i];
1100: PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
1101: } else PetscCall(VecCopy(ts->vec_sol, U));
1102: if (done) *done = PETSC_TRUE;
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: } else if (order == tab->order - 1) {
1105: if (!tab->bembedt) goto unavailable;
1106: if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
1107: PetscCall(VecCopy(ts->vec_sol, U));
1108: for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i];
1109: PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
1110: } else { /* Use rollback-and-recomplete formula (bembedt - bt) */
1111: for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
1112: PetscCall(VecCopy(ts->vec_sol, U));
1113: PetscCall(VecMAXPY(U, tab->s, w, ros->Y));
1114: }
1115: if (done) *done = PETSC_TRUE;
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1118: unavailable:
1119: PetscCheck(done, 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.",
1120: tab->name, tab->order, order);
1121: *done = PETSC_FALSE;
1122: PetscFunctionReturn(PETSC_SUCCESS);
1123: }
1125: static PetscErrorCode TSStep_RosW(TS ts)
1126: {
1127: TS_RosW *ros = (TS_RosW *)ts->data;
1128: RosWTableau tab = ros->tableau;
1129: const PetscInt s = tab->s;
1130: const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv;
1131: const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
1132: const PetscBool *GammaZeroDiag = tab->GammaZeroDiag;
1133: PetscScalar *w = ros->work;
1134: Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage;
1135: SNES snes;
1136: TSAdapt adapt;
1137: PetscInt i, j, its, lits;
1138: PetscInt rejections = 0;
1139: PetscBool stageok, accept = PETSC_TRUE;
1140: PetscReal next_time_step = ts->time_step;
1141: PetscInt lag;
1143: PetscFunctionBegin;
1144: if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev));
1146: ros->status = TS_STEP_INCOMPLETE;
1147: while (!ts->reason && ros->status != TS_STEP_COMPLETE) {
1148: const PetscReal h = ts->time_step;
1149: for (i = 0; i < s; i++) {
1150: ros->stage_time = ts->ptime + h * ASum[i];
1151: PetscCall(TSPreStage(ts, ros->stage_time));
1152: if (GammaZeroDiag[i]) {
1153: ros->stage_explicit = PETSC_TRUE;
1154: ros->scoeff = 1.;
1155: } else {
1156: ros->stage_explicit = PETSC_FALSE;
1157: ros->scoeff = 1. / Gamma[i * s + i];
1158: }
1160: PetscCall(VecCopy(ts->vec_sol, Zstage));
1161: for (j = 0; j < i; j++) w[j] = At[i * s + j];
1162: PetscCall(VecMAXPY(Zstage, i, w, Y));
1164: for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j];
1165: PetscCall(VecZeroEntries(Zdot));
1166: PetscCall(VecMAXPY(Zdot, i, w, Y));
1168: /* Initial guess taken from last stage */
1169: PetscCall(VecZeroEntries(Y[i]));
1171: if (!ros->stage_explicit) {
1172: PetscCall(TSGetSNES(ts, &snes));
1173: if (!ros->recompute_jacobian && !i) {
1174: PetscCall(SNESGetLagJacobian(snes, &lag));
1175: if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */
1176: PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */
1177: }
1178: }
1179: PetscCall(SNESSolve(snes, NULL, Y[i]));
1180: if (!ros->recompute_jacobian && i == s - 1 && lag == 1) PetscCall(SNESSetLagJacobian(snes, lag)); /* Set lag back to 1 so we know user did not set it */
1181: PetscCall(SNESGetIterationNumber(snes, &its));
1182: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1183: ts->snes_its += its;
1184: ts->ksp_its += lits;
1185: } else {
1186: Mat J, Jp;
1187: PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
1188: PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE));
1189: PetscCall(VecScale(Y[i], -1.0));
1190: PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/
1192: PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
1193: for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j];
1194: PetscCall(VecMAXPY(Zstage, i, w, Y));
1196: /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
1197: PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL));
1198: PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE));
1199: PetscCall(MatMult(J, Zstage, Zdot));
1200: PetscCall(VecAXPY(Y[i], -1.0, Zdot));
1201: ts->ksp_its += 1;
1203: PetscCall(VecScale(Y[i], h));
1204: }
1205: PetscCall(TSPostStage(ts, ros->stage_time, i, Y));
1206: PetscCall(TSGetAdapt(ts, &adapt));
1207: PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok));
1208: if (!stageok) goto reject_step;
1209: }
1211: ros->status = TS_STEP_INCOMPLETE;
1212: PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL));
1213: ros->status = TS_STEP_PENDING;
1214: PetscCall(TSGetAdapt(ts, &adapt));
1215: PetscCall(TSAdaptCandidatesClear(adapt));
1216: PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
1217: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
1218: ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
1219: if (!accept) { /* Roll back the current step */
1220: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
1221: ts->time_step = next_time_step;
1222: goto reject_step;
1223: }
1225: ts->ptime += ts->time_step;
1226: ts->time_step = next_time_step;
1227: break;
1229: reject_step:
1230: ts->reject++;
1231: accept = PETSC_FALSE;
1232: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
1233: ts->reason = TS_DIVERGED_STEP_REJECTED;
1234: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
1235: }
1236: }
1237: PetscFunctionReturn(PETSC_SUCCESS);
1238: }
1240: static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U)
1241: {
1242: TS_RosW *ros = (TS_RosW *)ts->data;
1243: PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j;
1244: PetscReal h;
1245: PetscReal tt, t;
1246: PetscScalar *bt;
1247: const PetscReal *Bt = ros->tableau->binterpt;
1248: const PetscReal *GammaInv = ros->tableau->GammaInv;
1249: PetscScalar *w = ros->work;
1250: Vec *Y = ros->Y;
1252: PetscFunctionBegin;
1253: PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name);
1255: switch (ros->status) {
1256: case TS_STEP_INCOMPLETE:
1257: case TS_STEP_PENDING:
1258: h = ts->time_step;
1259: t = (itime - ts->ptime) / h;
1260: break;
1261: case TS_STEP_COMPLETE:
1262: h = ts->ptime - ts->ptime_prev;
1263: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1264: break;
1265: default:
1266: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1267: }
1268: PetscCall(PetscMalloc1(s, &bt));
1269: for (i = 0; i < s; i++) bt[i] = 0;
1270: for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
1271: for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt;
1272: }
1274: /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
1275: /* U <- 0*/
1276: PetscCall(VecZeroEntries(U));
1277: /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
1278: for (j = 0; j < s; j++) w[j] = 0;
1279: for (j = 0; j < s; j++) {
1280: for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j];
1281: }
1282: PetscCall(VecMAXPY(U, i, w, Y));
1283: /* U <- y(t) + U */
1284: PetscCall(VecAXPY(U, 1, ros->vec_sol_prev));
1286: PetscCall(PetscFree(bt));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*------------------------------------------------------------*/
1292: static PetscErrorCode TSRosWTableauReset(TS ts)
1293: {
1294: TS_RosW *ros = (TS_RosW *)ts->data;
1295: RosWTableau tab = ros->tableau;
1297: PetscFunctionBegin;
1298: if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1299: PetscCall(VecDestroyVecs(tab->s, &ros->Y));
1300: PetscCall(PetscFree(ros->work));
1301: PetscFunctionReturn(PETSC_SUCCESS);
1302: }
1304: static PetscErrorCode TSReset_RosW(TS ts)
1305: {
1306: TS_RosW *ros = (TS_RosW *)ts->data;
1308: PetscFunctionBegin;
1309: PetscCall(TSRosWTableauReset(ts));
1310: PetscCall(VecDestroy(&ros->Ydot));
1311: PetscCall(VecDestroy(&ros->Ystage));
1312: PetscCall(VecDestroy(&ros->Zdot));
1313: PetscCall(VecDestroy(&ros->Zstage));
1314: PetscCall(VecDestroy(&ros->vec_sol_prev));
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1319: {
1320: TS_RosW *rw = (TS_RosW *)ts->data;
1322: PetscFunctionBegin;
1323: if (Ydot) {
1324: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1325: else *Ydot = rw->Ydot;
1326: }
1327: if (Zdot) {
1328: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1329: else *Zdot = rw->Zdot;
1330: }
1331: if (Ystage) {
1332: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1333: else *Ystage = rw->Ystage;
1334: }
1335: if (Zstage) {
1336: if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1337: else *Zstage = rw->Zstage;
1338: }
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1343: {
1344: PetscFunctionBegin;
1345: if (Ydot) {
1346: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1347: }
1348: if (Zdot) {
1349: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1350: }
1351: if (Ystage) {
1352: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1353: }
1354: if (Zstage) {
1355: if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1356: }
1357: PetscFunctionReturn(PETSC_SUCCESS);
1358: }
1360: static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx)
1361: {
1362: PetscFunctionBegin;
1363: PetscFunctionReturn(PETSC_SUCCESS);
1364: }
1366: static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
1367: {
1368: TS ts = (TS)ctx;
1369: Vec Ydot, Zdot, Ystage, Zstage;
1370: Vec Ydotc, Zdotc, Ystagec, Zstagec;
1372: PetscFunctionBegin;
1373: PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1374: PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1375: PetscCall(MatRestrict(restrct, Ydot, Ydotc));
1376: PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc));
1377: PetscCall(MatRestrict(restrct, Ystage, Ystagec));
1378: PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec));
1379: PetscCall(MatRestrict(restrct, Zdot, Zdotc));
1380: PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc));
1381: PetscCall(MatRestrict(restrct, Zstage, Zstagec));
1382: PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec));
1383: PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1384: PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx)
1389: {
1390: PetscFunctionBegin;
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
1395: {
1396: TS ts = (TS)ctx;
1397: Vec Ydot, Zdot, Ystage, Zstage;
1398: Vec Ydots, Zdots, Ystages, Zstages;
1400: PetscFunctionBegin;
1401: PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1402: PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));
1404: PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));
1405: PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));
1407: PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));
1408: PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));
1410: PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));
1411: PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));
1413: PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));
1414: PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));
1416: PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1417: PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts)
1422: {
1423: TS_RosW *ros = (TS_RosW *)ts->data;
1424: Vec Ydot, Zdot, Ystage, Zstage;
1425: PetscReal shift = ros->scoeff / ts->time_step;
1426: DM dm, dmsave;
1428: PetscFunctionBegin;
1429: PetscCall(SNESGetDM(snes, &dm));
1430: PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1431: PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */
1432: PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */
1433: dmsave = ts->dm;
1434: ts->dm = dm;
1435: PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE));
1436: ts->dm = dmsave;
1437: PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts)
1442: {
1443: TS_RosW *ros = (TS_RosW *)ts->data;
1444: Vec Ydot, Zdot, Ystage, Zstage;
1445: PetscReal shift = ros->scoeff / ts->time_step;
1446: DM dm, dmsave;
1448: PetscFunctionBegin;
1449: /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1450: PetscCall(SNESGetDM(snes, &dm));
1451: PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1452: dmsave = ts->dm;
1453: ts->dm = dm;
1454: PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE));
1455: ts->dm = dmsave;
1456: PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: static PetscErrorCode TSRosWTableauSetUp(TS ts)
1461: {
1462: TS_RosW *ros = (TS_RosW *)ts->data;
1463: RosWTableau tab = ros->tableau;
1465: PetscFunctionBegin;
1466: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y));
1467: PetscCall(PetscMalloc1(tab->s, &ros->work));
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: static PetscErrorCode TSSetUp_RosW(TS ts)
1472: {
1473: TS_RosW *ros = (TS_RosW *)ts->data;
1474: DM dm;
1475: SNES snes;
1476: TSRHSJacobianFn *rhsjacobian;
1478: PetscFunctionBegin;
1479: PetscCall(TSRosWTableauSetUp(ts));
1480: PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot));
1481: PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage));
1482: PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot));
1483: PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage));
1484: PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev));
1485: PetscCall(TSGetDM(ts, &dm));
1486: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1487: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1488: /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1489: PetscCall(TSGetSNES(ts, &snes));
1490: if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1491: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1492: if (rhsjacobian == TSComputeRHSJacobianConstant) {
1493: Mat Amat, Pmat;
1495: /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */
1496: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
1497: if (Amat && Amat == ts->Arhs) {
1498: if (Amat == Pmat) {
1499: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1500: PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL));
1501: } else {
1502: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1503: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
1504: if (Pmat && Pmat == ts->Brhs) {
1505: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
1506: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
1507: PetscCall(MatDestroy(&Pmat));
1508: }
1509: }
1510: PetscCall(MatDestroy(&Amat));
1511: }
1512: }
1513: PetscFunctionReturn(PETSC_SUCCESS);
1514: }
1515: /*------------------------------------------------------------*/
1517: static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject)
1518: {
1519: TS_RosW *ros = (TS_RosW *)ts->data;
1520: SNES snes;
1522: PetscFunctionBegin;
1523: PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options");
1524: {
1525: RosWTableauLink link;
1526: PetscInt count, choice;
1527: PetscBool flg;
1528: const char **namelist;
1530: for (link = RosWTableauList, count = 0; link; link = link->next, count++);
1531: PetscCall(PetscMalloc1(count, (char ***)&namelist));
1532: for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1533: PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg));
1534: if (flg) PetscCall(TSRosWSetType(ts, namelist[choice]));
1535: PetscCall(PetscFree(namelist));
1537: PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL));
1538: }
1539: PetscOptionsHeadEnd();
1540: /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1541: PetscCall(TSGetSNES(ts, &snes));
1542: if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer)
1547: {
1548: TS_RosW *ros = (TS_RosW *)ts->data;
1549: PetscBool isascii;
1551: PetscFunctionBegin;
1552: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1553: if (isascii) {
1554: RosWTableau tab = ros->tableau;
1555: TSRosWType rostype;
1556: char buf[512];
1557: PetscInt i;
1558: PetscReal abscissa[512];
1559: PetscCall(TSRosWGetType(ts, &rostype));
1560: PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype));
1561: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum));
1562: PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf));
1563: for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->GammaSum[i];
1564: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa));
1565: PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf));
1566: }
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer)
1571: {
1572: SNES snes;
1573: TSAdapt adapt;
1575: PetscFunctionBegin;
1576: PetscCall(TSGetAdapt(ts, &adapt));
1577: PetscCall(TSAdaptLoad(adapt, viewer));
1578: PetscCall(TSGetSNES(ts, &snes));
1579: PetscCall(SNESLoad(snes, viewer));
1580: /* function and Jacobian context for SNES when used with TS is always ts object */
1581: PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
1582: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
1583: PetscFunctionReturn(PETSC_SUCCESS);
1584: }
1586: /*@
1587: TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme
1589: Logically Collective
1591: Input Parameters:
1592: + ts - timestepping context
1593: - roswtype - type of Rosenbrock-W scheme
1595: Level: beginner
1597: .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3`
1598: @*/
1599: PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype)
1600: {
1601: PetscFunctionBegin;
1603: PetscAssertPointer(roswtype, 2);
1604: PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: /*@
1609: TSRosWGetType - Get the type of Rosenbrock-W scheme
1611: Logically Collective
1613: Input Parameter:
1614: . ts - timestepping context
1616: Output Parameter:
1617: . rostype - type of Rosenbrock-W scheme
1619: Level: intermediate
1621: .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()`
1622: @*/
1623: PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype)
1624: {
1625: PetscFunctionBegin;
1627: PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype));
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: /*@
1632: TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step.
1634: Logically Collective
1636: Input Parameters:
1637: + ts - timestepping context
1638: - flg - `PETSC_TRUE` to recompute the Jacobian at each stage
1640: Level: intermediate
1642: .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()`
1643: @*/
1644: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg)
1645: {
1646: PetscFunctionBegin;
1648: PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg));
1649: PetscFunctionReturn(PETSC_SUCCESS);
1650: }
1652: static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype)
1653: {
1654: TS_RosW *ros = (TS_RosW *)ts->data;
1656: PetscFunctionBegin;
1657: *rostype = ros->tableau->name;
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype)
1662: {
1663: TS_RosW *ros = (TS_RosW *)ts->data;
1664: PetscBool match;
1665: RosWTableauLink link;
1667: PetscFunctionBegin;
1668: if (ros->tableau) {
1669: PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match));
1670: if (match) PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1672: for (link = RosWTableauList; link; link = link->next) {
1673: PetscCall(PetscStrcmp(link->tab.name, rostype, &match));
1674: if (match) {
1675: if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts));
1676: ros->tableau = &link->tab;
1677: if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts));
1678: ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1681: }
1682: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype);
1683: }
1685: static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg)
1686: {
1687: TS_RosW *ros = (TS_RosW *)ts->data;
1689: PetscFunctionBegin;
1690: ros->recompute_jacobian = flg;
1691: PetscFunctionReturn(PETSC_SUCCESS);
1692: }
1694: static PetscErrorCode TSDestroy_RosW(TS ts)
1695: {
1696: PetscFunctionBegin;
1697: PetscCall(TSReset_RosW(ts));
1698: if (ts->dm) {
1699: PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1700: PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1701: }
1702: PetscCall(PetscFree(ts->data));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL));
1706: PetscFunctionReturn(PETSC_SUCCESS);
1707: }
1709: /* ------------------------------------------------------------ */
1710: /*MC
1711: TSROSW - ODE solver using Rosenbrock-W schemes
1713: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1714: nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1715: of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.
1717: Level: beginner
1719: Notes:
1720: This method currently only works with autonomous ODE and DAE.
1722: Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear.
1724: 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
1726: Developer Notes:
1727: Rosenbrock-W methods are typically specified for autonomous ODE
1728: .vb
1729: udot = f(u)
1730: .ve
1731: by the stage equations
1732: .vb
1733: k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j
1734: .ve
1735: and step completion formula
1736: .vb
1737: u_1 = u_0 + sum_j b_j k_j
1738: .ve
1739: with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u)
1740: and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner,
1741: we define new variables for the stage equations
1742: .vb
1743: y_i = gamma_ij k_j
1744: .ve
1745: The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define
1746: .vb
1747: A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1}
1748: .ve
1749: to rewrite the method as
1750: .vb
1751: [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j
1752: u_1 = u_0 + sum_j bt_j y_j
1753: .ve
1755: where we have introduced the mass matrix M. Continue by defining
1756: .vb
1757: ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j
1758: .ve
1759: or, more compactly in tensor notation
1760: .vb
1761: Ydot = 1/h (Gamma^{-1} \otimes I) Y .
1762: .ve
1763: Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current
1764: stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the
1765: equation
1766: .vb
1767: g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0
1768: .ve
1769: with initial guess y_i = 0.
1771: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`,
1772: `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType`
1773: M*/
1774: PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1775: {
1776: TS_RosW *ros;
1778: PetscFunctionBegin;
1779: PetscCall(TSRosWInitializePackage());
1781: ts->ops->reset = TSReset_RosW;
1782: ts->ops->destroy = TSDestroy_RosW;
1783: ts->ops->view = TSView_RosW;
1784: ts->ops->load = TSLoad_RosW;
1785: ts->ops->setup = TSSetUp_RosW;
1786: ts->ops->step = TSStep_RosW;
1787: ts->ops->interpolate = TSInterpolate_RosW;
1788: ts->ops->evaluatestep = TSEvaluateStep_RosW;
1789: ts->ops->setfromoptions = TSSetFromOptions_RosW;
1790: ts->ops->snesfunction = SNESTSFormFunction_RosW;
1791: ts->ops->snesjacobian = SNESTSFormJacobian_RosW;
1793: ts->usessnes = PETSC_TRUE;
1795: PetscCall(PetscNew(&ros));
1796: ts->data = (void *)ros;
1798: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW));
1799: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW));
1800: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW));
1802: PetscCall(TSRosWSetType(ts, TSRosWDefault));
1803: PetscFunctionReturn(PETSC_SUCCESS);
1804: }