Actual source code: rosw.c
petsc-3.3-p7 2013-05-11
1: /*
2: Code for timestepping with Rosenbrock W methods
4: Notes:
5: The general system is written as
7: G(t,X,Xdot) = F(t,X)
9: where G represents the stiff part of the physics and F represents the non-stiff part.
10: This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.
12: */
13: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
15: #include <../src/mat/blockinvert.h>
17: static const TSRosWType TSRosWDefault = TSROSWRA34PW2;
18: static PetscBool TSRosWRegisterAllCalled;
19: static PetscBool TSRosWPackageInitialized;
21: typedef struct _RosWTableau *RosWTableau;
22: struct _RosWTableau {
23: char *name;
24: PetscInt order; /* Classical approximation order of the method */
25: PetscInt s; /* Number of stages */
26: PetscInt pinterp; /* Interpolation order */
27: PetscReal *A; /* Propagation table, strictly lower triangular */
28: PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */
29: PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */
30: PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/
31: PetscReal *b; /* Step completion table */
32: PetscReal *bembed; /* Step completion table for embedded method of order one less */
33: PetscReal *ASum; /* Row sum of A */
34: PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */
35: PetscReal *At; /* Propagation table in transformed variables */
36: PetscReal *bt; /* Step completion table in transformed variables */
37: PetscReal *bembedt; /* Step completion table of order one less in transformed variables */
38: PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */
39: PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
40: PetscReal *binterpt; /* Dense output formula */
41: };
42: typedef struct _RosWTableauLink *RosWTableauLink;
43: struct _RosWTableauLink {
44: struct _RosWTableau tab;
45: RosWTableauLink next;
46: };
47: static RosWTableauLink RosWTableauList;
49: typedef struct {
50: RosWTableau tableau;
51: Vec *Y; /* States computed during the step, used to complete the step */
52: Vec Ydot; /* Work vector holding Ydot during residual evaluation */
53: Vec Ystage; /* Work vector for the state value at each stage */
54: Vec Zdot; /* Ydot = Zdot + shift*Y */
55: Vec Zstage; /* Y = Zstage + Y */
56: Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/
57: PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */
58: PetscReal shift;
59: PetscReal stage_time;
60: PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */
61: PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
62: TSStepStatus status;
63: } TS_RosW;
65: /*MC
66: TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method).
68: Only an approximate Jacobian is needed.
70: Level: intermediate
72: .seealso: TSROSW
73: M*/
75: /*MC
76: TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method).
78: Only an approximate Jacobian is needed.
80: Level: intermediate
82: .seealso: TSROSW
83: M*/
85: /*MC
86: TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme.
88: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P.
90: Level: intermediate
92: .seealso: TSROSW
93: M*/
95: /*MC
96: TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme.
98: Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M.
100: Level: intermediate
102: .seealso: TSROSW
103: M*/
105: /*MC
106: TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1.
108: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
110: 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: References:
113: Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005.
115: Level: intermediate
117: .seealso: TSROSW
118: M*/
120: /*MC
121: TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1.
123: Only an approximate Jacobian is needed. By default, it is only recomputed once per step.
125: 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: References:
128: Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005.
130: Level: intermediate
132: .seealso: TSROSW
133: M*/
135: /*MC
136: TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme
138: By default, the Jacobian is only recomputed once per step.
140: Both the third order and embedded second order methods are stiffly accurate and L-stable.
142: References:
143: Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.
145: Level: intermediate
147: .seealso: TSROSW, TSROSWSANDU3
148: M*/
150: /*MC
151: TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme
153: By default, the Jacobian is only recomputed once per step.
155: The third order method is L-stable, but not stiffly accurate.
156: The second order embedded method is strongly A-stable with R(infty) = 0.5.
157: The internal stages are L-stable.
158: This method is called ROS3 in the paper.
160: References:
161: Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997.
163: Level: intermediate
165: .seealso: TSROSW, TSROSWRODAS3
166: M*/
168: /*MC
169: TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages
171: By default, the Jacobian is only recomputed once per step.
173: A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3)
175: References:
176: Emil Constantinescu
178: Level: intermediate
180: .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP
181: M*/
183: /*MC
184: TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
186: By default, the Jacobian is only recomputed once per step.
188: L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
190: References:
191: Emil Constantinescu
193: Level: intermediate
195: .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP
196: M*/
198: /*MC
199: TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages
201: By default, the Jacobian is only recomputed once per step.
203: L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2)
205: References:
206: Emil Constantinescu
208: Level: intermediate
210: .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP
211: M*/
215: /*@C
216: TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW
218: Not Collective, but should be called by all processes which will need the schemes to be registered
220: Level: advanced
222: .keywords: TS, TSRosW, register, all
224: .seealso: TSRosWRegisterDestroy()
225: @*/
226: PetscErrorCode TSRosWRegisterAll(void)
227: {
231: if (TSRosWRegisterAllCalled) return(0);
232: TSRosWRegisterAllCalled = PETSC_TRUE;
234: {
235: const PetscReal
236: A = 0,
237: Gamma = 1,
238: b = 1,
239: binterpt=1;
241: TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);
242: }
244: {
245: const PetscReal
246: A= 0,
247: Gamma = 0.5,
248: b = 1,
249: binterpt=1;
250: TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);
251: }
253: {
254: const PetscReal g = 1. + 1./PetscSqrtReal(2.0);
255: const PetscReal
256: A[2][2] = {{0,0}, {1.,0}},
257: Gamma[2][2] = {{g,0}, {-2.*g,g}},
258: b[2] = {0.5,0.5},
259: b1[2] = {1.0,0.0};
260: PetscReal binterpt[2][2];
261: binterpt[0][0]=g-1.0;
262: binterpt[1][0]=2.0-g;
263: binterpt[0][1]=g-1.5;
264: binterpt[1][1]=1.5-g;
265: TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);
266: }
267: {
268: const PetscReal g = 1. - 1./PetscSqrtReal(2.0);
269: const PetscReal
270: A[2][2] = {{0,0}, {1.,0}},
271: Gamma[2][2] = {{g,0}, {-2.*g,g}},
272: b[2] = {0.5,0.5},
273: b1[2] = {1.0,0.0};
274: PetscReal binterpt[2][2];
275: binterpt[0][0]=g-1.0;
276: binterpt[1][0]=2.0-g;
277: binterpt[0][1]=g-1.5;
278: binterpt[1][1]=1.5-g;
279: TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);
280: }
281: {
282: const PetscReal g = 7.8867513459481287e-01;
283: PetscReal binterpt[3][2];
284: const PetscReal
285: A[3][3] = {{0,0,0},
286: {1.5773502691896257e+00,0,0},
287: {0.5,0,0}},
288: Gamma[3][3] = {{g,0,0},
289: {-1.5773502691896257e+00,g,0},
290: {-6.7075317547305480e-01,-1.7075317547305482e-01,g}},
291: b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01},
292: b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01};
294: binterpt[0][0]=-0.8094010767585034;
295: binterpt[1][0]=-0.5;
296: binterpt[2][0]=2.3094010767585034;
297: binterpt[0][1]=0.9641016151377548;
298: binterpt[1][1]=0.5;
299: binterpt[2][1]=-1.4641016151377548;
300: TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);
301: }
302: {
303: PetscReal binterpt[4][3];
304: const PetscReal g = 4.3586652150845900e-01;
305: const PetscReal
306: A[4][4] = {{0,0,0,0},
307: {8.7173304301691801e-01,0,0,0},
308: {8.4457060015369423e-01,-1.1299064236484185e-01,0,0},
309: {0,0,1.,0}},
310: Gamma[4][4] = {{g,0,0,0},
311: {-8.7173304301691801e-01,g,0,0},
312: {-9.0338057013044082e-01,5.4180672388095326e-02,g,0},
313: {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,g}},
314: b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01},
315: b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01};
317: binterpt[0][0]=1.0564298455794094;
318: binterpt[1][0]=2.296429974281067;
319: binterpt[2][0]=-1.307599564525376;
320: binterpt[3][0]=-1.045260255335102;
321: binterpt[0][1]=-1.3864882699759573;
322: binterpt[1][1]=-8.262611700275677;
323: binterpt[2][1]=7.250979895056055;
324: binterpt[3][1]=2.398120075195581;
325: binterpt[0][2]=0.5721822314575016;
326: binterpt[1][2]=4.742931142090097;
327: binterpt[2][2]=-4.398120075195578;
328: binterpt[3][2]=-0.9169932983520199;
330: TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
331: }
332: {
333: const PetscReal g = 0.5;
334: const PetscReal
335: A[4][4] = {{0,0,0,0},
336: {0,0,0,0},
337: {1.,0,0,0},
338: {0.75,-0.25,0.5,0}},
339: Gamma[4][4] = {{g,0,0,0},
340: {1.,g,0,0},
341: {-0.25,-0.25,g,0},
342: {1./12,1./12,-2./3,g}},
343: b[4] = {5./6,-1./6,-1./6,0.5},
344: b2[4] = {0.75,-0.25,0.5,0};
345: TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,PETSC_NULL);
346: }
347: {
348: const PetscReal g = 0.43586652150845899941601945119356;
349: const PetscReal
350: A[3][3] = {{0,0,0},
351: {g,0,0},
352: {g,0,0}},
353: Gamma[3][3] = {{g,0,0},
354: {-0.19294655696029095575009695436041,g,0},
355: {0,1.74927148125794685173529749738960,g}},
356: b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829},
357: b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619};
359: PetscReal binterpt[3][2];
360: binterpt[0][0]=3.793692883777660870425141387941;
361: binterpt[1][0]=-2.918692883777660870425141387941;
362: binterpt[2][0]=0.125;
363: binterpt[0][1]=-0.725741064379812106687651020584;
364: binterpt[1][1]=0.559074397713145440020984353917;
365: binterpt[2][1]=0.16666666666666666666666666666667;
367: TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);
368: }
369: {
370: const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0;
371: const PetscReal
372: A[3][3] = {{0,0,0},
373: {1,0,0},
374: {0.25,0.25,0}},
375: Gamma[3][3] = {{0,0,0},
376: {(-3.0-s3)/6.0,g,0},
377: {(-3.0-s3)/24.0,(-3.0-s3)/8.0,g}},
378: b[3] = {1./6.,1./6.,2./3.},
379: b2[3] = {1./4.,1./4.,1./2.};
381: PetscReal binterpt[3][2];
382: binterpt[0][0]=0.089316397477040902157517886164709;
383: binterpt[1][0]=-0.91068360252295909784248211383529;
384: binterpt[2][0]=1.8213672050459181956849642276706;
385: binterpt[0][1]=0.077350269189625764509148780501957;
386: binterpt[1][1]=1.077350269189625764509148780502;
387: binterpt[2][1]=-1.1547005383792515290182975610039;
388: TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);
389: }
391: {
392: const PetscReal
393: A[4][4] = {{0,0,0,0},
394: {1./2.,0,0,0},
395: {1./2.,1./2.,0,0},
396: {1./6.,1./6.,1./6.,0}},
397: Gamma[4][4] = {{1./2.,0,0,0},
398: {0.0,1./4.,0,0},
399: {-2.,-2./3.,2./3.,0},
400: {1./2.,5./36.,-2./9,0}},
401: b[4] = {1./6.,1./6.,1./6.,1./2.},
402: b2[4] = {1./8.,3./4.,1./8.,0};
403: PetscReal binterpt[4][3];
404: binterpt[0][0]=6.25;
405: binterpt[1][0]=-30.25;
406: binterpt[2][0]=1.75;
407: binterpt[3][0]=23.25;
408: binterpt[0][1]=-9.75;
409: binterpt[1][1]=58.75;
410: binterpt[2][1]=-3.25;
411: binterpt[3][1]=-45.75;
412: binterpt[0][2]=3.6666666666666666666666666666667;
413: binterpt[1][2]=-28.333333333333333333333333333333;
414: binterpt[2][2]=1.6666666666666666666666666666667;
415: binterpt[3][2]=23.;
416: TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
417: }
419: {
420: const PetscReal
421: A[4][4] = {{0,0,0,0},
422: {1./2.,0,0,0},
423: {1./2.,1./2.,0,0},
424: {1./6.,1./6.,1./6.,0}},
425: Gamma[4][4] = {{1./2.,0,0,0},
426: {0.0,3./4.,0,0},
427: {-2./3.,-23./9.,2./9.,0},
428: {1./18.,65./108.,-2./27,0}},
429: b[4] = {1./6.,1./6.,1./6.,1./2.},
430: b2[4] = {3./16.,10./16.,3./16.,0};
432: PetscReal binterpt[4][3];
433: binterpt[0][0]=1.6911764705882352941176470588235;
434: binterpt[1][0]=3.6813725490196078431372549019608;
435: binterpt[2][0]=0.23039215686274509803921568627451;
436: binterpt[3][0]=-4.6029411764705882352941176470588;
437: binterpt[0][1]=-0.95588235294117647058823529411765;
438: binterpt[1][1]=-6.2401960784313725490196078431373;
439: binterpt[2][1]=-0.31862745098039215686274509803922;
440: binterpt[3][1]=7.5147058823529411764705882352941;
441: binterpt[0][2]=-0.56862745098039215686274509803922;
442: binterpt[1][2]=2.7254901960784313725490196078431;
443: binterpt[2][2]=0.25490196078431372549019607843137;
444: binterpt[3][2]=-2.4117647058823529411764705882353;
445: TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
446: }
448: {
449: PetscReal A[4][4],Gamma[4][4],b[4],b2[4];
450: PetscReal binterpt[4][3];
452: Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816;
453: Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0;
454: Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476;
455: Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816;
456: Gamma[1][2]=0; Gamma[1][3]=0;
457: Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903;
458: Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131;
459: Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816;
460: Gamma[2][3]=0;
461: Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783;
462: Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984;
463: Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198;
464: Gamma[3][3]=0;
466: A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0;
467: A[1][0]=0.8717330430169179988320388950590125027645343373957631;
468: A[1][1]=0; A[1][2]=0; A[1][3]=0;
469: A[2][0]=0.5275890119763004115618079766722914408876108660811028;
470: A[2][1]=0.07241098802369958843819203208518599088698057726988732;
471: A[2][2]=0; A[2][3]=0;
472: A[3][0]=0.3990960076760701320627260685975778145384666450351314;
473: A[3][1]=-0.4375576546135194437228463747348862825846903771419953;
474: A[3][2]=1.038461646937449311660120300601880176655352737312713;
475: A[3][3]=0;
477: b[0]=0.1876410243467238251612921333138006734899663569186926;
478: b[1]=-0.5952974735769549480478230473706443582188442040780541;
479: b[2]=0.9717899277217721234705114616271378792182450260943198;
480: b[3]=0.4358665215084589994160194475295062513822671686978816;
482: b2[0]=0.2147402862233891404862383521089097657790734483804460;
483: b2[1]=-0.4851622638849390928209050538171743017757490232519684;
484: b2[2]=0.8687250025203875511662123688667549217531982787600080;
485: b2[3]=0.4016969751411624011684543450940068201770721128357014;
487: binterpt[0][0]=2.2565812720167954547104627844105;
488: binterpt[1][0]=1.349166413351089573796243820819;
489: binterpt[2][0]=-2.4695174540533503758652847586647;
490: binterpt[3][0]=-0.13623023131453465264142184656474;
491: binterpt[0][1]=-3.0826699111559187902922463354557;
492: binterpt[1][1]=-2.4689115685996042534544925650515;
493: binterpt[2][1]=5.7428279814696677152129332773553;
494: binterpt[3][1]=-0.19124650171414467146619437684812;
495: binterpt[0][2]=1.0137296634858471607430756831148;
496: binterpt[1][2]=0.52444768167155973161042570784064;
497: binterpt[2][2]=-2.3015205996945452158771370439586;
498: binterpt[3][2]=0.76334325453713832352363565300308;
500: TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);
501: }
503: return(0);
504: }
508: /*@C
509: TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister().
511: Not Collective
513: Level: advanced
515: .keywords: TSRosW, register, destroy
516: .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic()
517: @*/
518: PetscErrorCode TSRosWRegisterDestroy(void)
519: {
521: RosWTableauLink link;
524: while ((link = RosWTableauList)) {
525: RosWTableau t = &link->tab;
526: RosWTableauList = link->next;
527: PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);
528: PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);
529: PetscFree2(t->bembed,t->bembedt);
530: PetscFree(t->binterpt);
531: PetscFree(t->name);
532: PetscFree(link);
533: }
534: TSRosWRegisterAllCalled = PETSC_FALSE;
535: return(0);
536: }
540: /*@C
541: TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called
542: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW()
543: when using static libraries.
545: Input Parameter:
546: path - The dynamic library path, or PETSC_NULL
548: Level: developer
550: .keywords: TS, TSRosW, initialize, package
551: .seealso: PetscInitialize()
552: @*/
553: PetscErrorCode TSRosWInitializePackage(const char path[])
554: {
558: if (TSRosWPackageInitialized) return(0);
559: TSRosWPackageInitialized = PETSC_TRUE;
560: TSRosWRegisterAll();
561: PetscRegisterFinalize(TSRosWFinalizePackage);
562: return(0);
563: }
567: /*@C
568: TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is
569: called from PetscFinalize().
571: Level: developer
573: .keywords: Petsc, destroy, package
574: .seealso: PetscFinalize()
575: @*/
576: PetscErrorCode TSRosWFinalizePackage(void)
577: {
581: TSRosWPackageInitialized = PETSC_FALSE;
582: TSRosWRegisterDestroy();
583: return(0);
584: }
588: /*@C
589: TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
591: Not Collective, but the same schemes should be registered on all processes on which they will be used
593: Input Parameters:
594: + name - identifier for method
595: . order - approximation order of method
596: . s - number of stages, this is the dimension of the matrices below
597: . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular
598: . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal
599: . b - Step completion table (dimension s)
600: - bembed - Step completion table for a scheme of order one less (dimension s, PETSC_NULL if no embedded scheme is available)
601: . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt
602: . binterpt - Coefficients of the interpolation formula (dimension s*pinterp)
604: Notes:
605: Several Rosenbrock W methods are provided, this function is only needed to create new methods.
607: Level: advanced
609: .keywords: TS, register
611: .seealso: TSRosW
612: @*/
613: PetscErrorCode TSRosWRegister(const TSRosWType name,PetscInt order,PetscInt s,
614: const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[],
615: PetscInt pinterp,const PetscReal binterpt[])
616: {
618: RosWTableauLink link;
619: RosWTableau t;
620: PetscInt i,j,k;
621: PetscScalar *GammaInv;
630: PetscMalloc(sizeof(*link),&link);
631: PetscMemzero(link,sizeof(*link));
632: t = &link->tab;
633: PetscStrallocpy(name,&t->name);
634: t->order = order;
635: t->s = s;
636: PetscMalloc5(s*s,PetscReal,&t->A,s*s,PetscReal,&t->Gamma,s,PetscReal,&t->b,s,PetscReal,&t->ASum,s,PetscReal,&t->GammaSum);
637: PetscMalloc5(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s*s,PetscReal,&t->GammaInv,s,PetscBool,&t->GammaZeroDiag,s*s,PetscReal,&t->GammaExplicitCorr);
638: PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
639: PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));
640: PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));
641: PetscMemcpy(t->b,b,s*sizeof(b[0]));
642: if (bembed) {
643: PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);
644: PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
645: }
646: for (i=0; i<s; i++) {
647: t->ASum[i] = 0;
648: t->GammaSum[i] = 0;
649: for (j=0; j<s; j++) {
650: t->ASum[i] += A[i*s+j];
651: t->GammaSum[i] += Gamma[i*s+j];
652: }
653: }
654: PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv); /* Need to use Scalar for inverse, then convert back to Real */
655: for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i];
656: for (i=0; i<s; i++) {
657: if (Gamma[i*s+i] == 0.0) {
658: GammaInv[i*s+i] = 1.0;
659: t->GammaZeroDiag[i] = PETSC_TRUE;
660: } else {
661: t->GammaZeroDiag[i] = PETSC_FALSE;
662: }
663: }
665: switch (s) {
666: case 1: GammaInv[0] = 1./GammaInv[0]; break;
667: case 2: PetscKernel_A_gets_inverse_A_2(GammaInv,0); break;
668: case 3: PetscKernel_A_gets_inverse_A_3(GammaInv,0); break;
669: case 4: PetscKernel_A_gets_inverse_A_4(GammaInv,0); break;
670: case 5: {
671: PetscInt ipvt5[5];
672: MatScalar work5[5*5];
673: PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0); break;
674: }
675: case 6: PetscKernel_A_gets_inverse_A_6(GammaInv,0); break;
676: case 7: PetscKernel_A_gets_inverse_A_7(GammaInv,0); break;
677: default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s);
678: }
679: for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]);
680: PetscFree(GammaInv);
682: for (i=0; i<s; i++) {
683: for (k=0; k<i+1; k++) {
684: t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]);
685: for (j=k+1; j<i+1; j++) {
686: t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]);
687: }
688: }
689: }
691: for (i=0; i<s; i++) {
692: for (j=0; j<s; j++) {
693: t->At[i*s+j] = 0;
694: for (k=0; k<s; k++) {
695: t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j];
696: }
697: }
698: t->bt[i] = 0;
699: for (j=0; j<s; j++) {
700: t->bt[i] += t->b[j] * t->GammaInv[j*s+i];
701: }
702: if (bembed) {
703: t->bembedt[i] = 0;
704: for (j=0; j<s; j++) {
705: t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i];
706: }
707: }
708: }
709: t->ccfl = 1.0; /* Fix this */
711: t->pinterp = pinterp;
712: PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);
713: PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));
714: link->next = RosWTableauList;
715: RosWTableauList = link;
716: return(0);
717: }
721: /*
722: The step completion formula is
724: x1 = x0 + b^T Y
726: where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been
727: updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write
729: x1e = x0 + be^T Y
730: = x1 - b^T Y + be^T Y
731: = x1 + (be - b)^T Y
733: so we can evaluate the method of different order even after the step has been optimistically completed.
734: */
735: static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec X,PetscBool *done)
736: {
737: TS_RosW *ros = (TS_RosW*)ts->data;
738: RosWTableau tab = ros->tableau;
739: PetscScalar *w = ros->work;
740: PetscInt i;
744: if (order == tab->order) {
745: if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */
746: VecCopy(ts->vec_sol,X);
747: for (i=0; i<tab->s; i++) w[i] = tab->bt[i];
748: VecMAXPY(X,tab->s,w,ros->Y);
749: } else {VecCopy(ts->vec_sol,X);}
750: if (done) *done = PETSC_TRUE;
751: return(0);
752: } else if (order == tab->order-1) {
753: if (!tab->bembedt) goto unavailable;
754: if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */
755: VecCopy(ts->vec_sol,X);
756: for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i];
757: VecMAXPY(X,tab->s,w,ros->Y);
758: } else { /* Use rollback-and-recomplete formula (bembedt - bt) */
759: for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i];
760: VecCopy(ts->vec_sol,X);
761: VecMAXPY(X,tab->s,w,ros->Y);
762: }
763: if (done) *done = PETSC_TRUE;
764: return(0);
765: }
766: unavailable:
767: if (done) *done = PETSC_FALSE;
768: else SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
769: return(0);
770: }
774: static PetscErrorCode TSStep_RosW(TS ts)
775: {
776: TS_RosW *ros = (TS_RosW*)ts->data;
777: RosWTableau tab = ros->tableau;
778: const PetscInt s = tab->s;
779: const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv;
780: const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr;
781: const PetscBool *GammaZeroDiag = tab->GammaZeroDiag;
782: PetscScalar *w = ros->work;
783: Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage;
784: SNES snes;
785: TSAdapt adapt;
786: PetscInt i,j,its,lits,reject,next_scheme;
787: PetscReal next_time_step;
788: PetscBool accept;
789: PetscErrorCode ierr;
790: MatStructure str;
793: TSGetSNES(ts,&snes);
794: next_time_step = ts->time_step;
795: accept = PETSC_TRUE;
796: ros->status = TS_STEP_INCOMPLETE;
798: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
799: const PetscReal h = ts->time_step;
800: TSPreStep(ts);
801: VecCopy(ts->vec_sol,ros->VecSolPrev);/*move this at the end*/
802: for (i=0; i<s; i++) {
803: ros->stage_time = ts->ptime + h*ASum[i];
804: TSPreStage(ts,ros->stage_time);
805: if (GammaZeroDiag[i]) {
806: ros->stage_explicit = PETSC_TRUE;
807: ros->shift = 1./h;
808: } else {
809: ros->stage_explicit = PETSC_FALSE;
810: ros->shift = 1./(h*Gamma[i*s+i]);
811: }
813: VecCopy(ts->vec_sol,Zstage);
814: for (j=0; j<i; j++) w[j] = At[i*s+j];
815: VecMAXPY(Zstage,i,w,Y);
817: for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j];
818: VecZeroEntries(Zdot);
819: VecMAXPY(Zdot,i,w,Y);
821: /* Initial guess taken from last stage */
822: VecZeroEntries(Y[i]);
824: if (!ros->stage_explicit) {
825: if (!ros->recompute_jacobian && !i) {
826: SNESSetLagJacobian(snes,-2); /* Recompute the Jacobian on this solve, but not again */
827: }
828: SNESSolve(snes,PETSC_NULL,Y[i]);
829: SNESGetIterationNumber(snes,&its);
830: SNESGetLinearSolveIterations(snes,&lits);
831: ts->snes_its += its; ts->ksp_its += lits;
832: TSGetAdapt(ts,&adapt);
833: TSAdaptCheckStage(adapt,ts,&accept);
834: if (!accept) goto reject_step;
835: } else {
836: Mat J,Jp;
837: VecZeroEntries(Ydot); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */
838: TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);
839: VecScale(Y[i],-1.0);
840: VecAXPY(Y[i],-1.0,Zdot); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/
841:
842: VecZeroEntries(Zstage); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */
843: for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j];
844: VecMAXPY(Zstage,i,w,Y);
845: /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */
846: str = SAME_NONZERO_PATTERN;
847: TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);
848: TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);
849: MatMult(J,Zstage,Zdot);
851: VecAXPY(Y[i],-1.0,Zdot);
852: VecScale(Y[i],h);
853: ts->ksp_its += 1;
854: }
855: }
856: TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);
857: ros->status = TS_STEP_PENDING;
859: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
860: TSGetAdapt(ts,&adapt);
861: TSAdaptCandidatesClear(adapt);
862: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
863: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
864: if (accept) {
865: /* ignore next_scheme for now */
866: ts->ptime += ts->time_step;
867: ts->time_step = next_time_step;
868: ts->steps++;
869: ros->status = TS_STEP_COMPLETE;
870: break;
871: } else { /* Roll back the current step */
872: for (i=0; i<s; i++) w[i] = -tab->bt[i];
873: VecMAXPY(ts->vec_sol,s,w,Y);
874: ts->time_step = next_time_step;
875: ros->status = TS_STEP_INCOMPLETE;
876: }
877: reject_step: continue;
878: }
879: if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
880: return(0);
881: }
885: static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X)
886: {
887: TS_RosW *ros = (TS_RosW*)ts->data;
888: PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j;
889: PetscReal h;
890: PetscReal tt,t;
891: PetscScalar *bt;
892: const PetscReal *Bt = ros->tableau->binterpt;
894: const PetscReal *GammaInv = ros->tableau->GammaInv;
895: PetscScalar *w = ros->work;
896: Vec *Y = ros->Y;
899: if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name);
901: switch (ros->status) {
902: case TS_STEP_INCOMPLETE:
903: case TS_STEP_PENDING:
904: h = ts->time_step;
905: t = (itime - ts->ptime)/h;
906: break;
907: case TS_STEP_COMPLETE:
908: h = ts->time_step_prev;
909: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
910: break;
911: default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus");
912: }
913: PetscMalloc(s*sizeof(bt[0]),&bt);
914: for (i=0; i<s; i++) bt[i] = 0;
915: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
916: for (i=0; i<s; i++) {
917: bt[i] += Bt[i*pinterp+j] * tt;
918: }
919: }
921: /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */
922: /*X<-0*/
923: VecZeroEntries(X);
924:
925: /*X<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */
926: for (j=0; j<s; j++) w[j]=0;
927: for (j=0; j<s; j++) {
928: for (i=j; i<s; i++) {
929: w[j] += bt[i]*GammaInv[i*s+j];
930: }
931: }
932: VecMAXPY(X,i,w,Y);
934: /*X<-y(t) + X*/
935: VecAXPY(X,1.0,ros->VecSolPrev);
936:
937: PetscFree(bt);
939: return(0);
940: }
942: /*------------------------------------------------------------*/
945: static PetscErrorCode TSReset_RosW(TS ts)
946: {
947: TS_RosW *ros = (TS_RosW*)ts->data;
948: PetscInt s;
952: if (!ros->tableau) return(0);
953: s = ros->tableau->s;
954: VecDestroyVecs(s,&ros->Y);
955: VecDestroy(&ros->Ydot);
956: VecDestroy(&ros->Ystage);
957: VecDestroy(&ros->Zdot);
958: VecDestroy(&ros->Zstage);
959: VecDestroy(&ros->VecSolPrev);
960: PetscFree(ros->work);
961: return(0);
962: }
966: static PetscErrorCode TSDestroy_RosW(TS ts)
967: {
968: PetscErrorCode ierr;
971: TSReset_RosW(ts);
972: PetscFree(ts->data);
973: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);
974: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);
975: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);
976: return(0);
977: }
979: /*
980: This defines the nonlinear equation that is to be solved with SNES
981: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
982: */
985: static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts)
986: {
987: TS_RosW *ros = (TS_RosW*)ts->data;
991: VecWAXPY(ros->Ydot,ros->shift,X,ros->Zdot); /* Ydot = shift*X + Zdot */
992: VecWAXPY(ros->Ystage,1.0,X,ros->Zstage); /* Ystage = X + Zstage */
993: TSComputeIFunction(ts,ros->stage_time,ros->Ystage,ros->Ydot,F,PETSC_FALSE);
994: return(0);
995: }
999: static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts)
1000: {
1001: TS_RosW *ros = (TS_RosW*)ts->data;
1005: /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */
1006: TSComputeIJacobian(ts,ros->stage_time,ros->Ystage,ros->Ydot,ros->shift,A,B,str,PETSC_TRUE);
1007: return(0);
1008: }
1012: static PetscErrorCode TSSetUp_RosW(TS ts)
1013: {
1014: TS_RosW *ros = (TS_RosW*)ts->data;
1015: RosWTableau tab = ros->tableau;
1016: PetscInt s = tab->s;
1020: if (!ros->tableau) {
1021: TSRosWSetType(ts,TSRosWDefault);
1022: }
1023: VecDuplicateVecs(ts->vec_sol,s,&ros->Y);
1024: VecDuplicate(ts->vec_sol,&ros->Ydot);
1025: VecDuplicate(ts->vec_sol,&ros->Ystage);
1026: VecDuplicate(ts->vec_sol,&ros->Zdot);
1027: VecDuplicate(ts->vec_sol,&ros->Zstage);
1028: VecDuplicate(ts->vec_sol,&ros->VecSolPrev);
1029: PetscMalloc(s*sizeof(ros->work[0]),&ros->work);
1030: return(0);
1031: }
1032: /*------------------------------------------------------------*/
1036: static PetscErrorCode TSSetFromOptions_RosW(TS ts)
1037: {
1038: TS_RosW *ros = (TS_RosW*)ts->data;
1040: char rostype[256];
1043: PetscOptionsHead("RosW ODE solver options");
1044: {
1045: RosWTableauLink link;
1046: PetscInt count,choice;
1047: PetscBool flg;
1048: const char **namelist;
1049: SNES snes;
1051: PetscStrncpy(rostype,TSRosWDefault,sizeof rostype);
1052: for (link=RosWTableauList,count=0; link; link=link->next,count++) ;
1053: PetscMalloc(count*sizeof(char*),&namelist);
1054: for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
1055: PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);
1056: TSRosWSetType(ts,flg ? namelist[choice] : rostype);
1057: PetscFree(namelist);
1059: PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);
1061: /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */
1062: TSGetSNES(ts,&snes);
1063: if (!((PetscObject)snes)->type_name) {
1064: SNESSetType(snes,SNESKSPONLY);
1065: }
1066: SNESSetFromOptions(snes);
1067: }
1068: PetscOptionsTail();
1069: return(0);
1070: }
1074: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
1075: {
1077: PetscInt i;
1078: size_t left,count;
1079: char *p;
1082: for (i=0,p=buf,left=len; i<n; i++) {
1083: PetscSNPrintfCount(p,left,fmt,&count,x[i]);
1084: if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
1085: left -= count;
1086: p += count;
1087: *p++ = ' ';
1088: }
1089: p[i ? 0 : -1] = 0;
1090: return(0);
1091: }
1095: static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer)
1096: {
1097: TS_RosW *ros = (TS_RosW*)ts->data;
1098: RosWTableau tab = ros->tableau;
1099: PetscBool iascii;
1103: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1104: if (iascii) {
1105: const TSRosWType rostype;
1106: PetscInt i;
1107: PetscReal abscissa[512];
1108: char buf[512];
1109: TSRosWGetType(ts,&rostype);
1110: PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);
1111: PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ASum);
1112: PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);
1113: for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i];
1114: PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,abscissa);
1115: PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);
1116: }
1117: SNESView(ts->snes,viewer);
1118: return(0);
1119: }
1123: /*@C
1124: TSRosWSetType - Set the type of Rosenbrock-W scheme
1126: Logically collective
1128: Input Parameter:
1129: + ts - timestepping context
1130: - rostype - type of Rosenbrock-W scheme
1132: Level: beginner
1134: .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3
1135: @*/
1136: PetscErrorCode TSRosWSetType(TS ts,const TSRosWType rostype)
1137: {
1142: PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,rostype));
1143: return(0);
1144: }
1148: /*@C
1149: TSRosWGetType - Get the type of Rosenbrock-W scheme
1151: Logically collective
1153: Input Parameter:
1154: . ts - timestepping context
1156: Output Parameter:
1157: . rostype - type of Rosenbrock-W scheme
1159: Level: intermediate
1161: .seealso: TSRosWGetType()
1162: @*/
1163: PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *rostype)
1164: {
1169: PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,rostype));
1170: return(0);
1171: }
1175: /*@C
1176: TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step.
1178: Logically collective
1180: Input Parameter:
1181: + ts - timestepping context
1182: - flg - PETSC_TRUE to recompute the Jacobian at each stage
1184: Level: intermediate
1186: .seealso: TSRosWGetType()
1187: @*/
1188: PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg)
1189: {
1194: PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));
1195: return(0);
1196: }
1198: EXTERN_C_BEGIN
1201: PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *rostype)
1202: {
1203: TS_RosW *ros = (TS_RosW*)ts->data;
1207: if (!ros->tableau) {TSRosWSetType(ts,TSRosWDefault);}
1208: *rostype = ros->tableau->name;
1209: return(0);
1210: }
1213: PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType rostype)
1214: {
1215: TS_RosW *ros = (TS_RosW*)ts->data;
1216: PetscErrorCode ierr;
1217: PetscBool match;
1218: RosWTableauLink link;
1221: if (ros->tableau) {
1222: PetscStrcmp(ros->tableau->name,rostype,&match);
1223: if (match) return(0);
1224: }
1225: for (link = RosWTableauList; link; link=link->next) {
1226: PetscStrcmp(link->tab.name,rostype,&match);
1227: if (match) {
1228: TSReset_RosW(ts);
1229: ros->tableau = &link->tab;
1230: return(0);
1231: }
1232: }
1233: SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype);
1234: return(0);
1235: }
1239: PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg)
1240: {
1241: TS_RosW *ros = (TS_RosW*)ts->data;
1244: ros->recompute_jacobian = flg;
1245: return(0);
1246: }
1247: EXTERN_C_END
1249: /* ------------------------------------------------------------ */
1250: /*MC
1251: TSROSW - ODE solver using Rosenbrock-W schemes
1253: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1254: nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1255: of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
1257: Notes:
1258: This method currently only works with autonomous ODE and DAE.
1260: Developer notes:
1261: Rosenbrock-W methods are typically specified for autonomous ODE
1263: $ xdot = f(x)
1265: by the stage equations
1267: $ k_i = h f(x_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j
1269: and step completion formula
1271: $ x_1 = x_0 + sum_j b_j k_j
1273: with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(x)
1274: and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner,
1275: we define new variables for the stage equations
1277: $ y_i = gamma_ij k_j
1279: The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define
1281: $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i}
1283: to rewrite the method as
1285: $ [M/(h gamma_ii) - J] y_i = f(x_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j
1286: $ x_1 = x_0 + sum_j bt_j y_j
1288: where we have introduced the mass matrix M. Continue by defining
1290: $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j
1292: or, more compactly in tensor notation
1294: $ Ydot = 1/h (Gamma^{-1} \otimes I) Y .
1296: Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current
1297: stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the
1298: equation
1300: $ g(x_0 + sum_j a_ij y_j + y_i, ydot_i) = 0
1302: with initial guess y_i = 0.
1304: Level: beginner
1306: .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister()
1308: M*/
1309: EXTERN_C_BEGIN
1312: PetscErrorCode TSCreate_RosW(TS ts)
1313: {
1314: TS_RosW *ros;
1318: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1319: TSRosWInitializePackage(PETSC_NULL);
1320: #endif
1322: ts->ops->reset = TSReset_RosW;
1323: ts->ops->destroy = TSDestroy_RosW;
1324: ts->ops->view = TSView_RosW;
1325: ts->ops->setup = TSSetUp_RosW;
1326: ts->ops->step = TSStep_RosW;
1327: ts->ops->interpolate = TSInterpolate_RosW;
1328: ts->ops->evaluatestep = TSEvaluateStep_RosW;
1329: ts->ops->setfromoptions = TSSetFromOptions_RosW;
1330: ts->ops->snesfunction = SNESTSFormFunction_RosW;
1331: ts->ops->snesjacobian = SNESTSFormJacobian_RosW;
1333: PetscNewLog(ts,TS_RosW,&ros);
1334: ts->data = (void*)ros;
1336: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);
1337: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);
1338: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);
1339: return(0);
1340: }
1341: EXTERN_C_END