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