1 /*
2 Code for timestepping with Rosenbrock W methods
3
4 Notes:
5 The general system is written as
6
7 F(t,U,Udot) = G(t,U)
8
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.
11
12 */
13 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
14 #include <petscdm.h>
15
16 #include <petsc/private/kernels/blockinvert.h>
17
18 static TSRosWType TSRosWDefault = TSROSWRA34PW2;
19 static PetscBool TSRosWRegisterAllCalled;
20 static PetscBool TSRosWPackageInitialized;
21
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;
49
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;
65
66 /*MC
67 TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method).
68
69 Only an approximate Jacobian is needed.
70
71 Level: intermediate
72
73 .seealso: [](ch_ts), `TSROSW`
74 M*/
75
76 /*MC
77 TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method).
78
79 Only an approximate Jacobian is needed.
80
81 Level: intermediate
82
83 .seealso: [](ch_ts), `TSROSW`
84 M*/
85
86 /*MC
87 TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme.
88
89 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P.
90
91 Level: intermediate
92
93 .seealso: [](ch_ts), `TSROSW`
94 M*/
95
96 /*MC
97 TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme.
98
99 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M.
100
101 Level: intermediate
102
103 .seealso: [](ch_ts), `TSROSW`
104 M*/
105
106 /*MC
107 TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`
108
109 Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
110
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.
112
113 Level: intermediate
114
115 .seealso: [](ch_ts), `TSROSW`
116 M*/
117
118 /*MC
119 TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`.
120
121 Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
122
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.
124
125 Level: intermediate
126
127 .seealso: [](ch_ts), `TSROSW`
128 M*/
129
130 /*MC
131 TSROSWR34PRW - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`.
132
133 Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
134
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.
138
139 Level: intermediate
140
141 .seealso: [](ch_ts), `TSROSW`
142 M*/
143
144 /*MC
145 TSROSWR3PRL2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`.
146
147 Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
148
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.
151
152 Level: intermediate
153
154 .seealso: [](ch_ts), `TSROSW`
155 M*/
156
157 /*MC
158 TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`
159
160 By default, the Jacobian is only recomputed once per step.
161
162 Both the third order and embedded second order methods are stiffly accurate and L-stable.
163
164 Level: intermediate
165
166 .seealso: [](ch_ts), `TSROSW`, `TSROSWSANDU3`
167 M*/
168
169 /*MC
170 TSROSWRODASPR - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved`
171
172 By default, the Jacobian is only recomputed once per step.
173
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.
176
177 Level: intermediate
178
179 .seealso: [](ch_ts), `TSROSW`, `TSROSWR34PRW`, `TSROSWR3PRL2`
180 M*/
181
182 /*MC
183 TSROSWRODASPR2 - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved`
184
185 By default, the Jacobian is only recomputed once per step.
186
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.
190
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.
196
197 Level: intermediate
198
199 .seealso: [](ch_ts), `TSROSW`, `TSROSWRODASPR`
200 M*/
201
202 /*MC
203 TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme {cite}`sandu_1997`
204
205 By default, the Jacobian is only recomputed once per step.
206
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`.
211
212 Level: intermediate
213
214 .seealso: [](ch_ts), `TSROSW`, `TSROSWRODAS3`
215 M*/
216
217 /*MC
218 TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages
219
220 By default, the Jacobian is only recomputed once per step.
221
222 A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3)
223
224 Level: intermediate
225
226 .seealso: [](ch_ts), `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP`
227 M*/
228
229 /*MC
230 TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
231
232 By default, the Jacobian is only recomputed once per step.
233
234 L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
235
236 Level: intermediate
237
238 .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP`
239 M*/
240
241 /*MC
242 TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
243
244 By default, the Jacobian is only recomputed once per step.
245
246 L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
247
248 Level: intermediate
249
250 .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP`
251 M*/
252
253 /*MC
254 TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop {cite}`kaps1979generalized`
255
256 By default, the Jacobian is only recomputed once per step.
257
258 A(89.3 degrees)-stable, |R(infty)| = 0.454.
259
260 This method does not provide a dense output formula.
261
262 Level: intermediate
263
264 Note:
265 See Section 4 Table 7.2 in {cite}`wanner1996solving`
266
267 .seealso: [](ch_ts), `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`
268 M*/
269
270 /*MC
271 TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine {cite}`shampine1982implementation`
272
273 By default, the Jacobian is only recomputed once per step.
274
275 A-stable, |R(infty)| = 1/3.
276
277 This method does not provide a dense output formula.
278
279 Level: intermediate
280
281 Note:
282 See Section 4 Table 7.2 in {cite}`wanner1996solving`
283
284 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L`
285 M*/
286
287 /*MC
288 TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen {cite}`veldhuizen1984d`
289
290 By default, the Jacobian is only recomputed once per step.
291
292 A(89.5 degrees)-stable, |R(infty)| = 0.24.
293
294 This method does not provide a dense output formula.
295
296 Level: intermediate
297
298 Note:
299 See Section 4 Table 7.2 in {cite}`wanner1996solving`
300
301 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
302 M*/
303
304 /*MC
305 TSROSW4L - four stage, fourth order Rosenbrock (not W) method
306
307 By default, the Jacobian is only recomputed once per step.
308
309 A-stable and L-stable
310
311 This method does not provide a dense output formula.
312
313 Level: intermediate
314
315 Note:
316 See Section 4 Table 7.2 in {cite}`wanner1996solving`
317
318 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L`
319 M*/
320
321 /*@C
322 TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in `TSROSW`
323
324 Not Collective, but should be called by all MPI processes which will need the schemes to be registered
325
326 Level: advanced
327
328 .seealso: [](ch_ts), `TSROSW`, `TSRosWRegisterDestroy()`
329 @*/
TSRosWRegisterAll(void)330 PetscErrorCode TSRosWRegisterAll(void)
331 {
332 PetscFunctionBegin;
333 if (TSRosWRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
334 TSRosWRegisterAllCalled = PETSC_TRUE;
335
336 {
337 const PetscReal A = 0;
338 const PetscReal Gamma = 1;
339 const PetscReal b = 1;
340 const PetscReal binterpt = 1;
341
342 PetscCall(TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
343 }
344
345 {
346 const PetscReal A = 0;
347 const PetscReal Gamma = 0.5;
348 const PetscReal b = 1;
349 const PetscReal binterpt = 1;
350
351 PetscCall(TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt));
352 }
353
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;
371
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;
391
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};
409
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;
416
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};
436
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;
449
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};
468
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};
487
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};
506
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};
529
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};
552
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};
569
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;
577
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];
598
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;
605
606 PetscCall(TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0]));
607 }
608
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];
625
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.;
638
639 PetscCall(TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
640 }
641
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];
658
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;
671
672 PetscCall(TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0]));
673 }
674
675 {
676 PetscReal A[4][4], Gamma[4][4], b[4], b2[4];
677 PetscReal binterpt[4][3];
678
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;
695
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;
712
713 b[0] = 0.1876410243467238251612921333138006734899663569186926;
714 b[1] = -0.5952974735769549480478230473706443582188442040780541;
715 b[2] = 0.9717899277217721234705114616271378792182450260943198;
716 b[3] = 0.4358665215084589994160194475295062513822671686978816;
717
718 b2[0] = 0.2147402862233891404862383521089097657790734483804460;
719 b2[1] = -0.4851622638849390928209050538171743017757490232519684;
720 b2[2] = 0.8687250025203875511662123688667549217531982787600080;
721 b2[3] = 0.4016969751411624011684543450940068201770721128357014;
722
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;
735
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 }
744
745 /*@C
746 TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`.
747
748 Not Collective
749
750 Level: advanced
751
752 .seealso: [](ch_ts), `TSRosWRegister()`, `TSRosWRegisterAll()`
753 @*/
TSRosWRegisterDestroy(void)754 PetscErrorCode TSRosWRegisterDestroy(void)
755 {
756 RosWTableauLink link;
757
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 }
772
773 /*@C
774 TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called
775 from `TSInitializePackage()`.
776
777 Level: developer
778
779 .seealso: [](ch_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()`
780 @*/
TSRosWInitializePackage(void)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 }
790
791 /*@C
792 TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is
793 called from `PetscFinalize()`.
794
795 Level: developer
796
797 .seealso: [](ch_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()`
798 @*/
TSRosWFinalizePackage(void)799 PetscErrorCode TSRosWFinalizePackage(void)
800 {
801 PetscFunctionBegin;
802 TSRosWPackageInitialized = PETSC_FALSE;
803 PetscCall(TSRosWRegisterDestroy());
804 PetscFunctionReturn(PETSC_SUCCESS);
805 }
806
807 /*@C
808 TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
809
810 Not Collective, but the same schemes should be registered on all processes on which they will be used
811
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)
822
823 Level: advanced
824
825 Note:
826 Several Rosenbrock W methods are provided, this function is only needed to create new methods.
827
828 .seealso: [](ch_ts), `TSROSW`
829 @*/
TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[],PetscInt pinterp,const PetscReal binterpt[])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;
836
837 PetscFunctionBegin;
838 PetscAssertPointer(name, 1);
839 PetscAssertPointer(A, 4);
840 PetscAssertPointer(Gamma, 5);
841 PetscAssertPointer(b, 6);
842 if (bembed) PetscAssertPointer(bembed, 7);
843
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 }
878
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));
909
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 }
916
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 */
930
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 }
938
939 /*@C
940 TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices
941
942 Not Collective, but the same schemes should be registered on all processes on which they will be used
943
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
951
952 Level: developer
953
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.
958
959 .seealso: [](ch_ts), `TSROSW`, `TSRosWRegister()`
960 @*/
TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4)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];
968
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 */
974
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]);
992
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]);
1013
1014 /* Step 4: back-substitute */
1015 beta32 = beta32beta2p / beta2p;
1016 beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p;
1017
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;
1022
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;
1059
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 */
1065
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 }
1073
1074 /*
1075 The step completion formula is
1076
1077 x1 = x0 + b^T Y
1078
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
1081
1082 x1e = x0 + be^T Y
1083 = x1 - b^T Y + be^T Y
1084 = x1 + (be - b)^T Y
1085
1086 so we can evaluate the method of different order even after the step has been optimistically completed.
1087 */
TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool * done)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;
1094
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 }
1124
TSStep_RosW(TS ts)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;
1142
1143 PetscFunctionBegin;
1144 if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev));
1145
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 }
1159
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));
1163
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));
1167
1168 /* Initial guess taken from last stage */
1169 PetscCall(VecZeroEntries(Y[i]));
1170
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]*/
1191
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));
1195
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;
1202
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 }
1210
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 }
1224
1225 ts->ptime += ts->time_step;
1226 ts->time_step = next_time_step;
1227 break;
1228
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 }
1239
TSInterpolate_RosW(TS ts,PetscReal itime,Vec U)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;
1251
1252 PetscFunctionBegin;
1253 PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name);
1254
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 }
1273
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));
1285
1286 PetscCall(PetscFree(bt));
1287 PetscFunctionReturn(PETSC_SUCCESS);
1288 }
1289
TSRosWTableauReset(TS ts)1290 static PetscErrorCode TSRosWTableauReset(TS ts)
1291 {
1292 TS_RosW *ros = (TS_RosW *)ts->data;
1293 RosWTableau tab = ros->tableau;
1294
1295 PetscFunctionBegin;
1296 if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1297 PetscCall(VecDestroyVecs(tab->s, &ros->Y));
1298 PetscCall(PetscFree(ros->work));
1299 PetscFunctionReturn(PETSC_SUCCESS);
1300 }
1301
TSReset_RosW(TS ts)1302 static PetscErrorCode TSReset_RosW(TS ts)
1303 {
1304 TS_RosW *ros = (TS_RosW *)ts->data;
1305
1306 PetscFunctionBegin;
1307 PetscCall(TSRosWTableauReset(ts));
1308 PetscCall(VecDestroy(&ros->Ydot));
1309 PetscCall(VecDestroy(&ros->Ystage));
1310 PetscCall(VecDestroy(&ros->Zdot));
1311 PetscCall(VecDestroy(&ros->Zstage));
1312 PetscCall(VecDestroy(&ros->vec_sol_prev));
1313 PetscFunctionReturn(PETSC_SUCCESS);
1314 }
1315
TSRosWGetVecs(TS ts,DM dm,Vec * Ydot,Vec * Zdot,Vec * Ystage,Vec * Zstage)1316 static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1317 {
1318 TS_RosW *rw = (TS_RosW *)ts->data;
1319
1320 PetscFunctionBegin;
1321 if (Ydot) {
1322 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1323 else *Ydot = rw->Ydot;
1324 }
1325 if (Zdot) {
1326 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1327 else *Zdot = rw->Zdot;
1328 }
1329 if (Ystage) {
1330 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1331 else *Ystage = rw->Ystage;
1332 }
1333 if (Zstage) {
1334 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1335 else *Zstage = rw->Zstage;
1336 }
1337 PetscFunctionReturn(PETSC_SUCCESS);
1338 }
1339
TSRosWRestoreVecs(TS ts,DM dm,Vec * Ydot,Vec * Zdot,Vec * Ystage,Vec * Zstage)1340 static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage)
1341 {
1342 PetscFunctionBegin;
1343 if (Ydot) {
1344 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot));
1345 }
1346 if (Zdot) {
1347 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot));
1348 }
1349 if (Ystage) {
1350 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage));
1351 }
1352 if (Zstage) {
1353 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage));
1354 }
1355 PetscFunctionReturn(PETSC_SUCCESS);
1356 }
1357
DMCoarsenHook_TSRosW(DM fine,DM coarse,PetscCtx ctx)1358 static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, PetscCtx ctx)
1359 {
1360 PetscFunctionBegin;
1361 PetscFunctionReturn(PETSC_SUCCESS);
1362 }
1363
DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)1364 static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
1365 {
1366 TS ts = (TS)ctx;
1367 Vec Ydot, Zdot, Ystage, Zstage;
1368 Vec Ydotc, Zdotc, Ystagec, Zstagec;
1369
1370 PetscFunctionBegin;
1371 PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1372 PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1373 PetscCall(MatRestrict(restrct, Ydot, Ydotc));
1374 PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc));
1375 PetscCall(MatRestrict(restrct, Ystage, Ystagec));
1376 PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec));
1377 PetscCall(MatRestrict(restrct, Zdot, Zdotc));
1378 PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc));
1379 PetscCall(MatRestrict(restrct, Zstage, Zstagec));
1380 PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec));
1381 PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage));
1382 PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec));
1383 PetscFunctionReturn(PETSC_SUCCESS);
1384 }
1385
DMSubDomainHook_TSRosW(DM fine,DM coarse,PetscCtx ctx)1386 static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, PetscCtx ctx)
1387 {
1388 PetscFunctionBegin;
1389 PetscFunctionReturn(PETSC_SUCCESS);
1390 }
1391
DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)1392 static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
1393 {
1394 TS ts = (TS)ctx;
1395 Vec Ydot, Zdot, Ystage, Zstage;
1396 Vec Ydots, Zdots, Ystages, Zstages;
1397
1398 PetscFunctionBegin;
1399 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1400 PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));
1401
1402 PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));
1403 PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD));
1404
1405 PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));
1406 PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD));
1407
1408 PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));
1409 PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD));
1410
1411 PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));
1412 PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD));
1413
1414 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage));
1415 PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages));
1416 PetscFunctionReturn(PETSC_SUCCESS);
1417 }
1418
SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts)1419 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts)
1420 {
1421 TS_RosW *ros = (TS_RosW *)ts->data;
1422 Vec Ydot, Zdot, Ystage, Zstage;
1423 PetscReal shift = ros->scoeff / ts->time_step;
1424 DM dm, dmsave;
1425
1426 PetscFunctionBegin;
1427 PetscCall(SNESGetDM(snes, &dm));
1428 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1429 PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */
1430 PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */
1431 dmsave = ts->dm;
1432 ts->dm = dm;
1433 PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE));
1434 ts->dm = dmsave;
1435 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1436 PetscFunctionReturn(PETSC_SUCCESS);
1437 }
1438
SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat A,Mat B,TS ts)1439 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts)
1440 {
1441 TS_RosW *ros = (TS_RosW *)ts->data;
1442 Vec Ydot, Zdot, Ystage, Zstage;
1443 PetscReal shift = ros->scoeff / ts->time_step;
1444 DM dm, dmsave;
1445
1446 PetscFunctionBegin;
1447 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1448 PetscCall(SNESGetDM(snes, &dm));
1449 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1450 dmsave = ts->dm;
1451 ts->dm = dm;
1452 PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE));
1453 ts->dm = dmsave;
1454 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage));
1455 PetscFunctionReturn(PETSC_SUCCESS);
1456 }
1457
TSRosWTableauSetUp(TS ts)1458 static PetscErrorCode TSRosWTableauSetUp(TS ts)
1459 {
1460 TS_RosW *ros = (TS_RosW *)ts->data;
1461 RosWTableau tab = ros->tableau;
1462
1463 PetscFunctionBegin;
1464 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y));
1465 PetscCall(PetscMalloc1(tab->s, &ros->work));
1466 PetscFunctionReturn(PETSC_SUCCESS);
1467 }
1468
TSSetUp_RosW(TS ts)1469 static PetscErrorCode TSSetUp_RosW(TS ts)
1470 {
1471 TS_RosW *ros = (TS_RosW *)ts->data;
1472 DM dm;
1473 SNES snes;
1474 TSRHSJacobianFn *rhsjacobian;
1475
1476 PetscFunctionBegin;
1477 PetscCall(TSRosWTableauSetUp(ts));
1478 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot));
1479 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage));
1480 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot));
1481 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage));
1482 PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev));
1483 PetscCall(TSGetDM(ts, &dm));
1484 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1485 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1486 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1487 PetscCall(TSGetSNES(ts, &snes));
1488 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1489 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1490 if (rhsjacobian == TSComputeRHSJacobianConstant) {
1491 Mat Amat, Pmat;
1492
1493 /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */
1494 PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
1495 if (Amat && Amat == ts->Arhs) {
1496 if (Amat == Pmat) {
1497 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1498 PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL));
1499 } else {
1500 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
1501 PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
1502 if (Pmat && Pmat == ts->Brhs) {
1503 PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
1504 PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
1505 PetscCall(MatDestroy(&Pmat));
1506 }
1507 }
1508 PetscCall(MatDestroy(&Amat));
1509 }
1510 }
1511 PetscFunctionReturn(PETSC_SUCCESS);
1512 }
1513
TSSetFromOptions_RosW(TS ts,PetscOptionItems PetscOptionsObject)1514 static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject)
1515 {
1516 TS_RosW *ros = (TS_RosW *)ts->data;
1517 SNES snes;
1518
1519 PetscFunctionBegin;
1520 PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options");
1521 {
1522 RosWTableauLink link;
1523 PetscInt count, choice;
1524 PetscBool flg;
1525 const char **namelist;
1526
1527 for (link = RosWTableauList, count = 0; link; link = link->next, count++);
1528 PetscCall(PetscMalloc1(count, (char ***)&namelist));
1529 for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1530 PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg));
1531 if (flg) PetscCall(TSRosWSetType(ts, namelist[choice]));
1532 PetscCall(PetscFree(namelist));
1533
1534 PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL));
1535 }
1536 PetscOptionsHeadEnd();
1537 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1538 PetscCall(TSGetSNES(ts, &snes));
1539 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY));
1540 PetscFunctionReturn(PETSC_SUCCESS);
1541 }
1542
TSView_RosW(TS ts,PetscViewer viewer)1543 static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer)
1544 {
1545 TS_RosW *ros = (TS_RosW *)ts->data;
1546 PetscBool isascii;
1547
1548 PetscFunctionBegin;
1549 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1550 if (isascii) {
1551 RosWTableau tab = ros->tableau;
1552 TSRosWType rostype;
1553 char buf[512];
1554 PetscInt i;
1555 PetscReal abscissa[512];
1556 PetscCall(TSRosWGetType(ts, &rostype));
1557 PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype));
1558 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum));
1559 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf));
1560 for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->GammaSum[i];
1561 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa));
1562 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf));
1563 }
1564 PetscFunctionReturn(PETSC_SUCCESS);
1565 }
1566
TSLoad_RosW(TS ts,PetscViewer viewer)1567 static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer)
1568 {
1569 SNES snes;
1570 TSAdapt adapt;
1571
1572 PetscFunctionBegin;
1573 PetscCall(TSGetAdapt(ts, &adapt));
1574 PetscCall(TSAdaptLoad(adapt, viewer));
1575 PetscCall(TSGetSNES(ts, &snes));
1576 PetscCall(SNESLoad(snes, viewer));
1577 /* function and Jacobian context for SNES when used with TS is always ts object */
1578 PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
1579 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
1580 PetscFunctionReturn(PETSC_SUCCESS);
1581 }
1582
1583 /*@
1584 TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme
1585
1586 Logically Collective
1587
1588 Input Parameters:
1589 + ts - timestepping context
1590 - roswtype - type of Rosenbrock-W scheme
1591
1592 Level: beginner
1593
1594 .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3`
1595 @*/
TSRosWSetType(TS ts,TSRosWType roswtype)1596 PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype)
1597 {
1598 PetscFunctionBegin;
1599 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
1600 PetscAssertPointer(roswtype, 2);
1601 PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype));
1602 PetscFunctionReturn(PETSC_SUCCESS);
1603 }
1604
1605 /*@
1606 TSRosWGetType - Get the type of Rosenbrock-W scheme
1607
1608 Logically Collective
1609
1610 Input Parameter:
1611 . ts - timestepping context
1612
1613 Output Parameter:
1614 . rostype - type of Rosenbrock-W scheme
1615
1616 Level: intermediate
1617
1618 .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()`
1619 @*/
TSRosWGetType(TS ts,TSRosWType * rostype)1620 PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype)
1621 {
1622 PetscFunctionBegin;
1623 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
1624 PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype));
1625 PetscFunctionReturn(PETSC_SUCCESS);
1626 }
1627
1628 /*@
1629 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step.
1630
1631 Logically Collective
1632
1633 Input Parameters:
1634 + ts - timestepping context
1635 - flg - `PETSC_TRUE` to recompute the Jacobian at each stage
1636
1637 Level: intermediate
1638
1639 .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()`
1640 @*/
TSRosWSetRecomputeJacobian(TS ts,PetscBool flg)1641 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg)
1642 {
1643 PetscFunctionBegin;
1644 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
1645 PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg));
1646 PetscFunctionReturn(PETSC_SUCCESS);
1647 }
1648
TSRosWGetType_RosW(TS ts,TSRosWType * rostype)1649 static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype)
1650 {
1651 TS_RosW *ros = (TS_RosW *)ts->data;
1652
1653 PetscFunctionBegin;
1654 *rostype = ros->tableau->name;
1655 PetscFunctionReturn(PETSC_SUCCESS);
1656 }
1657
TSRosWSetType_RosW(TS ts,TSRosWType rostype)1658 static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype)
1659 {
1660 TS_RosW *ros = (TS_RosW *)ts->data;
1661 PetscBool match;
1662 RosWTableauLink link;
1663
1664 PetscFunctionBegin;
1665 if (ros->tableau) {
1666 PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match));
1667 if (match) PetscFunctionReturn(PETSC_SUCCESS);
1668 }
1669 for (link = RosWTableauList; link; link = link->next) {
1670 PetscCall(PetscStrcmp(link->tab.name, rostype, &match));
1671 if (match) {
1672 if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts));
1673 ros->tableau = &link->tab;
1674 if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts));
1675 ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1676 PetscFunctionReturn(PETSC_SUCCESS);
1677 }
1678 }
1679 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype);
1680 }
1681
TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg)1682 static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg)
1683 {
1684 TS_RosW *ros = (TS_RosW *)ts->data;
1685
1686 PetscFunctionBegin;
1687 ros->recompute_jacobian = flg;
1688 PetscFunctionReturn(PETSC_SUCCESS);
1689 }
1690
TSDestroy_RosW(TS ts)1691 static PetscErrorCode TSDestroy_RosW(TS ts)
1692 {
1693 PetscFunctionBegin;
1694 PetscCall(TSReset_RosW(ts));
1695 if (ts->dm) {
1696 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts));
1697 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts));
1698 }
1699 PetscCall(PetscFree(ts->data));
1700 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL));
1701 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL));
1702 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL));
1703 PetscFunctionReturn(PETSC_SUCCESS);
1704 }
1705
1706 /*MC
1707 TSROSW - ODE solver using Rosenbrock-W schemes
1708
1709 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1710 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1711 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`.
1712
1713 Level: beginner
1714
1715 Notes:
1716 This is an IMEX method.
1717
1718 This method currently only works with autonomous ODE and DAE.
1719
1720 Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear.
1721
1722 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`.
1723
1724 Developer Notes:
1725 Rosenbrock-W methods are typically specified for autonomous ODE
1726 $$
1727 \dot{u} = f(u)
1728 $$
1729 by the stage equations
1730 $$
1731 k_i = h f(u_0 + \sum_j \alpha_{ij} k_j) + h J \sum_j \gamma_{ij} k_j
1732 $$
1733 and step completion formula
1734 $$
1735 u_1 = u_0 + \sum_j b_j k_j
1736 $$
1737 with step size $h$ and coefficients $\alpha_{ij}$, $\gamma_{ij}$, and $b_i$. Implementing the method in this form would require $f(u)$
1738 and the Jacobian $J$ to be available, in addition to the shifted matrix $I - h \gamma_{ii} J$. Following Hairer and Wanner,
1739 we define new variables for the stage equations
1740 $$
1741 y_i = \gamma_{ij} k_j
1742 $$
1743 The $ k_j $ can be recovered because $\Gamma$ is invertible. Let $C$ be the lower triangular part of $\Gamma^{-1}$ and define
1744 $$
1745 A = \Alpha \Gamma^{-1}, bt^T = b^T \Gamma^{-1}
1746 $$
1747 to rewrite the method as
1748 $$
1749 [M/(h \gamma_ii) - J] y_i = f(u_0 + \sum_j a_{ij} y_j) + M \sum_j (c_{ij}/h) y_j \\
1750 u_1 = u_0 + \sum_j bt_j y_j
1751 $$
1752
1753 where we have introduced the mass matrix $M$. Continue by defining
1754 $$
1755 \dot{y}_i = 1/(h \gamma_ii) y_i - \sum_j (c_{ij}/h) y_j
1756 $$
1757 or, more compactly in tensor notation
1758 $$
1759 \dot{Y} = 1/h (Gamma^{-1} \otimes I) Y .
1760 $$
1761 Note that $\Gamma^{-1}$ is lower triangular. With this definition of $\dot{Y} in terms of known quantities and the current
1762 stage $y_i$, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the
1763 equation
1764 $$
1765 g(u_0 + \sum_j a_{ij} y_j + y_i, \dot{y}_i) = 0
1766 $$
1767 with initial guess $y_i = 0$.
1768
1769 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`,
1770 `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`,
1771 `TSROSW4L`, `TSType`
1772 M*/
TSCreate_RosW(TS ts)1773 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts)
1774 {
1775 TS_RosW *ros;
1776
1777 PetscFunctionBegin;
1778 PetscCall(TSRosWInitializePackage());
1779
1780 ts->ops->reset = TSReset_RosW;
1781 ts->ops->destroy = TSDestroy_RosW;
1782 ts->ops->view = TSView_RosW;
1783 ts->ops->load = TSLoad_RosW;
1784 ts->ops->setup = TSSetUp_RosW;
1785 ts->ops->step = TSStep_RosW;
1786 ts->ops->interpolate = TSInterpolate_RosW;
1787 ts->ops->evaluatestep = TSEvaluateStep_RosW;
1788 ts->ops->setfromoptions = TSSetFromOptions_RosW;
1789 ts->ops->snesfunction = SNESTSFormFunction_RosW;
1790 ts->ops->snesjacobian = SNESTSFormJacobian_RosW;
1791
1792 ts->usessnes = PETSC_TRUE;
1793
1794 PetscCall(PetscNew(&ros));
1795 ts->data = (void *)ros;
1796
1797 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW));
1798 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW));
1799 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW));
1800
1801 PetscCall(TSRosWSetType(ts, TSRosWDefault));
1802 PetscFunctionReturn(PETSC_SUCCESS);
1803 }
1804