Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk.c @ 4517

Last change on this file since 4517 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk
File size: 38.8 KB
Line 
1#define MAX(a,b) ( ((a) >= (b)) ?(a):(b)  )
2#define MIN(b,c) ( ((b) <  (c)) ?(b):(c)  )
3#define ABS(x)   ( ((x) >=  0 ) ?(x):(-x) )
4#define SQRT(d)  ( pow((d),0.5)  )
5#define SIGN(x,y)( ( (x*y) >= 0 ) ?(x):(-x) )/* Sign transfer function */
6#define MOD(A,B) (int)((A)%(B))
7
8/* ~~~> Numerical constants    */
9#define  ZERO     (KPP_REAL)0.0
10#define  ONE      (KPP_REAL)1.0
11
12/* ~~~>  Statistics on the work performed by the SDIRK method */
13#define  Nfun 1
14#define  Njac 2
15#define  Nstp 3
16#define  Nacc 4
17#define  Nrej 5
18#define  Ndec 6
19#define  Nsol 7
20#define  Nsng 8
21#define  Ntexit 1
22#define  Nhexit 2
23#define  Nhnew 3
24
25/*~~~>  SDIRK method coefficients, up to 5 stages    */
26#define Smax 5
27
28int S2A=1,
29        S2B=2,
30        S3A=3, 
31        S4A=4, 
32        S4B=5;
33
34int sdMethod, 
35        rkS;
36
37KPP_REAL rkGamma, 
38        rkA[Smax][Smax], 
39        rkB[Smax], 
40        rkELO, 
41        rkBhat[Smax],
42        rkC[Smax], 
43        rkD[Smax], 
44        rkE[Smax], 
45        rkTheta[Smax][Smax], 
46        rkAlpha[Smax][Smax];
47
48/*~~~> Function headers     */
49//void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[],
50//      int ISTATUS_U[], KPP_REAL RSTATUS_U[], int Ierr);
51void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT);
52int SDIRK(int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[], 
53        KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[], int ICNTRL[], 
54        KPP_REAL RSTATUS[], int ISTATUS[]);
55int SDIRK_Integrator(int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[], 
56        int Ierr, KPP_REAL Hstart, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Roundoff, 
57        KPP_REAL AbsTol[], KPP_REAL RelTol[], int ISTATUS[], KPP_REAL RSTATUS[],
58        int ITOL, int Max_no_steps, int StartNewton, KPP_REAL NewtonTol, 
59        KPP_REAL ThetaMin, KPP_REAL FacSafe, KPP_REAL FacMin, KPP_REAL FacMax,
60        KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int NewtonMaxit);
61void SDIRK_ErrorScale(int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[],
62        KPP_REAL Y[], KPP_REAL SCAL[]);
63KPP_REAL SDIRK_ErrorNorm(int N, KPP_REAL Y[], KPP_REAL SCAL[]);
64int SDIRK_ErrorMsg(int code, KPP_REAL T, KPP_REAL H, int Ierr);
65void SDIRK_PrepareMatrix(KPP_REAL H, KPP_REAL T, KPP_REAL Y[], KPP_REAL FJAC[], 
66        int SkipJac, int SkipLU, KPP_REAL E[], int IP[], int Reject, 
67        int ISING, int ISTATUS[]);
68void SDIRK_Solve(KPP_REAL H, int N, KPP_REAL E[], int IP[], int ISING, 
69        KPP_REAL RHS[], int ISTATUS[]);
70void Sdirk4a(void);
71void Sdirk4b(void);
72void Sdirk2a(void);
73void Sdirk2b(void);
74void Sdirk3a(void);
75void FUN_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL P[]);
76void JAC_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL JV[]);
77void Fun(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]);
78void Jac_SP(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]);
79void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY);
80void WSCAL(int N, KPP_REAL Alpha, KPP_REAL X[], int incX);
81KPP_REAL WLAMCH(char C);
82void WADD(int N, KPP_REAL Y[], KPP_REAL Z[], KPP_REAL TMP[]);
83void Set2Zero(int N, KPP_REAL Y[]);
84void KppSolve(KPP_REAL A[], KPP_REAL b[]);
85int KppDecomp(KPP_REAL A[]);
86void Update_SUN();
87void Update_RCONST();
88
89/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
90//void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[],
91//              int ISTATUS_U[], KPP_REAL RSTATUS_U[], int Ierr_U)
92void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT)
93/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
94{
95
96/* int Ntotal = 0; *//* Used for debug option below to print the number of steps */
97KPP_REAL RCNTRL[20], 
98        RSTATUS[20], 
99        T1, 
100        T2;
101       
102int ICNTRL[20], 
103        ISTATUS[20], 
104        Ierr;
105
106Ierr = 0;
107
108int i;
109for(i=0; i<20; i++) {
110        ICNTRL[i] = 0;
111        RCNTRL[i] = (KPP_REAL)0.0;
112        ISTATUS[i] = 0;
113        RSTATUS[i] = (KPP_REAL)0.0;
114} /* end for */
115
116/*~> fine-tune the integrator: */
117ICNTRL[1] = 0; /* 0 - vector tolerances, 1 - scalar tolerances */
118ICNTRL[5] = 0; /* starting values of N. iter.: interpolated 0), zero (1) */
119
120///* If optional parameters are given, and if they are >0,
121//      then they overwrite default settings. */
122//if(ICNTRL_U != NULL) { /* Check to see if ICNTRL_U is not NULL */
123//      for(i=0; i<20; i++) {
124//              if(ICNTRL_U[i] > 0) {
125//                      ICNTRL[i] = ICNTRL_U[i];
126//              } /* end if */
127//      } /* end for */
128//} /* end if */
129//
130//if(RCNTRL_U != NULL) { /* Check to see if RCNTRL_U is not NULL */
131//      for(i=0; i<20; i++) {
132//              if(RCNTRL_U[i] > 0) {
133//                      RCNTRL[i] = RCNTRL_U[i];
134//              } /* end if */
135//      } /* end for */
136//} /* end if */
137
138
139T1 = TIN;
140T2 = TOUT;
141Ierr = SDIRK( NVAR, T1, T2, VAR, RTOL, ATOL, RCNTRL, ICNTRL, RSTATUS,
142        ISTATUS);
143
144
145/*~~~> Debug option: print number of steps
146   Ntotal += ISTATUS[Nstp]; */
147
148if(Ierr < 0) {
149        printf("SDIRK: Unsuccessful exit at T=%f(Ierr=%d)", TIN, Ierr);
150}
151
152///*if optional parameters are given for output they to return information */
153//if(ISTATUS_U != NULL) { /* Check to see if ISTATUS_U is not NULL */
154//      for(i=0; i<20; i++) {
155//              ISTATUS_U[i] = ISTATUS[i];
156//      } /* end for */
157//} /* end if */
158//
159//if(RSTATUS_U != NULL) { /* Check to see if RSTATUS_U is not NULL */
160//      for(i=0; i<20; i++) {
161//              RSTATUS_U[i] = RSTATUS[i];
162//      } /* end for */
163//} /* end if */
164//
165//Ierr_U = Ierr;
166
167}
168
169/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
170int SDIRK(int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[],
171        KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[],
172        int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[])
173/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
174
175    Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
176    Runge-Kutta (SDIRK) method.
177
178    This implementation is based on the book and the code Sdirk4:
179
180    E. Hairer and G. Wanner
181    "Solving ODEs II. Stiff and differential-algebraic problems".
182    Springer series in computational mathematics, Springer-Verlag, 1996.
183    This code is based on the SDIRK4 routine in the above book.
184
185    Methods:
186        * Sdirk 2a, 2b: L-stable, 2 stages, order 2
187        * Sdirk 3a:     L-stable, 3 stages, order 2, adjoint-invariant
188        * Sdirk 4a, 4b: L-stable, 5 stages, order 4
189
190    (C)  Adrian Sandu, July 2005
191    Virginia Polytechnic Institute and State University
192    Contact: sandu@cs.vt.edu
193    Revised by Philipp Miehe and Adrian Sandu, May 2006
194    Translation F90 to C by Paul Eller and Nicholas Hobbs, July 2006
195    This implementation is part of KPP - the Kinetic PreProcessor
196~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
197
198~~~>  INPUT ARGUMENTS:
199
200   -     Y[NVAR]    = vector of initial conditions (at T=Tinitial)
201   -    [Tinitial,Tfinal]  = time range of integration
202        (if Tinitial>Tfinal the integration is performed backwards in time)
203   -    RelTol, AbsTol = user precribed accuracy
204   - SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
205        returns Ydot = Y' = F(T,Y)
206   - SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
207        returns Jcb = dF/dY
208   -    ICNTRL[1:20]    = integer inputs parameters
209   -    RCNTRL[1:20]    = real inputs parameters
210~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211
212~~~>  OUTPUT ARGUMENTS:
213
214   -    Y[NVAR]         -> vector of final states (at T->Tfinal)
215   -    ISTATUS[1:20]   -> integer output parameters
216   -    RSTATUS[1:20]   -> real output parameters
217   -    Ierr            -> job status upon return
218                           success (positive value) or
219                           failure (negative value)
220~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221~~~>  INPUT PARAMETERS:
222    Note: For input parameters equal to zero the default values of the
223       corresponding variables are used.
224
225    Note: For input parameters equal to zero the default values of the
226       corresponding variables are used.
227~~~>
228    ICNTRL[0] = not used
229    ICNTRL[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors
230                                           = 1: AbsTol, RelTol are scalars
231    ICNTRL[2] = Method
232    ICNTRL[3]  -> maximum number of integration steps
233            For ICNTRL[3]=0 the default value of 100000 is used
234    ICNTRL[4]  -> maximum number of Newton iterations
235        For ICNTRL(4)=0 the default value of 8 is used
236    ICNTRL[5]  -> starting values of Newton iterations:
237    ICNTRL[5]=0 : starting values are interpolated (the default)
238    ICNTRL[5]=1 : starting values are zero
239
240~~~>  Real parameters
241
242    RCNTRL[0]  -> Hmin, lower bound for the integration step size
243                  It is strongly recommended to keep Hmin = ZERO
244    RCNTRL[1]  -> Hmax, upper bound for the integration step size
245    RCNTRL[2]  -> Hstart, starting value for the integration step size
246    RCNTRL[3]  -> FacMin, lower bound on step decrease factor (default=0.2)
247    RCNTRL[4]  -> FacMax, upper bound on step increase factor (default=6)
248    RCNTRL[5]  -> FacRej, step decrease factor after multiple rejections
249                  (default=0.1)
250    RCNTRL[6]  -> FacSafe, by which the new step is slightly smaller
251                  than the predicted value  (default=0.9)
252    RCNTRL[7]  -> ThetaMin. If Newton convergence rate smaller
253                  than ThetaMin the Jacobian is not recomputed;
254                  (default=0.001)
255    RCNTRL[8]  -> NewtonTol, stopping criterion for Newton's method
256                  (default=0.03)
257    RCNTRL[9]  -> Qmin
258    RCNTRL[10] -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
259                  step size is kept constant and the LU factorization
260                  reused (default Qmin=1, Qmax=1.2)
261~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
262~~~>  OUTPUT PARAMETERS:
263   Note: each call to Rosenbrock adds the current no. of fcn calls
264   to previous value of ISTATUS(1), and similar for the other params.
265   Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
266  ISTATUS[0] = No. of function calls
267  ISTATUS[1] = No. of jacobian calls
268  ISTATUS[2] = No. of steps
269  ISTATUS[3] = No. of accepted steps
270  ISTATUS[4] = No. of rejected steps (except at the beginning)
271  ISTATUS[5] = No. of LU decompositions
272  ISTATUS[6] = No. of forward/backward substitutions
273  ISTATUS[7] = No. of singular matrix decompositions
274  RSTATUS[0]  -> Texit, the time corresponding to the computed Y upon return
275  RSTATUS[1]  -> Hexit,last accepted step before return
276  RSTATUS[2]  -> Hnew, last predicted step before return
277  For multiple restarts, use Hnew as Hstart in the following run
278
279~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
280{ 
281
282int Max_no_steps=0;
283
284/*~~~> Local variables      */
285int StartNewton;
286
287KPP_REAL  Hmin=0,
288        Hmax=0,
289        Hstart=0,
290        Roundoff,
291        FacMin=0,
292        FacMax=0,
293        FacSafe=0,
294        FacRej=0,
295        ThetaMin,
296        NewtonTol,
297        Qmin,
298        Qmax;
299
300int ITOL,
301        NewtonMaxit,
302        i,
303        Ierr = 0;
304
305/*~~~>  For Scalar tolerances (ICNTRL[1] !=0 )  the code uses AbsTol[1] and RelTol[1)
306    For Vector tolerances (ICNTRL[1] == 0) the code uses AbsTol[1:NVAR] and RelTol[1:NVAR] */
307if (ICNTRL[1]==0){
308        ITOL = 1;
309}
310else {
311        ITOL = 0;
312}  /* end if */
313
314/*~~~> ICNTRL[3] - method selection   */
315
316switch (ICNTRL[2]) {
317
318case 0: Sdirk2a();
319        break;
320
321case 1: Sdirk2a();
322        break;
323
324case 2: Sdirk2b();
325        break;
326
327case 3: Sdirk3a();
328        break;
329
330case 4: Sdirk4a();
331        break;
332
333case 5: Sdirk4b();
334        break;
335
336default: Sdirk2a();
337
338} /* end switch */
339
340/*~~~>   The maximum number of time steps admitted   */
341if (ICNTRL[3] == 0) {
342        Max_no_steps = 200000;
343}
344else if (ICNTRL[3] > 0) {
345        Max_no_steps = ICNTRL[3];
346}
347else {
348        printf("User-selected ICNTRL(4)=%d", ICNTRL[3]);
349        SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr);
350} /*end if */
351
352/*~~~>The maximum number of Newton iterations admitted */
353if(ICNTRL[4]==0) {
354        NewtonMaxit = 8;
355}
356else {
357        NewtonMaxit=ICNTRL[4];
358        if(NewtonMaxit <=0) {
359                printf("User-selected ICNTRL(5)=%d", ICNTRL[4] );
360                SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr);
361        }
362}
363
364/*~~~> StartNewton: Extrapolate for starting values of Newton iterations */
365if (ICNTRL[5] == 0) {
366        StartNewton = 1;
367}
368else {
369        StartNewton = 0;
370} /* end if */
371
372/*~~~>  Unit roundoff (1+Roundoff>1) */
373Roundoff = WLAMCH('E');
374
375/*~~~>  Lower bound on the step size: (positive value) */
376if (RCNTRL[0] == ZERO) {
377        Hmin = ZERO;
378}
379else if (RCNTRL[0] > ZERO) {
380        Hmin = RCNTRL[0];
381}
382else {
383        printf("User-selected RCNTRL[0]=%f", RCNTRL[0]);
384        SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr);
385} /* end if */
386
387/*~~~>  Upper bound on the step size: (positive value) */
388if (RCNTRL[1] == ZERO) {
389        Hmax = ABS(Tfinal-Tinitial);
390}
391else if (RCNTRL[1] > ZERO) {
392        Hmax = MIN( ABS(RCNTRL[1]), ABS(Tfinal-Tinitial) );
393}
394else {
395        printf("User-selected RCNTRL[1]=%f", RCNTRL[1]);
396        SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr);
397} /* end if */
398
399/*~~~>  Starting step size: (positive value) */
400if (RCNTRL[2] == ZERO) {
401        Hstart = MAX( Hmin, Roundoff);
402}
403else if (RCNTRL[2] > ZERO) {
404        Hstart = MIN( ABS(RCNTRL[2]), ABS(Tfinal-Tinitial) );
405}
406else {
407        printf("User-selected Hstart: RCNTRL[2]=%f", RCNTRL[2]);
408        SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr);
409} /* end if */
410
411/*~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax */
412if (RCNTRL[3] == ZERO) {
413        FacMin = (KPP_REAL)0.2;
414}
415else if (RCNTRL[3] > ZERO) {
416        FacMin = RCNTRL[3];
417}
418else {
419        printf("User-selected FacMin: RCNTRL[3]=%f", RCNTRL[3]);
420        SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr);
421} /* end if */
422
423if (RCNTRL[4] == ZERO) {
424        FacMax = (KPP_REAL)10.0;
425}
426else if (RCNTRL[4] > ZERO) {
427        FacMax = RCNTRL[4];
428}
429else {
430        printf("User-selected FacMax: RCNTRL[4]=%f", RCNTRL[4]);
431        SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr);
432} /* end if */
433
434/*~~~>   FacRej: Factor to decrease step after 2 succesive rejections */
435if (RCNTRL[5] == ZERO) {
436        FacRej = (KPP_REAL)0.1;
437}
438else if (RCNTRL[5] > ZERO) {
439        FacRej = RCNTRL[5];
440}
441else {
442        printf("User-selected FacRej: RCNTRL[5]=%f", RCNTRL[5]);
443        SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr);
444} /*end if */
445
446/* ~~~>   FacSafe: Safety Factor in the computation of new step size  */
447if (RCNTRL[6] == ZERO) {
448        FacSafe = (KPP_REAL)0.9;
449}
450else if (RCNTRL[6] > ZERO) {
451        FacSafe = RCNTRL[6];
452}
453else {
454        printf("User-selected FacSafe: RCNTRL[6]=%f", RCNTRL[6]);
455        SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr);
456} /* end if */
457
458/*~~~> ThetaMin: decides whether the Jacobian should be recomputed */
459if (RCNTRL[7] == ZERO) {
460        ThetaMin = (KPP_REAL)1.0e-03;
461}
462else {
463        ThetaMin = RCNTRL[7];
464} /* end if */
465
466/*~~~> Stopping criterion for Newton's method */
467if (RCNTRL[8] == ZERO) {
468        NewtonTol = (KPP_REAL)3.0e-02;
469}
470else {
471        NewtonTol = RCNTRL[8];
472} /* end if */
473
474/* ~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST. */
475if (RCNTRL[9] == ZERO) {
476        Qmin = ONE;
477}
478else {
479        Qmin = RCNTRL[9];
480} /* end if */
481
482if (RCNTRL[10] == ZERO) {
483        Qmax = (KPP_REAL)1.2;
484}
485else {
486        Qmax = RCNTRL [10];
487} /* end if */
488
489/* ~~~>  Check if tolerances are reasonable */
490if (ITOL == 0) {
491        if ((AbsTol[0]<=ZERO || RelTol[0])<=(((KPP_REAL)10.0)*Roundoff)) {
492                SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr);
493        } /* end internal if */
494}
495else {
496        for (i = 0; i < N; i++) {
497                if((AbsTol[i]<=ZERO)||(RelTol[i]<=((KPP_REAL)10.0)*Roundoff)){
498                        SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr);
499                } /* end internal if */
500        } /* end for */
501} /* end if */
502
503if (Ierr < 0) {
504        return Ierr;
505} /*end if */
506
507Ierr = SDIRK_Integrator(N, Tinitial, Tfinal, Y, Ierr, Hstart, Hmin, Hmax,
508        Roundoff, AbsTol, RelTol, ISTATUS, RSTATUS, ITOL, Max_no_steps,
509        StartNewton, NewtonTol, ThetaMin, FacSafe, FacMin, FacMax, FacRej,
510        Qmin, Qmax, NewtonMaxit);
511
512return Ierr;
513
514}  /* end of main SDIRK function */
515
516/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
517int SDIRK_Integrator(int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[], 
518        int Ierr, KPP_REAL Hstart, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Roundoff,
519        KPP_REAL AbsTol[], KPP_REAL RelTol[], int ISTATUS[], KPP_REAL RSTATUS[],
520        int ITOL, int Max_no_steps, int StartNewton, KPP_REAL NewtonTol,
521        KPP_REAL ThetaMin, KPP_REAL FacSafe, KPP_REAL FacMin, KPP_REAL FacMax,
522        KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int NewtonMaxit)
523/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
524{
525
526/*~~~> Local variables:    */
527KPP_REAL Z[Smax][NVAR],
528        G[NVAR],
529        TMP[NVAR],
530        NewtonRate,
531        SCAL[NVAR],
532        RHS[NVAR],
533        T,
534        H,
535        Theta=0,
536        Hratio,
537        NewtonPredictedErr,
538        Qnewton,
539        Err=0,
540        Fac,
541        Hnew,
542        Tdirection,
543        NewtonIncrement=0,
544        NewtonIncrementOld=0;
545
546int IER=0,
547        istage,
548        NewtonIter,
549        IP[NVAR],
550        Reject,
551        FirstStep,
552        SkipJac,
553        SkipLU,
554        NewtonDone,
555        CycleTloop,
556        i,
557        j;
558
559#ifdef FULL_ALGEBRA
560        KPP_REAL FJAC[NVAR][NVAR];
561        KPP_REAL E[NVAR][NVAR];
562#else
563        KPP_REAL FJAC[LU_NONZERO];
564        KPP_REAL E[LU_NONZERO];
565#endif
566
567/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
568/*~~~~>   Initializations              */
569/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
570T = Tinitial;
571Tdirection = SIGN(ONE, Tfinal-Tinitial);
572H = MAX(ABS(Hmin), ABS(Hstart));
573
574if(ABS(H) <= ((KPP_REAL)10.0 * Roundoff)) {
575        H = (KPP_REAL)(1.0e-06);
576} /* end if */
577
578H = MIN(ABS(H), Hmax);
579H = SIGN(H, Tdirection);
580SkipLU = 0;
581SkipJac = 0;
582Reject = 0;
583FirstStep = 1;
584CycleTloop = 0;
585
586SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL);
587
588/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
589/*~~~>  Time loop begins                */
590/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
591while((Tfinal-T)*Tdirection - Roundoff > ZERO) {  /* Tloop */
592
593/*~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition */
594        if(SkipLU == 0) { /* This time around skip the Jac update and LU */
595                SDIRK_PrepareMatrix(H, T, Y, FJAC, SkipJac, SkipLU, E, IP, 
596                        Reject, IER, ISTATUS);
597                if(IER != 0) {
598                        SDIRK_ErrorMsg(-8, T, H, Ierr);
599                return Ierr;
600                } /* end if */
601        } /* end if */
602
603        if(ISTATUS[Nstp] > Max_no_steps) {
604                SDIRK_ErrorMsg(-6, T, H, Ierr);
605                return Ierr;
606        } /* end if */
607
608        if((T + ((KPP_REAL)0.1) * H == T) || (ABS(H) <= Roundoff)) {
609                SDIRK_ErrorMsg(-7, T, H, Ierr);
610                return Ierr;
611        } /* end if */
612
613/*stages*/ for(istage=0; istage < rkS; istage++) {
614                /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
615                /*~~~>  Simplified Newton iterations          */
616                /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
617
618                /*~~~>  Starting values for Newton iterations */
619                Set2Zero(N, &Z[istage][0]);
620
621        /*~~~>   Prepare the loop-independent part of the right-hand side */
622                Set2Zero(N, G);
623                        if(istage > 0) {
624                                for(j=0; j < istage; j++) {
625                                         WAXPY(N, rkTheta[j][istage], 
626                                                &Z[j][0], 1, G, 1);
627                                        if(StartNewton == 1) {
628                                                WAXPY(N, rkAlpha[j][istage],
629                                                        &Z[j][0], 1, 
630                                                        &Z[istage][0], 1);
631                                        } /* end if */
632                                } /* end for */
633                        } /* end if */
634
635        /*~~~>  Initializations for Newton iteration */
636                NewtonDone = 0; /* false */
637                Fac = (KPP_REAL)0.5; /* Step reduction factor */
638
639/*NewtonLoop*/  for(NewtonIter=0; NewtonIter<NewtonMaxit; NewtonIter++ ) {
640
641        /*~~~>  Prepare the loop-dependent part of the right-hand side */
642                        WADD(N, Y, &Z[istage][0], TMP);
643                        FUN_CHEM(T+rkC[istage]*H, TMP, RHS);
644                        ISTATUS[Nfun]++;
645                        WSCAL(N, H*rkGamma, RHS, 1);
646                        WAXPY(N, -ONE, &Z[istage][0], 1, RHS, 1);
647                        WAXPY(N, ONE, G, 1, RHS, 1 );
648
649                /*~~~>   Solve the linear system  */
650                        SDIRK_Solve(H, N, E, IP, IER, RHS, ISTATUS);
651
652                /*~~~>   Check convergence of Newton iterations */
653                        NewtonIncrement = SDIRK_ErrorNorm(N, RHS, SCAL);
654
655                        if(NewtonIter == 0) {
656                                Theta = ABS(ThetaMin);
657                                NewtonRate = (KPP_REAL)2.0;
658                        }
659                        else {
660                                Theta = NewtonIncrement/NewtonIncrementOld;
661
662                                if(Theta < (KPP_REAL)0.99) {
663                                        NewtonRate = Theta/(ONE-Theta);
664                        /* Predict error at the end of Newton process */
665                                        NewtonPredictedErr =   
666                                                (NewtonIncrement*pow(Theta,
667                                                (NewtonMaxit - (NewtonIter +
668                                                1)) / (ONE - Theta)));
669                                        if(NewtonPredictedErr >= NewtonTol) {
670                              /* Non-convergence of Newton:
671                                predicted error too large*/
672                                                Qnewton = MIN((KPP_REAL)10.0,
673                                                        NewtonPredictedErr/
674                                                        NewtonTol);
675                                                Fac = (KPP_REAL)0.8 * pow
676                                                        (Qnewton, (-ONE / (1
677                                                        + NewtonMaxit -
678                                                        NewtonIter + 1)));
679                                                break;
680                                        } /* end internal if */
681                                }
682                                else  /* Non-convergence of Newton:
683                                        Theta too large */ {
684                                        break;
685                                } /* end internal if else */
686                        } /* end if else */
687
688                        NewtonIncrementOld = NewtonIncrement;
689
690                /* Update solution: Z(:) <-- Z(:)+RHS(:) */
691                        WAXPY(N, ONE, RHS, 1, &Z[istage][0], 1);
692
693                /* Check error in Newton iterations */
694                        NewtonDone=(NewtonRate*NewtonIncrement<=NewtonTol);
695
696                        if(NewtonDone == 1) {
697                                break;
698                        }
699                } /* end NewtonLoop for */
700
701                if(NewtonDone == 0) {
702                /* CALL RK_ErrorMsg(-12,T,H,Ierr); */
703                        H = Fac*H;
704                        Reject = 1; /* true */
705                        SkipJac = 1;/* true */
706                        SkipLU = 0;/* false */
707                        CycleTloop = 1; /* cycle Tloop */
708                } /* end if */
709
710                if(CycleTloop == 1) {
711                        CycleTloop=0;
712                        break;
713                }
714        } /* end stages for */
715
716        if(CycleTloop==0) {
717
718/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
719/*~~~>  Error estimation      */
720/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
721        ISTATUS[Nstp]++;
722        Set2Zero(N, TMP);
723
724        for(i=0; i<rkS; i++) {
725                if(rkE[i] != ZERO) {
726                        WAXPY(N, rkE[i], &Z[i][0], 1, TMP, 1);
727                } /* end if */
728        } /* end for */
729
730        SDIRK_Solve(H, N, E, IP, IER, TMP, ISTATUS);
731        Err = SDIRK_ErrorNorm(N, TMP, SCAL);
732
733/*~~~~> Computation of new step size Hnew */
734        Fac = FacSafe * pow((Err), (-ONE/rkELO));
735        Fac = MAX(FacMin, MIN(FacMax, Fac));
736        Hnew = H*Fac;
737
738/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
739/*~~~>  Accept/Reject step    */
740/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
741
742        if(Err < ONE) { /*~~~> Step is accepted */
743                FirstStep = 0; /* false */
744                ISTATUS[Nacc]++;
745
746        /*~~~> Update time and solution */
747                T = T + H;
748
749        /* Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) */
750                for(i=0; i<rkS; i++) {
751                        if(rkD[i] != ZERO) {
752                                WAXPY(N, rkD[i], &Z[i][0], 1, Y, 1);
753                        } /* end if */
754                } /* end for */
755
756        /*~~~> Update scaling coefficients */
757                SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL);
758
759        /*~~~> Next time step */
760                Hnew = Tdirection*MIN(ABS(Hnew), Hmax);
761
762        /* Last T and H */
763                RSTATUS[Ntexit] = T;
764                RSTATUS[Nhexit] = H;
765                RSTATUS[Nhnew] = Hnew;
766
767        /* No step increase after a rejection */
768                if(Reject==1) {
769                        Hnew = Tdirection*MIN(ABS(Hnew), ABS(H));
770                } /* end if */
771
772                Reject = 0; /* false */
773
774                if((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) {
775                        H = Tfinal-T;
776                }
777                else {
778                        Hratio = Hnew/H;
779                /* If step not changed too much keep Jacobian and reuse LU */
780                        SkipLU = ((Theta <= ThetaMin) && (Hratio >= Qmin) && 
781                                (Hratio <= Qmax));
782
783                        if(SkipLU==0) {
784                                H = Hnew;
785                        } /* end internal if */
786                } /* end if else */
787
788                SkipJac = (Theta <= ThetaMin);
789                SkipJac = 0; /* false */
790                }
791                else { /*~~~> Step is rejected */
792                        if((FirstStep==1) || (Reject==1)) {
793                                H = FacRej * H;
794                        }
795                        else {
796                                H = Hnew;
797                        } /* end internal if */
798
799                        Reject = 1;
800                        SkipJac = 1;
801                        SkipLU = 0;
802
803                        if(ISTATUS[Nacc] >=1) {
804                                ISTATUS[Nrej]++;
805                        } /* end if */
806                } /* end if else */
807        } /* end CycleTloop if */
808} /* end Tloop */
809
810/* Successful return */
811Ierr = 1;
812return Ierr;
813
814} /* end SDIRK_Integrator */
815
816/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
817void SDIRK_ErrorScale(int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[],
818                        KPP_REAL Y[],KPP_REAL SCAL[])
819/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
820{
821
822int i;
823if (ITOL == 0){
824        for (i = 0; i < NVAR; i++){
825                SCAL[i] = ONE / (AbsTol[0]+RelTol[0]*ABS(Y[i]) );
826        } /* end for */
827}
828else {
829        for (i = 0; i < NVAR; i++){
830                SCAL[i] = ONE / (AbsTol[i]+RelTol[i]*ABS(Y[i]) );
831        } /*  end for  */
832} /*  end if  */
833
834} /*  end SDIRK_ErrorScale  */
835
836/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
837KPP_REAL SDIRK_ErrorNorm(int N, KPP_REAL Y[], KPP_REAL SCAL[])
838/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
839{
840
841int i;
842KPP_REAL Err = ZERO;
843
844for (i = 0; i < N; i++) {
845        Err = Err + pow( (Y[i]*SCAL[i]), 2);
846} /* end for */
847
848Err = MAX( SQRT(Err/(KPP_REAL)N), (KPP_REAL)1.0e-10);
849
850return Err;
851
852} /*  end SDIRK_ErrorNorm  */
853
854/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
855int SDIRK_ErrorMsg(int code, KPP_REAL T, KPP_REAL H, int Ierr)
856/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
857 *    Handles all error messages
858 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
859{
860
861Ierr = code;
862
863printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
864printf("\nForced exit from Sdirk due to the following error:\n");
865
866switch (code) {
867
868case -1:
869        printf("--> Improper value for maximal no of steps");
870        break;
871case -2:
872        printf("--> Selected Rosenbrock method not implemented");
873        break;
874case -3:
875        printf("--> Hmin/Hmax/Hstart must be positive");
876        break;
877case -4:
878        printf("--> FacMin/FacMax/FacRej must be positive");
879        break;
880case -5:
881        printf("--> Improper tolerance values");
882        break;
883case -6:
884        printf("--> No of steps exceeds maximum bound");
885        break;
886case -7:
887        printf("--> Step size too small (T + H/10 = T) or H < Roundoff");
888        break;
889case -8:
890        printf("--> Matrix is repeatedly singular");
891        break;
892default: /* causing an error */
893        printf("Unknown Error code: %d", code);
894
895} /*  end switch   */
896
897printf("\n   Time = %f and H = %f", T, H );
898printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
899
900return code;
901
902} /*  end SDIRK_ErrorMsg   */
903
904/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
905void SDIRK_PrepareMatrix(KPP_REAL H, KPP_REAL T, KPP_REAL Y[], KPP_REAL FJAC[],
906                        int SkipJac, int SkipLU, KPP_REAL E[], int IP[], 
907                        int Reject, int ISING, int ISTATUS[] )
908/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
909 *  Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition
910 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
911{
912
913KPP_REAL HGammaInv;
914int i, j;
915int ConsecutiveSng = 0;
916ISING = 1;
917
918
919while( ISING != 0) {
920        HGammaInv = ONE/(H*rkGamma);
921
922/*~~~>  Compute the Jacobian */
923        if(SkipJac==0) {
924                JAC_CHEM(T,Y,FJAC);
925                ISTATUS[Njac]++;
926        } /* end if */
927
928        #ifdef FULL_ALGEBRA
929                for(j=0; j<NVAR; j++) {
930                        for(i=0; i<NVAR; i++) {
931                                E[j][i] = -FJAC[j][i];
932                        } /* end for */
933
934                        E[j][j] = E[j][j] + HGammaInv;
935                } /* end for */
936
937                DGETRF(NVAR, NVAR, E, NVAR, IP, ISING);
938        #else
939                for(i=0; i<LU_NONZERO; i++) {
940                        E[i] = -FJAC[i];
941                } /* end for */
942
943                for(i=0; i<NVAR; i++) {
944                        j = LU_DIAG[i];
945                        E[j]=E[j] + HGammaInv;
946                } /* end for */
947
948                ISING = KppDecomp(E);
949                IP[0] = 1;
950        #endif
951
952        ISTATUS[Ndec]++;
953
954        if(ISING != 0) {
955                ISTATUS[Nsng]++;
956                ConsecutiveSng++;
957
958                if(ConsecutiveSng >= 6) {
959                        return; /* Failure */
960                } /* end internal if */
961
962                H = (KPP_REAL)(0.5) * H;
963                SkipJac = 1; /* true */
964                SkipLU = 0; /* false */
965                Reject = 1; /* true */
966        } /* end if */
967} /* end while */
968
969} /* end SDIRK_PrepareMatrix */
970
971/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
972void SDIRK_Solve( KPP_REAL H, int N, KPP_REAL E[], int IP[], int ISING, 
973                KPP_REAL RHS[], int ISTATUS[] )
974/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
975 *  Solves the system (H*Gamma-Jac)*x = RHS
976 *      using the LU decomposition of E = I - 1/(H*Gamma)*Jac
977 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
978{
979
980KPP_REAL HGammaInv = ONE/(H * rkGamma);
981
982WSCAL(N, HGammaInv, RHS, 1);
983
984#ifdef FULL_ALGEBRA
985        DGETRS('N', N, 1, E, N, IP, RHS, N, ISING);
986#else
987        KppSolve(E, RHS);
988#endif
989        ISTATUS[Nsol]++;
990
991} /* end SDIRK_Solve */
992
993/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
994void Sdirk4a()
995/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
996{
997
998sdMethod = S4A;
999
1000/* Number of stages */
1001rkS = 5;
1002
1003/* Method Coefficients */
1004rkGamma = (KPP_REAL)0.2666666666666666666666666666666667;
1005
1006rkA[0][0] = (KPP_REAL)0.2666666666666666666666666666666667;
1007rkA[0][1] = (KPP_REAL)0.5000000000000000000000000000000000;
1008rkA[1][1] = (KPP_REAL)0.2666666666666666666666666666666667;
1009rkA[0][2] = (KPP_REAL)0.3541539528432732316227461858529820;
1010rkA[1][2] = (KPP_REAL)(-0.5415395284327323162274618585298197e-01);
1011rkA[2][2] = (KPP_REAL)0.2666666666666666666666666666666667;
1012rkA[0][3] = (KPP_REAL)0.8515494131138652076337791881433756e-01;
1013rkA[1][3] = (KPP_REAL)(-0.6484332287891555171683963466229754e-01);
1014rkA[2][3] = (KPP_REAL)0.7915325296404206392428857585141242e-01;
1015rkA[3][3] = (KPP_REAL)0.2666666666666666666666666666666667;
1016rkA[0][4] = (KPP_REAL)2.100115700566932777970612055999074;
1017rkA[1][4] = (KPP_REAL)(-0.7677800284445976813343102185062276);
1018rkA[2][4] = (KPP_REAL)2.399816361080026398094746205273880;
1019rkA[3][4] = (KPP_REAL)(-2.998818699869028161397714709433394);
1020rkA[4][4] = (KPP_REAL)0.2666666666666666666666666666666667;
1021rkB[0] = (KPP_REAL)2.100115700566932777970612055999074;
1022rkB[1] = (KPP_REAL)(-0.7677800284445976813343102185062276);
1023rkB[2] = (KPP_REAL)2.399816361080026398094746205273880;
1024rkB[3] = (KPP_REAL)(-2.998818699869028161397714709433394);
1025rkB[4] = (KPP_REAL)0.2666666666666666666666666666666667;
1026
1027rkBhat[0] = (KPP_REAL)2.885264204387193942183851612883390;
1028rkBhat[1] = (KPP_REAL)(-0.1458793482962771337341223443218041);
1029rkBhat[2] = (KPP_REAL)2.390008682465139866479830743628554;
1030rkBhat[3] = (KPP_REAL)(-4.129393538556056674929560012190140);
1031rkBhat[4] = ZERO;
1032
1033rkC[0] = (KPP_REAL)0.2666666666666666666666666666666667;
1034rkC[1] = (KPP_REAL)0.7666666666666666666666666666666667;
1035rkC[2] = (KPP_REAL)0.5666666666666666666666666666666667;
1036rkC[3] = (KPP_REAL)0.3661315380631796996374935266701191;
1037rkC[4] = ONE;
1038
1039/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1040rkD[0] = ZERO;
1041rkD[1] = ZERO;
1042rkD[2] = ZERO;
1043rkD[3] = ZERO;
1044rkD[4] = ONE;
1045
1046/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1047rkE[0] = (KPP_REAL)(-0.6804000050475287124787034884002302);
1048rkE[1] = (KPP_REAL)(1.558961944525217193393931795738823);
1049rkE[2] = (KPP_REAL)(-13.55893003128907927748632408763868);
1050rkE[3] = (KPP_REAL)(15.48522576958521253098585004571302);
1051rkE[4] = ONE;
1052
1053/* Local order of Err estimate */
1054rkELO = 4;
1055
1056/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1057rkTheta[0][1] = (KPP_REAL)1.875000000000000000000000000000000;
1058rkTheta[0][2] = (KPP_REAL)1.708847304091539528432732316227462;
1059rkTheta[1][2] = (KPP_REAL)(-0.2030773231622746185852981969486824);
1060rkTheta[0][3] = (KPP_REAL)0.2680325578937783958847157206823118;
1061rkTheta[1][3] = (KPP_REAL)(-0.1828840955527181631794050728644549);
1062rkTheta[2][3] = (KPP_REAL)0.2968246986151577397160821594427966;
1063rkTheta[0][4] = (KPP_REAL)0.9096171815241460655379433581446771;
1064rkTheta[1][4] = (KPP_REAL)(-3.108254967778352416114774430509465);
1065rkTheta[2][4] = (KPP_REAL)12.33727431701306195581826123274001;
1066rkTheta[3][4] = (KPP_REAL)(-11.24557012450885560524143016037523);
1067
1068/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
1069rkAlpha[0][1] = (KPP_REAL)2.875000000000000000000000000000000;
1070rkAlpha[0][2] = (KPP_REAL)0.8500000000000000000000000000000000;
1071rkAlpha[1][2] = (KPP_REAL)0.4434782608695652173913043478260870;
1072rkAlpha[0][3] = (KPP_REAL)0.7352046091658870564637910527807370;
1073rkAlpha[1][3] = (KPP_REAL)(-0.9525565003057343527941920657462074e-01);
1074rkAlpha[2][3] = (KPP_REAL)0.4290111305453813852259481840631738;
1075rkAlpha[0][4] = (KPP_REAL)(-16.10898993405067684831655675112808);
1076rkAlpha[1][4] = (KPP_REAL)6.559571569643355712998131800797873;
1077rkAlpha[2][4] = (KPP_REAL)(-15.90772144271326504260996815012482);
1078rkAlpha[3][4] = (KPP_REAL)25.34908987169226073668861694892683;
1079
1080rkELO = (KPP_REAL)4.0;
1081
1082} /* end Sdirk4a */
1083
1084/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1085void Sdirk4b()
1086/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1087{
1088
1089sdMethod = S4B;
1090
1091/* Number of stages */
1092rkS = 5;
1093
1094/* Method coefficients */
1095rkGamma = (KPP_REAL)0.25;
1096
1097rkA[0][0] = (KPP_REAL)0.25;
1098rkA[0][1] = (KPP_REAL)0.5;
1099rkA[1][1] = (KPP_REAL)0.25;
1100rkA[0][2] = (KPP_REAL)0.34;
1101rkA[1][2] = (KPP_REAL)(-0.40e-01);
1102rkA[2][2] = (KPP_REAL)0.25;
1103rkA[0][3] = (KPP_REAL)0.2727941176470588235294117647058824;
1104rkA[1][3] = (KPP_REAL)(-0.5036764705882352941176470588235294e-01);
1105rkA[2][3] = (KPP_REAL)0.2757352941176470588235294117647059e-01;
1106rkA[3][3] = (KPP_REAL)0.25;
1107rkA[0][4] = (KPP_REAL)1.041666666666666666666666666666667;
1108rkA[1][4] = (KPP_REAL)(-1.020833333333333333333333333333333);
1109rkA[2][4] = (KPP_REAL)7.812500000000000000000000000000000;
1110rkA[3][4] = (KPP_REAL)(-7.083333333333333333333333333333333);
1111rkA[4][4] = (KPP_REAL)0.25;
1112
1113rkB[0] = (KPP_REAL)1.041666666666666666666666666666667;
1114rkB[1] = (KPP_REAL)(-1.020833333333333333333333333333333);
1115rkB[2] = (KPP_REAL)7.812500000000000000000000000000000;
1116rkB[3] = (KPP_REAL)(-7.083333333333333333333333333333333);
1117rkB[4] = (KPP_REAL)0.250000000000000000000000000000000;
1118
1119rkBhat[0] = (KPP_REAL)1.069791666666666666666666666666667;
1120rkBhat[1] = (KPP_REAL)(-0.894270833333333333333333333333333);
1121rkBhat[2] = (KPP_REAL)7.695312500000000000000000000000000;
1122rkBhat[3] = (KPP_REAL)(-7.083333333333333333333333333333333);
1123rkBhat[4] = (KPP_REAL)0.212500000000000000000000000000000;
1124
1125rkC[0] = (KPP_REAL)0.25;
1126rkC[1] = (KPP_REAL)0.75;
1127rkC[2] = (KPP_REAL)0.55;
1128rkC[3] = (KPP_REAL)0.5;
1129rkC[4] = ONE;
1130
1131/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1132rkD[0] = ZERO;
1133rkD[1] = ZERO;
1134rkD[2] = ZERO;
1135rkD[3] = ZERO;
1136rkD[4] = ONE;
1137
1138/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1139rkE[0] = (KPP_REAL)0.5750;
1140rkE[1] = (KPP_REAL)0.2125;
1141rkE[2] = (KPP_REAL)(-4.6875);
1142rkE[3] = (KPP_REAL)4.2500;
1143rkE[4] = (KPP_REAL)0.1500;
1144
1145/* Local order of Err estimate */
1146rkELO = 4;
1147
1148/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1149rkTheta[0][1] = (KPP_REAL)2.0;
1150rkTheta[0][2] = (KPP_REAL)1.680000000000000000000000000000000;
1151rkTheta[1][2] = (KPP_REAL)(-0.1600000000000000000000000000000000);
1152rkTheta[0][3] = (KPP_REAL)1.308823529411764705882352941176471;
1153rkTheta[1][3] = (KPP_REAL)(-0.1838235294117647058823529411764706);
1154rkTheta[2][3] = (KPP_REAL)0.1102941176470588235294117647058824;
1155rkTheta[0][4] = (KPP_REAL)(-3.083333333333333333333333333333333);
1156rkTheta[1][4] = (KPP_REAL)(-4.291666666666666666666666666666667);
1157rkTheta[2][4] = (KPP_REAL)34.37500000000000000000000000000000;
1158rkTheta[3][4] = (KPP_REAL)(-28.3333333333333333333333333333);
1159
1160/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
1161rkAlpha[0][1] = (KPP_REAL)3.0;
1162rkAlpha[0][2] = (KPP_REAL)0.8800000000000000000000000000000000;
1163rkAlpha[1][2] = (KPP_REAL)0.4400000000000000000000000000000000;
1164rkAlpha[0][3] = (KPP_REAL)0.1666666666666666666666666666666667;
1165rkAlpha[1][3] = (KPP_REAL)(-0.8333333333333333333333333333333333e-01);
1166rkAlpha[2][3] = (KPP_REAL)0.9469696969696969696969696969696970;
1167rkAlpha[0][4] = (KPP_REAL)(-6.0);
1168rkAlpha[1][4] = (KPP_REAL)9.0;
1169rkAlpha[2][4] = (KPP_REAL)(-56.81818181818181818181818181818182);
1170rkAlpha[3][4] = (KPP_REAL)54.0;
1171
1172} /* end Sdirk4b */
1173
1174/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1175void Sdirk2a()
1176/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1177{
1178
1179sdMethod = S2A;
1180
1181/* ~~~> Number of stages    */
1182rkS = 2;
1183
1184/* ~~~> Method coefficients    */
1185rkGamma = (KPP_REAL)0.2928932188134524755991556378951510;
1186rkA[0][0] = (KPP_REAL)0.2928932188134524755991556378951510;
1187rkA[0][1] = (KPP_REAL)0.7071067811865475244008443621048490;
1188rkA[1][1] = (KPP_REAL)0.2928932188134524755991556378951510;
1189rkB[0] = (KPP_REAL)0.7071067811865475244008443621048490;
1190rkB[1] = (KPP_REAL)0.2928932188134524755991556378951510;
1191rkBhat[0] = (KPP_REAL)0.6666666666666666666666666666666667;
1192rkBhat[1] = (KPP_REAL)0.3333333333333333333333333333333333;
1193rkC[0] = (KPP_REAL)0.292893218813452475599155637895151;
1194rkC[1] = ONE;
1195
1196/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
1197rkD[0] = ZERO;
1198rkD[1] = ONE;
1199
1200/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
1201rkE[0] = (KPP_REAL)0.4714045207910316829338962414032326;
1202rkE[1] = (KPP_REAL)(-0.1380711874576983496005629080698993);
1203
1204/* ~~~> Local order of Err estimate    */
1205rkELO = 2;
1206
1207/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
1208rkTheta[0][1] = (KPP_REAL)2.414213562373095048801688724209698;
1209
1210/* ~~~> Starting value for Newton iterations */
1211rkAlpha[0][1] = (KPP_REAL)3.414213562373095048801688724209698;
1212
1213} /* end Sdirk2a */
1214
1215/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1216void Sdirk2b()
1217/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1218{
1219
1220sdMethod = S2B;
1221
1222/* ~~~> Number of stages    */
1223rkS      = 2;
1224
1225/* ~~~> Method coefficients    */
1226rkGamma = (KPP_REAL)1.707106781186547524400844362104849;
1227rkA[0][0] = (KPP_REAL)1.707106781186547524400844362104849;
1228rkA[0][1] = (KPP_REAL)(-0.707106781186547524400844362104849);
1229rkA[1][1] = (KPP_REAL)1.707106781186547524400844362104849;
1230rkB[0] = (KPP_REAL)(-0.707106781186547524400844362104849);
1231rkB[1] = (KPP_REAL)1.707106781186547524400844362104849;
1232rkBhat[0] = (KPP_REAL)0.6666666666666666666666666666666667;
1233rkBhat[1] = (KPP_REAL)0.3333333333333333333333333333333333;
1234rkC[0] = (KPP_REAL)1.707106781186547524400844362104849;
1235rkC[1] = ONE;
1236
1237/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
1238rkD[0] = ZERO;
1239rkD[1] = ONE;
1240
1241/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
1242rkE[0] = (KPP_REAL)(-0.4714045207910316829338962414032326);
1243rkE[1] = (KPP_REAL)0.8047378541243650162672295747365659;
1244
1245/* ~~~> Local order of Err estimate    */
1246rkELO = 2;
1247
1248/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
1249rkTheta[0][1] = (KPP_REAL)(-0.414213562373095048801688724209698);
1250
1251/* ~~~> Starting value for Newton iterations */
1252rkAlpha[0][1] = (KPP_REAL)0.5857864376269049511983112757903019;
1253
1254} /* end Sdirk2b */
1255
1256/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1257void Sdirk3a()
1258/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1259{
1260
1261sdMethod = S3A;
1262
1263/* ~~~> Number of stages    */
1264rkS = 3;
1265
1266/* ~~~> Method coefficients    */
1267rkGamma = (KPP_REAL)0.2113248654051871177454256097490213;
1268rkA[0][0] = (KPP_REAL)0.2113248654051871177454256097490213;
1269rkA[0][1] = (KPP_REAL)0.2113248654051871177454256097490213;
1270rkA[1][1] = (KPP_REAL)0.2113248654051871177454256097490213;
1271rkA[0][2] = (KPP_REAL)0.2113248654051871177454256097490213;
1272rkA[1][2] = (KPP_REAL)0.5773502691896257645091487805019573;
1273rkA[2][2] = (KPP_REAL)0.2113248654051871177454256097490213;
1274rkB[0] = (KPP_REAL)0.2113248654051871177454256097490213;
1275rkB[1] = (KPP_REAL)0.5773502691896257645091487805019573;
1276rkB[2] = (KPP_REAL)0.2113248654051871177454256097490213;
1277rkBhat[0]= (KPP_REAL)0.2113248654051871177454256097490213;
1278rkBhat[1]= (KPP_REAL)0.6477918909913548037576239837516312;
1279rkBhat[2]= (KPP_REAL)0.1408832436034580784969504064993475;
1280rkC[0] = (KPP_REAL)0.2113248654051871177454256097490213;
1281rkC[1] = (KPP_REAL)0.4226497308103742354908512194980427;
1282rkC[2] = ONE;
1283
1284/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
1285rkD[0] = ZERO;
1286rkD[1] = ZERO;
1287rkD[2] = ONE;
1288
1289/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
1290rkE[0] = (KPP_REAL)0.9106836025229590978424821138352906;
1291rkE[1] = (KPP_REAL)(-1.244016935856292431175815447168624);
1292rkE[2] = (KPP_REAL)0.3333333333333333333333333333333333;
1293
1294/* ~~~> Local order of Err estimate    */
1295rkELO    = 2;
1296
1297/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
1298rkTheta[0][1] =  ONE;
1299rkTheta[0][2] = (KPP_REAL)(-1.732050807568877293527446341505872);
1300rkTheta[1][2] = (KPP_REAL)2.732050807568877293527446341505872;
1301
1302/* ~~~> Starting value for Newton iterations */
1303rkAlpha[0][1] = (KPP_REAL)2.0;
1304rkAlpha[0][2] = (KPP_REAL)(-12.92820323027550917410978536602349);
1305rkAlpha[1][2] = (KPP_REAL)8.83012701892219323381861585376468;
1306
1307} /* end Sdirk3a */
1308
1309/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1310void FUN_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL P[])
1311/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1312{
1313
1314KPP_REAL Told;
1315
1316Told = TIME;
1317TIME = T;
1318Update_SUN();
1319Update_RCONST();
1320Fun( Y, FIX, RCONST, P );
1321TIME = Told;
1322
1323
1324}  /*  end FUN_CHEM  */
1325
1326/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1327void JAC_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL JV[])
1328/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1329{
1330
1331KPP_REAL Told;
1332
1333#ifdef FULL_ALGEBRA
1334        KPP_REAL JS[LU_NONZERO];
1335        int i,j;
1336#endif
1337
1338Told = TIME;
1339TIME = T;
1340Update_SUN();
1341Update_RCONST();
1342
1343#ifdef FULL_ALGEBRA
1344        Jac_SP( Y, FIX, RCONST, JS);
1345
1346        for(j=0; j<NVAR; j++) {
1347                for(i=0; i<NVAR; i++) {
1348                        JV[j][i] = (KPP_REAL)0.0;
1349                } /* end for */
1350        } /* end for */
1351
1352        for(i=0; i<LU_NONZERO; i++) {
1353                JV[LU_ICOL[i]][LU_IROW[i]] = JS[i];
1354        } /* end for */
1355#else
1356        Jac_SP(Y, FIX, RCONST, JV);
1357#endif
1358
1359TIME = Told;
1360
1361} /* end JAC_CHEM  */
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |