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_adj.c @ 4446

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

Merge of branch palm4u into trunk
File size: 55.0 KB
Line 
1
2#define MAX(a,b) ( ((a) >= (b)) ? (a):(b)  )
3#define MIN(b,c) ( ((b) < (c))  ? (b):(c)  )
4#define ABS(x)   ( ((x) >= 0 )  ? (x):(-x) )
5#define SIGN(x,y)( ( (x*y) >= 0 ) ?(x):(-x) )
6#define SQRT(d)  ( pow((d),0.5)  )
7
8/* Numerical Constants */
9#define ZERO        (KPP_REAL)0.0
10#define ONE         (KPP_REAL)1.0
11#define HALF        (KPP_REAL)0.5
12#define DeltaMin    (KPP_REAL)1.0e-5
13enum boolean { FALSE=0, TRUE=1 };
14 
15/*~~~>  Statistics on the work performed by the SDIRK method */
16enum statistics { Nfun=1, Njac=2, Nstp=3, Nacc=4, Nrej=5, Ndec=6, Nsol=7, 
17                  Nsng=8, Ntexit=1, Nhexit=2, Nhnew=3 };
18
19/*~~~>  SDIRK method coefficients, up to 5 stages */
20int Smax = 5;
21enum ros_Params { S2A=1, S2B=2, S3A=3, S4A=4, S4B=5 };
22KPP_REAL rkGamma, rkA[5][5], rkB[5], rkC[5], rkD[5], rkE[5], 
23  rkBhat[5], rkELO, rkAlpha[5][5], rkTheta[5][5];
24int sdMethod, rkS; /* The number of stages */
25
26/*~~~>  Checkpoints in memory buffers */
27int stack_ptr = -1; /* last written entry in checkpoint */
28KPP_REAL *chk_H, *chk_T;
29KPP_REAL **chk_Y;
30KPP_REAL ***chk_Z;
31int **chk_P;
32#ifdef FULL_ALGEBRA
33  KPP_REAL ***chk_J;
34#else
35  KPP_REAL **chk_J;
36#endif
37
38/* Function Headers */
39void INTEGRATE_ADJ( int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], 
40                    KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], 
41                    KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], 
42                    KPP_REAL RCNTRL_U[], int ISTATUS_U[], 
43                    KPP_REAL RSTATUS_U[] );
44int SDIRKADJ( int N, int NADJ, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[],
45              KPP_REAL Lambda[][NVAR], KPP_REAL RelTol[], KPP_REAL AbsTol[], 
46              KPP_REAL RelTol_adj[][NVAR], KPP_REAL AbsTol_adj[][NVAR], 
47              KPP_REAL RCNTRL[], int ICNTRL[], KPP_REAL RSTATUS[], 
48              int ISTATUS[]);
49int SDIRK_FwdInt( int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[], 
50                  KPP_REAL Hmax, KPP_REAL Hmin, KPP_REAL Hstart, 
51                  KPP_REAL Roundoff, KPP_REAL AbsTol[], KPP_REAL RelTol[], 
52                  int ISTATUS[], KPP_REAL RSTATUS[], int Max_no_steps, 
53                  int NewtonMaxit, KPP_REAL NewtonTol, KPP_REAL ThetaMin, 
54                  KPP_REAL FacSafe, KPP_REAL FacMin, KPP_REAL FacMax, 
55                  KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int ITOL, 
56                  int SaveLU );
57int SDIRK_DadjInt( int N, int NADJ, KPP_REAL Lambda[][NVAR], int SaveLU, 
58                   int ISTATUS[], int ITOL, KPP_REAL AbsTol_adj[][NVAR], 
59                   KPP_REAL RelTol_adj[][NVAR], int NewtonMaxit,
60                   KPP_REAL ThetaMin, KPP_REAL NewtonTol, int DirectADJ );
61void SDIRK_AllocBuffers( int Max_no_steps, int rkS, int SaveLU );
62void SDIRK_FreeBuffers( int Max_no_steps, int SaveLU );
63void SDIRK_Push( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL Z[][NVAR], 
64                 KPP_REAL E[], int P[], int Max_no_steps, int SaveLU );
65void SDIRK_Pop( KPP_REAL* T, KPP_REAL* H, KPP_REAL* Y, KPP_REAL* Z, 
66                KPP_REAL* E, int* P, int SaveLU );
67void SDIRK_ErrorScale( int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[],
68                       KPP_REAL Y[],KPP_REAL SCAL[]);
69KPP_REAL SDIRK_ErrorNorm( int N, KPP_REAL Y[], KPP_REAL SCAL[] );
70int SDIRK_ErrorMsg( int code, KPP_REAL T, KPP_REAL H );
71void SDIRK_PrepareMatrix( KPP_REAL H, KPP_REAL T, KPP_REAL Y[], KPP_REAL FJAC[],
72                          int SkipJac, int SkipLU, KPP_REAL E[], int IP[], 
73                          int Reject, int ISING, int ISTATUS[] );
74void SDIRK_Solve ( char Transp, KPP_REAL H, int N, KPP_REAL E[], 
75                   int IP[], int ISING, KPP_REAL RHS[], int ISTATUS[] );
76void Sdirk4a();
77void Sdirk4b();
78void Sdirk2a();
79void Sdirk2b();
80void Sdirk3a();
81void FUN_CHEM( KPP_REAL T, KPP_REAL Y[], KPP_REAL P[] );
82void JAC_CHEM( KPP_REAL T, KPP_REAL Y[], KPP_REAL JV[] );
83KPP_REAL WLAMCH( char C );
84void Set2Zero( int N, KPP_REAL Y[] );
85void WAXPY( int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], 
86            int incY );
87void WADD( int N, KPP_REAL Y[], KPP_REAL Z[], KPP_REAL TMP[] );
88void WSCAL( int N, KPP_REAL Alpha, KPP_REAL X[], int incX );
89void JacTR_SP_Vec( KPP_REAL Jac[], KPP_REAL Fcn[], KPP_REAL K[] );
90void KppSolve( KPP_REAL A[], KPP_REAL b[] );
91void KppSolveTR( KPP_REAL JVS[], KPP_REAL X[], KPP_REAL XX[] );
92int KppDecomp( KPP_REAL A[] );
93void Update_SUN();
94void Update_RCONST();
95void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] );
96void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] );
97void Set2Zero( int N, KPP_REAL Y[] );
98
99/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
100void INTEGRATE_ADJ( int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], 
101                    KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], 
102                    KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], 
103                    KPP_REAL RCNTRL_U[], int ISTATUS_U[], 
104                    KPP_REAL RSTATUS_U[] ) {
105/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
106
107  //int Ntotal = 0;
108  KPP_REAL RCNTRL[20], RSTATUS[20], T1, T2;
109  int ICNTRL[20], ISTATUS[20], Ierr, i;
110
111  for(i=0; i<20; i++) {
112    ICNTRL[i] = 0;
113    RCNTRL[i] = (KPP_REAL)0.0;
114    ISTATUS[i] = 0;
115    RSTATUS[i] = (KPP_REAL)0.0;
116  }
117
118  /*~~~> fine-tune the integrator: */
119  ICNTRL[4] = 8;    /* Max no. of Newton iterations */
120  ICNTRL[6] = 1;   /* Adjoint solution by: 0=Newton, 1=direct */
121  ICNTRL[7] = 1;    /* Save fwd LU factorization: 0 = do *not* save, 1 = save */
122
123  /* If optional parameters are given, and if they are >0,
124     then they overwrite default settings. */ 
125  //if(ICNTRL_U != NULL) { /* Check to see if ICNTRL_U is not NULL */
126  //  for(i=0; i<20; i++) {
127  //    if(ICNTRL_U[i] > 0)
128  //      ICNTRL[i] = ICNTRL_U[i];
129  //  }
130  //}
131  //
132  //if(RCNTRL_U != NULL) { /* Check to see if RCNTRL_U is not NULL */
133  //  for(i=0; i<20; i++) {
134  //    if(RCNTRL_U[i] > 0)
135  //      RCNTRL[i] = RCNTRL_U[i];
136  //  }
137  //}
138 
139  T1 = TIN; 
140  T2 = TOUT;
141  Ierr = SDIRKADJ( NVAR, NADJ, T1, T2, Y, Lambda, RTOL, ATOL, ATOL_adj, 
142                   RTOL_adj, RCNTRL, ICNTRL, RSTATUS, ISTATUS );
143
144  /*~~~> Debug option: number of steps */
145  // Ntotal = Ntotal + ISTATUS(Nstp)
146  // printf( "NSTEP=%d Ntotal=%d O3=%e NO2=%e\n", ISTATUS(Nstp), Ntotal,
147  //          VAR(ind_O3), VAR(ind_NO2) );   
148
149  if (Ierr < 0)
150    printf("SDIRK: Unsuccessful exit at T=%f (Ierr=%d)\n", TIN, Ierr );
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  //}
157  //
158  //if(RSTATUS_U != NULL) { /* Check to see if RSTATUS_U is not NULL */
159  //    for(i=0; i<20; i++)
160  //      RSTATUS_U[i] = RSTATUS[i];
161  //}
162  //
163  //Ierr_U = Ierr;
164 
165} /* END of SUBROUTINE INTEGRATE_ADJ */
166
167/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
168int SDIRKADJ( int N, int NADJ, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[],
169              KPP_REAL Lambda[][NVAR], KPP_REAL RelTol[], KPP_REAL AbsTol[], 
170              KPP_REAL RelTol_adj[][NVAR], KPP_REAL AbsTol_adj[][NVAR], 
171              KPP_REAL RCNTRL[], int ICNTRL[], KPP_REAL RSTATUS[], 
172              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, May 2007
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[0:19]   -> integer output parameters
216-    RSTATUS[0:19]   -> real output parameters
217-    Ierr            -> job status upon return
218                        success (positive value) or
219                        failure (negative value)
220~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221
222~~~>     INPUT PARAMETERS:
223
224    Note: For input parameters equal to zero the default values of the
225       corresponding variables are used.
226
227    Note: For input parameters equal to zero the default values of the
228          corresponding variables are used.
229~~~> 
230    ICNTRL[0] = not used
231
232    ICNTRL[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors
233              = 1: AbsTol, RelTol are scalars
234
235    ICNTRL[2] = Method
236
237    ICNTRL[3]  -> maximum number of integration steps
238        For ICNTRL[3]=0 the default value of 1500 is used
239        Note: use a conservative estimate, since the checkpoint
240              buffers are allocated to hold Max_no_steps
241
242    ICNTRL[4]  -> maximum number of Newton iterations
243        For ICNTRL[4]=0 the default value of 8 is used
244
245    ICNTRL[5]  -> starting values of Newton iterations:
246        ICNTRL[5]=0 : starting values are interpolated (the default)
247        ICNTRL[5]=1 : starting values are zero
248
249    ICNTRL[6]  -> method to solve ADJ equations:
250        ICNTRL[6]=0 : modified Newton re-using LU (the default)
251        ICNTRL[6]=1 : direct solution(additional one LU factorization per stage)
252
253    ICNTRL[7]  -> checkpointing the LU factorization at each step:
254        ICNTRL[7]=0 : do *not* save LU factorization (the default)
255        ICNTRL[7]=1 : save LU factorization
256        Note: if ICNTRL[6]=1 the LU factorization is *not* saved
257
258~~~>  Real parameters
259
260    RCNTRL[0]  -> Hmin, lower bound for the integration step size
261                  It is strongly recommended to keep Hmin = ZERO
262    RCNTRL[1]  -> Hmax, upper bound for the integration step size
263    RCNTRL[2]  -> Hstart, starting value for the integration step size
264
265    RCNTRL[3]  -> FacMin, lower bound on step decrease factor (default=0.2)
266    RCNTRL[4]  -> FacMax, upper bound on step increase factor (default=6)
267    RCNTRL[5]  -> FacRej, step decrease factor after multiple rejections
268                 (default=0.1)
269    RCNTRL[6]  -> FacSafe, by which the new step is slightly smaller
270                  than the predicted value  (default=0.9)
271    RCNTRL[7]  -> ThetaMin. If Newton convergence rate smaller
272                  than ThetaMin the Jacobian is not recomputed;
273                  (default=0.001)
274    RCNTRL[8]  -> NewtonTol, stopping criterion for Newton's method
275                  (default=0.03)
276    RCNTRL[9] -> Qmin
277    RCNTRL[10] -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
278                  step size is kept constant and the LU factorization
279                  reused (default Qmin=1, Qmax=1.2)
280
281~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282
283~~~>     OUTPUT PARAMETERS:
284
285    Note: each call to Rosenbrock adds the current no. of fcn calls
286      to previous value of ISTATUS[0], and similar for the other params.
287      Set ISTATUS[1:10] = 0 before call to avoid this accumulation.
288
289    ISTATUS[0] = No. of function calls
290    ISTATUS[1] = No. of jacobian calls
291    ISTATUS[2] = No. of steps
292    ISTATUS[3] = No. of accepted steps
293    ISTATUS[4] = No. of rejected steps (except at the beginning)
294    ISTATUS[5] = No. of LU decompositions
295    ISTATUS[6] = No. of forward/backward substitutions
296    ISTATUS[7] = No. of singular matrix decompositions
297
298    RSTATUS[0]  -> Texit, the time corresponding to the
299                   computed Y upon return
300    RSTATUS[1]  -> Hexit,last accepted step before return
301    RSTATUS[2]  -> Hnew, last predicted step before return
302        For multiple restarts, use Hnew as Hstart in the following run
303
304~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
305
306/* Local variables */     
307  KPP_REAL Hmin=0.0, Hmax=0.0, Hstart=0.0, Roundoff, FacMin=0.0, FacMax=0.0, 
308    FacSafe=0.0, FacRej=0.0, ThetaMin, NewtonTol,Qmin, Qmax;
309  int SaveLU, DirectADJ; /* Boolean variables */                 
310  int ITOL, NewtonMaxit, Max_no_steps=0, i, Ierr=0;
311
312  stack_ptr = -1;
313
314/*~~~>  Initialize statistics */
315  for(i=0; i<20; i++) {
316    ISTATUS[i] = 0;
317    RSTATUS[i] = ZERO;
318  }
319
320/*~~~>  For Scalar tolerances (ICNTRL[1] != 0)  the code uses AbsTol[0]
321            and RelTol[0]
322        For Vector tolerances (ICNTRL[1] == 0) the code uses AbsTol[0:NVAR-1]
323            and RelTol[0:NVAR-1]*/
324  if (ICNTRL[1]==0)
325    ITOL = 1;
326  else
327    ITOL = 0;
328
329/*~~~> ICNTRL(3) - method selection */
330  switch (ICNTRL[2]) {
331    case 0:
332    case 1:
333      Sdirk2a();
334      break;
335    case 2:
336      Sdirk2b();
337      break;
338    case 3:
339      Sdirk3a();
340      break;
341    case 4:
342      Sdirk4a();
343      break;
344    case 5:
345      Sdirk4b();
346      break;
347    default:
348      Sdirk2a();
349  }
350     
351/*~~~>   The maximum number of time steps admitted */
352  if (ICNTRL[3] == 0)
353    Max_no_steps = 200000;
354  else if (ICNTRL[3] > 0)
355    Max_no_steps = ICNTRL[3];
356  else {
357    printf("User-selected ICNTRL(4)=%d", ICNTRL[3]);
358    Ierr = SDIRK_ErrorMsg(-1,Tinitial,ZERO);
359  }
360
361/*~~~>The maximum number of Newton iterations admitted */
362  if(ICNTRL[4]==0)
363    NewtonMaxit = 8;
364  else {
365    NewtonMaxit=ICNTRL[4];
366    if(NewtonMaxit <=0) {
367      printf("User-selected ICNTRL(5)=%d", ICNTRL[4] );
368      Ierr = SDIRK_ErrorMsg(-2,Tinitial,ZERO);
369    }
370  }
371
372/*~~~> Solve ADJ equations directly or by Newton iterations */
373  DirectADJ = (ICNTRL[6] == 1);
374 
375/*~~~> Save or not the forward LU factorization */
376  SaveLU = ((ICNTRL[7] != 0) && (DirectADJ == 0));
377
378/*~~~>  Unit roundoff (1+Roundoff>1) */
379  Roundoff = WLAMCH('E');
380
381/*~~~>  Lower bound on the step size: (positive value) */
382  if (RCNTRL[0] == ZERO)
383    Hmin = ZERO;
384  else if (RCNTRL[0] > ZERO)
385    Hmin = RCNTRL[0];
386  else {
387    printf("User-selected RCNTRL[0]=%f", RCNTRL[0]);
388    Ierr = SDIRK_ErrorMsg(-3,Tinitial,ZERO);
389  }
390   
391/*~~~>  Upper bound on the step size: (positive value) */
392  if (RCNTRL[1] == ZERO)
393    Hmax = ABS(Tfinal-Tinitial);
394  else if (RCNTRL[1] > ZERO)
395    Hmax = MIN( ABS(RCNTRL[1]), ABS(Tfinal-Tinitial) );
396  else {
397    printf("User-selected RCNTRL[1]=%f", RCNTRL[1]);
398    Ierr = SDIRK_ErrorMsg(-3,Tinitial,ZERO);
399  }
400
401/*~~~>  Starting step size: (positive value) */
402  if (RCNTRL[2] == ZERO)
403    Hstart = MAX( Hmin, Roundoff);
404  else if (RCNTRL[2] > ZERO)
405    Hstart = MIN( ABS(RCNTRL[2]), ABS(Tfinal-Tinitial) );
406  else {
407    printf("User-selected Hstart: RCNTRL[2]=%f", RCNTRL[2]);
408    Ierr = SDIRK_ErrorMsg(-3,Tinitial,ZERO);
409  }
410
411/*~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax */
412  if (RCNTRL[3] == ZERO)
413    FacMin = (KPP_REAL)0.2;
414  else if (RCNTRL[3] > ZERO)
415    FacMin = RCNTRL[3];
416  else {
417    printf("User-selected FacMin: RCNTRL[3]=%f", RCNTRL[3]);
418    Ierr = SDIRK_ErrorMsg(-4,Tinitial,ZERO);
419  }
420   
421  if (RCNTRL[4] == ZERO)
422    FacMax = (KPP_REAL)10.0;
423  else if (RCNTRL[4] > ZERO)
424    FacMax = RCNTRL[4];
425  else {
426    printf("User-selected FacMax: RCNTRL[4]=%f", RCNTRL[4]);
427    Ierr = SDIRK_ErrorMsg(-4,Tinitial,ZERO);
428  }
429
430/*~~~>   FacRej: Factor to decrease step after 2 succesive rejections */
431  if (RCNTRL[5] == ZERO)
432    FacRej = (KPP_REAL)0.1;
433  else if (RCNTRL[5] > ZERO)
434    FacRej = RCNTRL[5];
435  else {
436    printf("User-selected FacRej: RCNTRL[5]=%f", RCNTRL[5]);
437    Ierr = SDIRK_ErrorMsg(-4,Tinitial,ZERO);
438  }
439
440/* ~~~>   FacSafe: Safety Factor in the computation of new step size  */
441  if (RCNTRL[6] == ZERO)
442    FacSafe = (KPP_REAL)0.9;
443  else if (RCNTRL[6] > ZERO)
444    FacSafe = RCNTRL[6];
445  else {
446    printf("User-selected FacSafe: RCNTRL[6]=%f", RCNTRL[6]);
447    Ierr = SDIRK_ErrorMsg(-4,Tinitial,ZERO);
448  }
449
450/*~~~> ThetaMin: decides whether the Jacobian should be recomputed */
451  if (RCNTRL[7] == ZERO)
452    ThetaMin = (KPP_REAL)1.0e-03;
453  else
454    ThetaMin = RCNTRL[7];
455
456/*~~~> Stopping criterion for Newton's method */
457  if (RCNTRL[8] == ZERO)
458    NewtonTol = (KPP_REAL)3.0e-02;
459  else
460    NewtonTol = RCNTRL[8];
461
462/* ~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST. */
463  if (RCNTRL[9] == ZERO)
464    Qmin = ONE;
465  else
466    Qmin = RCNTRL[9];
467
468  if (RCNTRL[10] == ZERO)
469    Qmax = (KPP_REAL)1.2;
470  else
471    Qmax = RCNTRL [10];
472
473/* ~~~>  Check if tolerances are reasonable */
474  if (ITOL == 0) {
475    if ((AbsTol[0]<=ZERO || RelTol[0])<=(((KPP_REAL)10.0)*Roundoff))
476      Ierr = SDIRK_ErrorMsg(-5,Tinitial,ZERO);
477  }
478  else {
479    for (i = 0; i < N; i++) {
480      if ((AbsTol[i]<=ZERO)||(RelTol[i]<=((KPP_REAL)10.0)*Roundoff))
481        Ierr = SDIRK_ErrorMsg(-5,Tinitial,ZERO);
482    }
483  }
484
485  if (Ierr < 0)
486    return Ierr;
487   
488/*~~~>  Allocate memory buffers */
489  SDIRK_AllocBuffers(Max_no_steps, rkS, SaveLU);
490
491/*~~~>  Call forward integration */
492  Ierr = SDIRK_FwdInt( N, Tinitial, Tfinal, Y, Hmax, Hmin, Hstart, Roundoff, 
493                       AbsTol, RelTol, ISTATUS, RSTATUS, Max_no_steps, 
494                       NewtonMaxit, NewtonTol, ThetaMin, FacSafe, FacMin, 
495                       FacMax, FacRej, Qmin, Qmax, ITOL, SaveLU );
496
497/*~~~>  Call adjoint integration */ 
498  Ierr = SDIRK_DadjInt( N, NADJ, Lambda, SaveLU, ISTATUS, ITOL, AbsTol_adj, 
499                        RelTol_adj, NewtonMaxit, ThetaMin, NewtonTol, 
500                        DirectADJ );
501
502/*~~~>  Free memory buffers */ 
503  SDIRK_FreeBuffers(Max_no_steps, SaveLU);
504
505  return Ierr;
506
507} /* End of main SDIRK_ADJ */
508
509/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
510int SDIRK_FwdInt( int N, KPP_REAL Tinitial, KPP_REAL Tfinal, KPP_REAL Y[], 
511                  KPP_REAL Hmax, KPP_REAL Hmin, KPP_REAL Hstart, 
512                  KPP_REAL Roundoff, KPP_REAL AbsTol[], KPP_REAL RelTol[], 
513                  int ISTATUS[], KPP_REAL RSTATUS[], int Max_no_steps, 
514                  int NewtonMaxit, KPP_REAL NewtonTol, KPP_REAL ThetaMin, 
515                  KPP_REAL FacSafe, KPP_REAL FacMin, KPP_REAL FacMax, 
516                  KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int ITOL, 
517                  int SaveLU ) {
518/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
519     
520/*~~~> Local variables: */   
521  KPP_REAL Z[Smax][NVAR], G[NVAR], TMP[NVAR], NewtonRate, SCAL[NVAR], RHS[NVAR],
522    T, H, Theta=0.0, Hratio, NewtonPredictedErr, Qnewton, Err, Fac, Hnew, 
523    Tdirection, NewtonIncrement, NewtonIncrementOld=0.0;
524  int i, j, IER=0, istage, NewtonIter, IP[NVAR];
525
526  /*Boolean Variables*/
527  int Reject, FirstStep, SkipJac, SkipLU, NewtonDone, CycleTloop;
528     
529#ifdef FULL_ALGEBRA     
530  KPP_REAL FJAC[NVAR][NVAR], E[NVAR][NVAR];
531#else     
532  KPP_REAL FJAC[LU_NONZERO], E[LU_NONZERO];
533#endif
534
535/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
536/*~~~~>   Initializations              */
537/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
538  T = Tinitial;
539  Tdirection = SIGN(ONE, Tfinal-Tinitial);
540  H = MAX(ABS(Hmin), ABS(Hstart));
541  if (ABS(H) <= ((KPP_REAL)10.0 * Roundoff))
542    H = (KPP_REAL)(1.0e-06);
543  H = MIN(ABS(H), Hmax);
544  H = SIGN(H, Tdirection);
545  SkipLU = 0;
546  SkipJac = 0;
547  Reject = 0;
548  FirstStep = 1;
549  CycleTloop = 0;
550
551  SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL);
552
553/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
554/*~~~>  Time loop begins                */
555/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
556  while((Tfinal-T)*Tdirection - Roundoff > ZERO) {  /* Tloop */
557
558  /*~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition */
559    if (SkipLU == 0) { /* This time around skip the Jac update and LU */
560      SDIRK_PrepareMatrix( H, T, Y, FJAC, SkipJac, SkipLU, E, IP, 
561                           Reject, IER, ISTATUS);
562      if (IER != 0)
563        return SDIRK_ErrorMsg(-8, T, H);
564    }
565
566    if (ISTATUS[Nstp] > Max_no_steps)
567      return SDIRK_ErrorMsg(-6, T, H);
568
569    if ((T + ((KPP_REAL)0.1) * H == T) || (ABS(H) <= Roundoff)) {
570      return SDIRK_ErrorMsg(-7, T, H);
571    }
572
573    for (istage=0; istage < rkS; istage++) {  /*stages*/
574/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
575/*~~~>  Simplified Newton iterations          */
576/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
577
578    /*~~~>  Starting values for Newton iterations */
579      Set2Zero(N, &Z[istage][0]);
580
581    /*~~~>   Prepare the loop-independent part of the right-hand side */
582      Set2Zero(N, G);
583      if (istage > 0) {
584        for (j=0; j < istage; j++) {
585        /* Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) */
586          WAXPY(N, rkTheta[j][istage], &Z[j][0], 1, G, 1);
587        /* Zi(:) = sum_j Alpha(i,j)*Zj(:) */
588          WAXPY(N, rkAlpha[j][istage],&Z[j][0], 1, &Z[istage][0], 1);
589        }
590      }
591
592    /*~~~>  Initializations for Newton iteration */
593      NewtonDone = 0; /* false */
594      Fac = (KPP_REAL)0.5; /* Step reduction factor */
595
596      for (NewtonIter=0; NewtonIter<NewtonMaxit; NewtonIter++ ) { /*NewtonLoop*/
597      /*~~~>  Prepare the loop-dependent part of the right-hand side */
598        WADD(N, Y, &Z[istage][0], TMP); /* TMP <- Y + Zi */
599        FUN_CHEM(T+rkC[istage]*H, TMP, RHS); /* RHS <- Run(Y+Zi) */
600        ISTATUS[Nfun]++;
601      /* RHS[0:N-1] = G[0:N-1] - Z[istage][0:N-1] + (H*rkGamma)*RHS[1:N] */
602        WSCAL(N, H*rkGamma, RHS, 1);
603        WAXPY(N, -ONE, &Z[istage][0], 1, RHS, 1);
604        WAXPY(N, ONE, G, 1, RHS, 1 );
605
606      /*~~~>   Solve the linear system  */
607        SDIRK_Solve('N', H, N, E, IP, IER, RHS, ISTATUS);
608           
609      /*~~~>   Check convergence of Newton iterations */
610        NewtonIncrement = SDIRK_ErrorNorm(N, RHS, SCAL);
611        if (NewtonIter == 0) {
612          Theta = ABS(ThetaMin);
613          NewtonRate = (KPP_REAL)2.0;
614        }
615        else {
616          Theta = NewtonIncrement/NewtonIncrementOld;
617          if (Theta < (KPP_REAL)0.99) {
618            NewtonRate = Theta/(ONE-Theta);
619          /* Predict error at the end of Newton process */
620            NewtonPredictedErr = (NewtonIncrement*pow(Theta,
621                                 (NewtonMaxit - (NewtonIter+1))/(ONE - Theta)));
622            if(NewtonPredictedErr >= NewtonTol) {
623            /* Non-convergence of Newton: predicted error too large*/
624              Qnewton = MIN((KPP_REAL)10.0, NewtonPredictedErr/NewtonTol);
625              Fac = (KPP_REAL)0.8 * pow(Qnewton, (-ONE / (1 + NewtonMaxit - 
626                                                          NewtonIter + 1)));
627              break;
628            }
629          }
630          else  /* Non-convergence of Newton: Theta too large */ {
631            break;
632          }
633        }
634
635        NewtonIncrementOld = NewtonIncrement;
636
637      /* Update solution: Z(:) <-- Z(:)+RHS(:) */
638        WAXPY(N, ONE, RHS, 1, &Z[istage][0], 1);
639
640      /* Check error in Newton iterations */
641        NewtonDone=(NewtonRate*NewtonIncrement<=NewtonTol);
642        if (NewtonDone == 1)
643          break;
644      } /* end NewtonLoop for */
645                     
646      if(NewtonDone == 0) {
647        H = Fac*H;
648        Reject = 1;
649        SkipJac = 1;
650        SkipLU = 0;
651        CycleTloop = 1;
652      } /* end if */
653
654      if (CycleTloop == 1) {
655        CycleTloop=0;
656        break;
657      }
658    /* End of implified Newton iterations */
659    } /* end stages for */
660
661    if (CycleTloop==0) {
662/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
663/*~~~>  Error estimation      */
664/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
665      ISTATUS[Nstp]++;
666      Set2Zero(N, TMP);
667
668      for (i=0; i<rkS; i++) {
669        if (rkE[i] != ZERO)
670          WAXPY(N, rkE[i], &Z[i][0], 1, TMP, 1);
671      }
672
673      SDIRK_Solve('N', H, N, E, IP, IER, TMP, ISTATUS);
674      Err = SDIRK_ErrorNorm(N, TMP, SCAL);
675
676    /*~~~~> Computation of new step size Hnew */
677      Fac = FacSafe * pow((Err), (-ONE/rkELO));
678      Fac = MAX(FacMin, MIN(FacMax, Fac));
679      Hnew = H*Fac;
680
681/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
682/*~~~>  Accept/Reject step    */
683/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
684      if (Err < ONE) { /*~~~> Step is accepted */
685        FirstStep = 0; /* false */
686        ISTATUS[Nacc]++;
687
688      /* Checkpoint solution */
689        SDIRK_Push( T, H, Y, Z, E, IP, Max_no_steps, SaveLU );
690
691      /*~~~> Update time and solution */
692        T = T + H;
693      /* Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) */
694        for(i=0; i<rkS; i++) {
695          if(rkD[i] != ZERO)
696            WAXPY(N, rkD[i], &Z[i][0], 1, Y, 1);
697        }
698
699      /*~~~> Update scaling coefficients */
700        SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL);
701
702      /*~~~> Next time step */
703        Hnew = Tdirection*MIN(ABS(Hnew), Hmax);
704
705      /* Last T and H */
706        RSTATUS[Ntexit] = T;
707        RSTATUS[Nhexit] = H;
708        RSTATUS[Nhnew] = Hnew;
709
710      /* No step increase after a rejection */
711        if (Reject==1)
712          Hnew = Tdirection*MIN(ABS(Hnew), ABS(H));
713        Reject = 0; /* false */
714        if ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO)
715          H = Tfinal-T;
716        else {
717          Hratio = Hnew/H;
718        /* If step not changed too much keep Jacobian and reuse LU */
719          SkipLU = ((Theta <= ThetaMin) && (Hratio >= Qmin) 
720                    && (Hratio <= Qmax));
721          if (SkipLU==0)
722            H = Hnew;
723        }
724
725        /* If convergence is fast enough, do not update Jacobian */
726        /* SkipJac = (Theta <= ThetaMin); */
727        SkipJac = 0;
728      }
729      else { /*~~~> Step is rejected */
730        if ((FirstStep==1) || (Reject==1))
731          H = FacRej * H;
732        else 
733          H = Hnew;
734
735        Reject = 1;
736        SkipJac = 1;
737        SkipLU = 0;
738        if (ISTATUS[Nacc] >=1)
739          ISTATUS[Nrej]++;
740
741      }
742    } /* end CycleTloop if */
743  } /* end Tloop */
744
745  /* Successful return */
746  return 1;
747
748} /* end SDIRK_FwdInt */
749         
750/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
751int SDIRK_DadjInt( int N, int NADJ, KPP_REAL Lambda[][NVAR], int SaveLU, 
752                   int ISTATUS[], int ITOL, KPP_REAL AbsTol_adj[][NVAR], 
753                   KPP_REAL RelTol_adj[][NVAR], int NewtonMaxit, 
754                   KPP_REAL ThetaMin, KPP_REAL NewtonTol, int DirectADJ) {
755/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
756
757/*~~~> Local variables: */
758  KPP_REAL Y[NVAR];
759  KPP_REAL Z[Smax][NVAR], U[Smax][NADJ][NVAR], TMP[NVAR], G[NVAR], NewtonRate, 
760    SCAL[NVAR], DU[NVAR], T, H, Theta, NewtonPredictedErr, NewtonIncrement, 
761    NewtonIncrementOld=0.0;
762  int i, j, IER=0, istage, iadj, NewtonIter, IP[NVAR], IP_adj[NVAR];
763  int Reject=0, SkipJac, SkipLU, NewtonDone; /* Boolean */
764     
765#ifdef FULL_ALGEBRA     
766  KPP_REAL E[NVAR][NVAR], Jac[NVAR][NVAR], E_adj[NVAR][NVAR];
767#else     
768  KPP_REAL E[LU_NONZERO], Jac[LU_NONZERO], E_adj[LU_NONZERO];
769#endif     
770
771/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
772~~~>  Time loop begins
773~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
774  while ( stack_ptr > -1 ) { /* Tloop */
775       
776  /*~~~>  Recover checkpoints for stage values and vectors */
777    SDIRK_Pop( &T, &H, &Y[0], &Z[0][0], &E[0], &IP[0], SaveLU );
778
779  /*~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition */
780    if (!SaveLU) {
781      SkipJac = FALSE; 
782      SkipLU = FALSE;
783      SDIRK_PrepareMatrix ( H, T, Y, Jac, SkipJac, SkipLU, E, IP, Reject, IER, 
784                            ISTATUS );
785      if (IER != 0)
786        return SDIRK_ErrorMsg(-8,T,H); 
787    } 
788
789    for( istage = rkS-1; istage >=0; istage--) { /* Stages Loop */
790    /*~~~>  Jacobian of the current stage solution */
791      for(i=0; i<N; i++)
792        TMP[i] = Y[i] + Z[istage][i];
793      JAC_CHEM(T+rkC[istage]*H,TMP,Jac);
794      ISTATUS[Njac]++;
795       
796      if (DirectADJ) {
797#ifdef FULL_ALGEBRA
798        for(i=0; i<N; i++) {
799          for(j=0; j<N; j++)
800            E_adj[i][j] = -Jac[i][j];
801        }
802        for(i=0; i<N; i++)
803          E_adj[i][i] = E_adj[i][i] + ONE/(H*rkGamma);
804
805        DGETRF( N, N, E_adj, N, IP_adj, IER );
806#else
807        for(i=0; i<LU_NONZERO; i++)
808          E_adj[i] = -Jac[i];
809        for(i=0; i<NVAR; i++) {
810          j = LU_DIAG[i];
811          E_adj[j] = E_adj[j] + ONE/(H*rkGamma);
812        }
813        IER = KppDecomp (E_adj);
814#endif
815        ISTATUS[Ndec]++;
816        if (IER != 0) {
817          printf("At stage %d the matrix used in adjoint computation is " 
818                 "singular\n", istage);
819          return SDIRK_ErrorMsg(-8,T,H);
820        }
821      }
822
823      for(iadj=0; iadj<NADJ; iadj++) { /* adj loop */
824       
825      /*~~~> Update scaling coefficients */
826        for(i=0; i<NVAR; i++)
827          SDIRK_ErrorScale(N, ITOL, &AbsTol_adj[iadj][i], &RelTol_adj[iadj][i], 
828                           &Lambda[iadj][i], SCAL);
829     
830      /*~~~>   Prepare the loop-independent part of the right-hand side
831               G(:) = H*Jac^T*( B(i)*Lambda + sum_j A(j,i)*Uj(:) ) */
832        for(i=0; i<N; i++)
833          G[i] = rkB[istage]*Lambda[iadj][i];
834        if (istage+1 < rkS) {
835          for (j=istage+1; j<rkS; j++)
836            WAXPY(N,rkA[istage][j],&U[j][iadj][0],1,G,1);
837        }
838#ifdef FULL_ALGEBRA 
839        TMP = MATMUL(TRANSPOSE(Jac),G); /* DZ <- Jac(Y+Z)*Y_tlm */
840#else     
841        JacTR_SP_Vec ( Jac, G, TMP );   
842#endif     
843        for(i=0; i<N; i++)
844          G[i] = H*TMP[i];
845       
846        if (DirectADJ) {
847          SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, IER, G, ISTATUS );
848          for(i=0; i<N; i++)
849            U[istage][iadj][i] = G[i];
850        } else {
851          /*~~~>  Initializations for Newton iteration */
852          Set2Zero(N,&U[istage][iadj][0]);
853          NewtonDone = FALSE;
854           
855          /* Newton Loop */
856          for( NewtonIter=0; NewtonIter<NewtonMaxit; NewtonIter++) {
857
858          /*~~~>   Prepare the loop-dependent part of the right-hand side */
859#ifdef FULL_ALGEBRA 
860            for(i=0; i<N; i++)
861              TMP = MATMUL(TRANSPOSE(Jac),U[istage][iadj][i]);   
862#else     
863            for(i=0; i<N; i++)
864              JacTR_SP_Vec ( Jac, &U[istage][iadj][i], TMP );   
865#endif     
866            for(i=0; i<N; i++)
867              DU[i] = U[istage][iadj][i] - (H*rkGamma)*TMP[i] - G[i];
868
869          /*~~~>   Solve the linear system */
870            SDIRK_Solve ( 'T', H, N, E, IP, IER, DU, ISTATUS );
871           
872          /*~~~>   Check convergence of Newton iterations */
873            NewtonIncrement = SDIRK_ErrorNorm(N, DU, SCAL);
874            if ( NewtonIter == 0 ) {
875              Theta = ABS(ThetaMin);
876              NewtonRate = (KPP_REAL)2.0;
877            }
878            else {
879              Theta = NewtonIncrement/NewtonIncrementOld;
880              if (Theta < (KPP_REAL)0.99) {
881                NewtonRate = Theta/(ONE-Theta);
882              /* Predict error at the end of Newton process */ 
883                NewtonPredictedErr = NewtonIncrement*
884                  pow(Theta,(NewtonMaxit-NewtonIter)) / (ONE-Theta);
885              /* Non-convergence of Newton: predicted error too large */
886                if (NewtonPredictedErr >= NewtonTol)
887                  break; /* Exit NewtonLoop */
888              } else { /* Non-convergence of Newton: Theta too large */
889                break;
890              }
891            }
892            NewtonIncrementOld = NewtonIncrement;
893          /* Update solution */
894            for(i=0; i<N; i++)
895              U[istage][iadj][i] -= DU[i];
896           
897          /* Check error in Newton iterations */
898            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol);
899          /* AbsTol is often inappropriate for adjoints -
900              we do at least 4 Newton iterations to ensure convergence
901              of all adjoint components */
902            if ((NewtonIter>=4) && NewtonDone)
903              break; /* Exit Newton Loop */ 
904          }
905           
906        /*~~~> If Newton iterations fail employ the direct solution */
907          if (!NewtonDone) {
908            printf("Problems with Newton Adjoint!!!\n");
909#ifdef FULL_ALGEBRA
910            for(i=0; i<N; i++)
911              E_adj[i][j] = -Jac[i][j];
912            for(j=0; j<N; j++)
913              E_adj[j][j] += ONE/(H*rkGamma);
914            DGETRF( N, N, E_adj, N, IP_adj, IER );
915#else
916            for(i=0; i<LU_NONZERO; i++)
917              E_adj[i] = -Jac[i];
918            for(i=0; i<NVAR; i++) {
919              j = LU_DIAG[i];
920              E_adj[j] += ONE/(H*rkGamma);
921            }
922            IER = KppDecomp ( E_adj);
923#endif
924            ISTATUS[Ndec]++;
925            if (IER != 0) {
926              printf("At stage %d the matrix used in adjoint computation is "
927                     "singular", istage);
928              return SDIRK_ErrorMsg(-8,T,H);
929            }
930            SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, IER, G, ISTATUS );
931            for(i=0; i<N; i++)
932              U[istage][iadj][i] = G[i]; 
933          }
934
935        /*~~~>  End of implified Newton iterations */
936        } /* End of DirADJ */
937
938      } /* End of adj */
939 
940    } /* End of stages */
941
942    /*~~~> Update adjoint solution
943           Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) */
944    for (istage=0; istage<rkS; istage++) {
945      for (iadj=0; iadj<NADJ; iadj++) {
946        for(i=0; i<N; i++)
947          Lambda[iadj][i] += U[istage][iadj][i];
948      }
949    }
950  } /* End of Tloop */
951
952  /* Successful return */
953  return 1;
954 
955} /* End of SDIRK_DadjInt */
956
957/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
958void SDIRK_AllocBuffers(int Max_no_steps, int rkS, int SaveLU) {
959/*~~~>  Allocate buffer space for checkpointing
960~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
961  int i,j;
962
963  chk_H = (KPP_REAL*) malloc(Max_no_steps * sizeof(KPP_REAL));
964  if (chk_H == NULL) {
965    printf("Failed allocation of buffer H\n");
966    exit(0);
967  }
968
969  chk_T = (KPP_REAL*) malloc(Max_no_steps * sizeof(KPP_REAL));
970  if (chk_T == NULL) {
971    printf("Failed allocation of buffer T\n");
972    exit(0);
973  }
974
975  chk_Y = (KPP_REAL**) malloc(Max_no_steps * sizeof(KPP_REAL*));
976  if (chk_Y == NULL) {
977    printf("Failed allocation of buffer Y\n");
978    exit(0);
979  }
980  for(i=0; i<Max_no_steps; i++) {
981    chk_Y[i] = (KPP_REAL*) malloc(NVAR * sizeof(KPP_REAL));
982    if (chk_Y[i] == NULL) {
983      printf("Failed allocation of buffer Y\n");
984      exit(0);
985    }
986  }
987   
988  chk_Z = (KPP_REAL***) malloc(Max_no_steps * sizeof(KPP_REAL**));
989  if (chk_Z == NULL) {
990    printf("Failed allocation of buffer Z\n");
991    exit(0);
992  }
993  for(i=0; i<Max_no_steps; i++) {
994    chk_Z[i] = (KPP_REAL**) malloc(rkS * sizeof(KPP_REAL*));
995    if (chk_Z[i] == NULL) {
996      printf("Failed allocation of buffer Z\n");
997      exit(0);
998    }
999    for(j=0; j<rkS; j++) {
1000      chk_Z[i][j] = (KPP_REAL*) malloc(NVAR * sizeof(KPP_REAL));
1001      if (chk_Z[i][j] == NULL) {
1002        printf("Failed allocation of buffer Z\n");
1003        exit(0);
1004      }
1005    }
1006  }
1007
1008  if (SaveLU) {
1009#ifdef FULL_ALGEBRA
1010    chk_J = (KPP_REAL***) malloc(Max_no_steps * sizeof(KPP_REAL**));
1011    if (chk_J == NULL) {
1012      printf("Failed allocation of buffer J\n");
1013      exit(0);
1014    }
1015    for(i=0; i<Max_no_steps; i++) {
1016      chk_J[i] = (KPP_REAL**) malloc(NVAR * sizeof(KPP_REAL*));
1017      if (chk_J[i] == NULL) {
1018        printf("Failed allocation of buffer J\n");
1019        exit(0);
1020      }
1021      for(j=0; j<NVAR; j++) {
1022        chk_J[i][j] = (KPP_REAL*) malloc(NVAR * sizeof(KPP_REAL));
1023        if (chk_J[i][j] == NULL) {
1024          printf("Failed allocation of buffer J\n");
1025          exit(0);
1026        }
1027      }
1028    }
1029#else
1030    chk_J = (KPP_REAL**) malloc(Max_no_steps * sizeof(KPP_REAL*));
1031    if (chk_J == NULL) {
1032      printf("Failed allocation of buffer J\n");
1033      exit(0);
1034    }
1035    for(i=0; i<Max_no_steps; i++) {
1036      chk_J[i] = (KPP_REAL*) malloc(LU_NONZERO * sizeof(KPP_REAL));
1037      if (chk_J[i] == NULL) {
1038        printf("Failed allocation of buffer J\n");
1039        exit(0);
1040      }
1041    }
1042#endif
1043 
1044    chk_P = (int**) malloc(Max_no_steps * sizeof(int*));
1045    if (chk_P == NULL) {
1046      printf("Failed allocation of buffer P\n");
1047      exit(0);
1048    }
1049    for(i=0; i<Max_no_steps; i++) {
1050      chk_P[i] = (int*) malloc(NVAR * sizeof(int));
1051      if (chk_P[i] == NULL) {
1052        printf("Failed allocation of buffer P\n");
1053        exit(0);
1054      }
1055    }
1056  }
1057} /* End of SDIRK_AllocBuffers */
1058
1059/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1060void SDIRK_FreeBuffers(int Max_no_steps, int SaveLU) {
1061/*~~~>  Deallocate buffer space for discrete adjoint
1062~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1063  int i,j;
1064
1065  free(chk_H);
1066  free(chk_T);
1067
1068  for(i=0; i<Max_no_steps; i++)
1069    free(chk_Y[i]);
1070  free(chk_Y);
1071
1072  for(i=0; i<Max_no_steps; i++) {
1073    for(j=0; j<rkS; j++) {
1074      free(chk_Z[i][j]);
1075    }
1076    free(chk_Z[i]);
1077  }
1078  free(chk_Z);
1079
1080  if(SaveLU) {
1081#ifdef FULL_ALGEBRA
1082    for(i=0; i<Max_no_steps; i++) {
1083      for(j=0; j<rkS; j++) {
1084        free(chk_J[i][j]);
1085      }
1086      free(chk_J[i]);
1087    }
1088    free(chk_Z);
1089#else
1090    for(i=0; i<Max_no_steps; i++)
1091      free(chk_J[i]);
1092    free(chk_J);
1093#endif
1094
1095    for(i=0; i<Max_no_steps; i++)
1096      free(chk_P[i]);
1097    free(chk_P);
1098  }
1099} /* End of SDIRK_FreeBuffers */
1100
1101/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1102void SDIRK_Push( KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL Z[][NVAR], 
1103                 KPP_REAL E[], int P[], int Max_no_steps, int SaveLU ) {
1104/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1105~~~> Saves the next trajectory snapshot for discrete adjoints
1106~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1107 
1108  int i,j;
1109
1110  stack_ptr++;
1111  if( stack_ptr > Max_no_steps ) {
1112    printf( "Push failed: buffer overflow");
1113    exit(0);
1114  }
1115
1116  chk_H[ stack_ptr ] = H;
1117  chk_T[ stack_ptr ] = T;
1118  for(i=0; i<NVAR; i++) {
1119    chk_Y[stack_ptr][i] = Y[i];
1120    for(j=0; j<rkS; j++)
1121      chk_Z[stack_ptr][j][i] = Z[j][i];
1122  }
1123
1124  if (SaveLU) {
1125#ifdef FULL_ALGEBRA
1126    for(i=0; i<NVAR; i++) {
1127      for(j=0; j<NVAR; j++)
1128        chk_J[stack_ptr][i][j] = E[i][j];
1129      chk_P[stack_ptr][i] = P[i];
1130    }
1131#else   
1132    for(i=0; i<LU_NONZERO; i++)
1133      chk_J[stack_ptr][i] = E[i];
1134#endif
1135  }
1136 
1137} /* End of SDIRK_Push */
1138 
1139   
1140/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1141void SDIRK_Pop( KPP_REAL* T, KPP_REAL* H, KPP_REAL* Y, KPP_REAL* Z, KPP_REAL* E,
1142                int* P, int SaveLU ) {
1143/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1144~~~> Retrieves the next trajectory snapshot for discrete adjoints
1145~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1146   
1147  int i,j;
1148
1149  if ( stack_ptr < 0 ) {
1150    printf( "Pop failed: empty buffer\n" );
1151    exit(0);
1152  }
1153
1154  *H = chk_H[ stack_ptr ];
1155  *T = chk_T[ stack_ptr ];
1156  for(i=0; i<NVAR; i++) {
1157    Y[i] = chk_Y[stack_ptr][i];
1158    for(j=0; j<rkS; j++)
1159      Z[(j * NVAR) + i] = chk_Z[stack_ptr][j][i];
1160  }
1161
1162  if (SaveLU) {
1163#ifdef FULL_ALGEBRA
1164    for(i=0; i<NVAR; i++) {
1165      for(j=0; j<NVAR; j++)
1166        E[(j*NVAR)+i] = chk_J[stack_ptr][j][i];
1167      P[i] = chk_P[stack_ptr][i];
1168    }
1169#else
1170    for(i=0; i<LU_NONZERO; i++)
1171      E[i] = chk_J[stack_ptr][i];
1172#endif
1173  }
1174
1175  stack_ptr--;
1176 
1177} /* End of SDIRK_Pop */
1178
1179/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1180void SDIRK_ErrorScale( int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[],
1181                       KPP_REAL Y[], KPP_REAL SCAL[]) {
1182/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1183
1184  int i;
1185  if (ITOL == 0){
1186    for (i = 0; i < NVAR; i++)
1187      SCAL[i] = ONE / (AbsTol[0]+RelTol[0]*ABS(Y[i]) );
1188  }
1189  else {
1190    for (i = 0; i < NVAR; i++)
1191      SCAL[i] = ONE / (AbsTol[i]+RelTol[i]*ABS(Y[i]) );
1192  }
1193
1194} /*  End of SDIRK_ErrorScale */
1195     
1196/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1197KPP_REAL SDIRK_ErrorNorm( int N, KPP_REAL Y[], KPP_REAL SCAL[] ) {
1198/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1199
1200  int i;
1201  KPP_REAL Err = ZERO;
1202
1203  for (i = 0; i < N; i++)
1204    Err = Err + pow( (Y[i]*SCAL[i]), 2);
1205  Err = MAX( SQRT(Err/(KPP_REAL)N), (KPP_REAL)1.0e-10);
1206
1207  return Err;
1208
1209} /*  End of SDIRK_ErrorNorm  */
1210
1211/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1212int SDIRK_ErrorMsg(int code, KPP_REAL T, KPP_REAL H) {
1213/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1214    Handles all error messages
1215~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1216
1217  printf("\nForced exit from Sdirk due to the following error:\n");
1218
1219  switch (code) {
1220    case -1:
1221      printf("--> Improper value for maximal no of steps\n");
1222      break;
1223    case -2:
1224      printf("--> Selected Rosenbrock method not implemented\n");
1225      break;
1226    case -3:
1227      printf("--> Hmin/Hmax/Hstart must be positive\n");
1228      break;
1229    case -4:
1230      printf("--> FacMin/FacMax/FacRej must be positive\n");
1231      break;
1232    case -5:
1233      printf("--> Improper tolerance values\n");
1234      break;
1235    case -6:
1236      printf("--> No of steps exceeds maximum bound\n");
1237      break;
1238    case -7:
1239      printf("--> Step size too small T + 10*H = T or H < Roundoff\n");
1240      break;
1241    case -8:
1242      printf("--> Matrix is repeatedly singular\n");
1243      break;
1244    default: /* causing an error */
1245      printf("Unknown Error code: %d\n", code);
1246  }
1247
1248  printf("\nTime = %f and H = %f\n", T, H );
1249  return code;
1250
1251} /*  end SDIRK_ErrorMsg   */
1252     
1253/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1254void SDIRK_PrepareMatrix ( KPP_REAL H, KPP_REAL T, KPP_REAL Y[], 
1255                           KPP_REAL FJAC[], int SkipJac, int SkipLU, 
1256                           KPP_REAL E[], int IP[], int Reject, int ISING, 
1257                           int ISTATUS[] ) {
1258/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1259~~~>  Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition
1260~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1261
1262  KPP_REAL HGammaInv;
1263  int i, j, ConsecutiveSng = 0;
1264
1265  ISING = 1;
1266
1267  while (ISING != 0) {
1268    HGammaInv = ONE/(H*rkGamma);
1269
1270  /*~~~>  Compute the Jacobian */
1271    if (SkipJac==0) {
1272      JAC_CHEM(T,Y,FJAC);
1273      ISTATUS[Njac]++;
1274    }
1275
1276#ifdef FULL_ALGEBRA
1277    for(j=0; j<NVAR; j++) {
1278      for(i=0; i<NVAR; i++)
1279        E[j][i] = -FJAC[j][i];
1280      E[j][j] = E[j][j] + HGammaInv;
1281    }
1282    DGETRF(NVAR, NVAR, E, NVAR, IP, ISING);
1283#else
1284    for(i=0; i<LU_NONZERO; i++)
1285      E[i] = -FJAC[i];
1286    for(i=0; i<NVAR; i++) {
1287      j = LU_DIAG[i];
1288      E[j]=E[j] + HGammaInv;
1289    }
1290
1291    ISING = KppDecomp(E);
1292    IP[0] = 1;
1293#endif
1294     
1295    ISTATUS[Ndec]++;
1296
1297    if (ISING != 0) {
1298      printf("MATRIX IS SINGULAR, ISING=%d T=%e H=%e\n", ISING, T, H);
1299      ISTATUS[Nsng]++;
1300      ConsecutiveSng++;
1301      if (ConsecutiveSng >= 6) 
1302        return; /* Failure */
1303      H = (KPP_REAL)(0.5)*H;
1304      SkipJac = 0; /* False */
1305      SkipLU  = 0; /* False */
1306      Reject  = 1; /* True */
1307    }
1308  }
1309} /* End of SDIRK_PrepareMatrix */
1310
1311/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1312void SDIRK_Solve ( char Transp, KPP_REAL H, int N, KPP_REAL E[], int IP[], 
1313                   int ISING, KPP_REAL RHS[], int ISTATUS[] ) {
1314/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1315~~~>  Solves the system (H*Gamma-Jac)*x = R
1316      using the LU decomposition of E = I - 1/(H*Gamma)*Jac
1317~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1318
1319  KPP_REAL HGammaInv;
1320     
1321  HGammaInv = ONE/(H*rkGamma);
1322  WSCAL(N,HGammaInv,RHS,1);
1323  switch (Transp) {
1324    case 'N':
1325#ifdef FULL_ALGEBRA 
1326      DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING );
1327#else
1328      KppSolve(E, RHS);
1329#endif
1330      break;
1331    case 'T':
1332#ifdef FULL_ALGEBRA 
1333      DGETRS( 'T', N, 1, E, N, IP, RHS, N, ISING );
1334#else
1335      KppSolveTR(E, RHS, RHS);
1336#endif
1337      break;
1338    default:
1339      printf("Error in SDIRK_Solve. Unknown Transp argument: %c\n", Transp);
1340      exit(0);
1341  }
1342  ISTATUS[Nsol]++;
1343} /* End of SDIRK_Solve */
1344
1345/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1346void Sdirk4a()
1347/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1348{
1349  sdMethod = S4A;
1350
1351/* Number of stages */
1352  rkS = 5;
1353
1354/* Method Coefficients */
1355  rkGamma = (KPP_REAL)0.2666666666666666666666666666666667;
1356
1357  rkA[0][0] = (KPP_REAL)0.2666666666666666666666666666666667;
1358  rkA[0][1] = (KPP_REAL)0.5000000000000000000000000000000000;
1359  rkA[1][1] = (KPP_REAL)0.2666666666666666666666666666666667;
1360  rkA[0][2] = (KPP_REAL)0.3541539528432732316227461858529820;
1361  rkA[1][2] = (KPP_REAL)(-0.5415395284327323162274618585298197e-01);
1362  rkA[2][2] = (KPP_REAL)0.2666666666666666666666666666666667;
1363  rkA[0][3] = (KPP_REAL)0.8515494131138652076337791881433756e-01;
1364  rkA[1][3] = (KPP_REAL)(-0.6484332287891555171683963466229754e-01);
1365  rkA[2][3] = (KPP_REAL)0.7915325296404206392428857585141242e-01;
1366  rkA[3][3] = (KPP_REAL)0.2666666666666666666666666666666667;
1367  rkA[0][4] = (KPP_REAL)2.100115700566932777970612055999074;
1368  rkA[1][4] = (KPP_REAL)(-0.7677800284445976813343102185062276);
1369  rkA[2][4] = (KPP_REAL)2.399816361080026398094746205273880;
1370  rkA[3][4] = (KPP_REAL)(-2.998818699869028161397714709433394);
1371  rkA[4][4] = (KPP_REAL)0.2666666666666666666666666666666667;
1372  rkB[0] = (KPP_REAL)2.100115700566932777970612055999074;
1373  rkB[1] = (KPP_REAL)(-0.7677800284445976813343102185062276);
1374  rkB[2] = (KPP_REAL)2.399816361080026398094746205273880;
1375  rkB[3] = (KPP_REAL)(-2.998818699869028161397714709433394);
1376  rkB[4] = (KPP_REAL)0.2666666666666666666666666666666667;
1377
1378  rkBhat[0] = (KPP_REAL)2.885264204387193942183851612883390;
1379  rkBhat[1] = (KPP_REAL)(-0.1458793482962771337341223443218041);
1380  rkBhat[2] = (KPP_REAL)2.390008682465139866479830743628554;
1381  rkBhat[3] = (KPP_REAL)(-4.129393538556056674929560012190140);
1382  rkBhat[4] = ZERO;
1383
1384  rkC[0] = (KPP_REAL)0.2666666666666666666666666666666667;
1385  rkC[1] = (KPP_REAL)0.7666666666666666666666666666666667;
1386  rkC[2] = (KPP_REAL)0.5666666666666666666666666666666667;
1387  rkC[3] = (KPP_REAL)0.3661315380631796996374935266701191;
1388  rkC[4] = ONE;
1389
1390/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1391  rkD[0] = ZERO;
1392  rkD[1] = ZERO;
1393  rkD[2] = ZERO;
1394  rkD[3] = ZERO;
1395  rkD[4] = ONE;
1396
1397/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1398  rkE[0] = (KPP_REAL)(-0.6804000050475287124787034884002302);
1399  rkE[1] = (KPP_REAL)(1.558961944525217193393931795738823);
1400  rkE[2] = (KPP_REAL)(-13.55893003128907927748632408763868);
1401  rkE[3] = (KPP_REAL)(15.48522576958521253098585004571302);
1402  rkE[4] = ONE;
1403
1404/* Local order of Err estimate */
1405  rkELO = 4;
1406
1407/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1408  rkTheta[0][1] = (KPP_REAL)1.875000000000000000000000000000000;
1409  rkTheta[0][2] = (KPP_REAL)1.708847304091539528432732316227462;
1410  rkTheta[1][2] = (KPP_REAL)(-0.2030773231622746185852981969486824);
1411  rkTheta[0][3] = (KPP_REAL)0.2680325578937783958847157206823118;
1412  rkTheta[1][3] = (KPP_REAL)(-0.1828840955527181631794050728644549);
1413  rkTheta[2][3] = (KPP_REAL)0.2968246986151577397160821594427966;
1414  rkTheta[0][4] = (KPP_REAL)0.9096171815241460655379433581446771;
1415  rkTheta[1][4] = (KPP_REAL)(-3.108254967778352416114774430509465);
1416  rkTheta[2][4] = (KPP_REAL)12.33727431701306195581826123274001;
1417  rkTheta[3][4] = (KPP_REAL)(-11.24557012450885560524143016037523);
1418
1419/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
1420  rkAlpha[0][1] = (KPP_REAL)2.875000000000000000000000000000000;
1421  rkAlpha[0][2] = (KPP_REAL)0.8500000000000000000000000000000000;
1422  rkAlpha[1][2] = (KPP_REAL)0.4434782608695652173913043478260870;
1423  rkAlpha[0][3] = (KPP_REAL)0.7352046091658870564637910527807370;
1424  rkAlpha[1][3] = (KPP_REAL)(-0.9525565003057343527941920657462074e-01);
1425  rkAlpha[2][3] = (KPP_REAL)0.4290111305453813852259481840631738;
1426  rkAlpha[0][4] = (KPP_REAL)(-16.10898993405067684831655675112808);
1427  rkAlpha[1][4] = (KPP_REAL)6.559571569643355712998131800797873;
1428  rkAlpha[2][4] = (KPP_REAL)(-15.90772144271326504260996815012482);
1429  rkAlpha[3][4] = (KPP_REAL)25.34908987169226073668861694892683;
1430
1431  rkELO = (KPP_REAL)4.0;
1432
1433} /* end Sdirk4a */
1434
1435/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1436void Sdirk4b() {
1437/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1438
1439  sdMethod = S4B;
1440
1441/* Number of stages */
1442  rkS = 5;
1443
1444/* Method coefficients */
1445  rkGamma = (KPP_REAL)0.25;
1446
1447  rkA[0][0] = (KPP_REAL)0.25;
1448  rkA[0][1] = (KPP_REAL)0.5;
1449  rkA[1][1] = (KPP_REAL)0.25;
1450  rkA[0][2] = (KPP_REAL)0.34;
1451  rkA[1][2] = (KPP_REAL)(-0.40e-01);
1452  rkA[2][2] = (KPP_REAL)0.25;
1453  rkA[0][3] = (KPP_REAL)0.2727941176470588235294117647058824;
1454  rkA[1][3] = (KPP_REAL)(-0.5036764705882352941176470588235294e-01);
1455  rkA[2][3] = (KPP_REAL)0.2757352941176470588235294117647059e-01;
1456  rkA[3][3] = (KPP_REAL)0.25;
1457  rkA[0][4] = (KPP_REAL)1.041666666666666666666666666666667;
1458  rkA[1][4] = (KPP_REAL)(-1.020833333333333333333333333333333);
1459  rkA[2][4] = (KPP_REAL)7.812500000000000000000000000000000;
1460  rkA[3][4] = (KPP_REAL)(-7.083333333333333333333333333333333);
1461  rkA[4][4] = (KPP_REAL)0.25;
1462
1463  rkB[0] = (KPP_REAL)1.041666666666666666666666666666667;
1464  rkB[1] = (KPP_REAL)(-1.020833333333333333333333333333333);
1465  rkB[2] = (KPP_REAL)7.812500000000000000000000000000000;
1466  rkB[3] = (KPP_REAL)(-7.083333333333333333333333333333333);
1467  rkB[4] = (KPP_REAL)0.250000000000000000000000000000000;
1468
1469  rkBhat[0] = (KPP_REAL)1.069791666666666666666666666666667;
1470  rkBhat[1] = (KPP_REAL)(-0.894270833333333333333333333333333);
1471  rkBhat[2] = (KPP_REAL)7.695312500000000000000000000000000;
1472  rkBhat[3] = (KPP_REAL)(-7.083333333333333333333333333333333);
1473  rkBhat[4] = (KPP_REAL)0.212500000000000000000000000000000;
1474
1475  rkC[0] = (KPP_REAL)0.25;
1476  rkC[1] = (KPP_REAL)0.75;
1477  rkC[2] = (KPP_REAL)0.55;
1478  rkC[3] = (KPP_REAL)0.5;
1479  rkC[4] = ONE;
1480
1481/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1482  rkD[0] = ZERO;
1483  rkD[1] = ZERO;
1484  rkD[2] = ZERO;
1485  rkD[3] = ZERO;
1486  rkD[4] = ONE;
1487
1488/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1489  rkE[0] = (KPP_REAL)0.5750;
1490  rkE[1] = (KPP_REAL)0.2125;
1491  rkE[2] = (KPP_REAL)(-4.6875);
1492  rkE[3] = (KPP_REAL)4.2500;
1493  rkE[4] = (KPP_REAL)0.1500;
1494
1495/* Local order of Err estimate */
1496  rkELO = 4;
1497
1498/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1499  rkTheta[0][1] = (KPP_REAL)2.0;
1500  rkTheta[0][2] = (KPP_REAL)1.680000000000000000000000000000000;
1501  rkTheta[1][2] = (KPP_REAL)(-0.1600000000000000000000000000000000);
1502  rkTheta[0][3] = (KPP_REAL)1.308823529411764705882352941176471;
1503  rkTheta[1][3] = (KPP_REAL)(-0.1838235294117647058823529411764706);
1504  rkTheta[2][3] = (KPP_REAL)0.1102941176470588235294117647058824;
1505  rkTheta[0][4] = (KPP_REAL)(-3.083333333333333333333333333333333);
1506  rkTheta[1][4] = (KPP_REAL)(-4.291666666666666666666666666666667);
1507  rkTheta[2][4] = (KPP_REAL)34.37500000000000000000000000000000;
1508  rkTheta[3][4] = (KPP_REAL)(-28.3333333333333333333333333333);
1509
1510/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
1511  rkAlpha[0][1] = (KPP_REAL)3.0;
1512  rkAlpha[0][2] = (KPP_REAL)0.8800000000000000000000000000000000;
1513  rkAlpha[1][2] = (KPP_REAL)0.4400000000000000000000000000000000;
1514  rkAlpha[0][3] = (KPP_REAL)0.1666666666666666666666666666666667;
1515  rkAlpha[1][3] = (KPP_REAL)(-0.8333333333333333333333333333333333e-01);
1516  rkAlpha[2][3] = (KPP_REAL)0.9469696969696969696969696969696970;
1517  rkAlpha[0][4] = (KPP_REAL)(-6.0);
1518  rkAlpha[1][4] = (KPP_REAL)9.0;
1519  rkAlpha[2][4] = (KPP_REAL)(-56.81818181818181818181818181818182);
1520  rkAlpha[3][4] = (KPP_REAL)54.0;
1521
1522} /* end Sdirk4b */
1523
1524/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1525void Sdirk2a() {
1526/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1527
1528  sdMethod = S2A;
1529
1530/* ~~~> Number of stages */
1531  rkS = 2;
1532
1533/* ~~~> Method coefficients */
1534  rkGamma = (KPP_REAL)0.2928932188134524755991556378951510;
1535  rkA[0][0] = (KPP_REAL)0.2928932188134524755991556378951510;
1536  rkA[0][1] = (KPP_REAL)0.7071067811865475244008443621048490;
1537  rkA[1][1] = (KPP_REAL)0.2928932188134524755991556378951510;
1538  rkB[0] = (KPP_REAL)0.7071067811865475244008443621048490;
1539  rkB[1] = (KPP_REAL)0.2928932188134524755991556378951510;
1540  rkBhat[0] = (KPP_REAL)0.6666666666666666666666666666666667;
1541  rkBhat[1] = (KPP_REAL)0.3333333333333333333333333333333333;
1542  rkC[0] = (KPP_REAL)0.292893218813452475599155637895151;
1543  rkC[1] = ONE;
1544
1545/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1546  rkD[0] = ZERO;
1547  rkD[1] = ONE;
1548
1549/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1550  rkE[0] = (KPP_REAL)0.4714045207910316829338962414032326;
1551  rkE[1] = (KPP_REAL)(-0.1380711874576983496005629080698993);
1552
1553/* ~~~> Local order of Err estimate */
1554  rkELO = 2;
1555
1556/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1557  rkTheta[0][1] = (KPP_REAL)2.414213562373095048801688724209698;
1558
1559/* ~~~> Starting value for Newton iterations */
1560  rkAlpha[0][1] = (KPP_REAL)3.414213562373095048801688724209698;
1561
1562} /* end Sdirk2a */
1563
1564/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1565void Sdirk2b() {
1566/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1567
1568  sdMethod = S2B;
1569
1570/* ~~~> Number of stages */
1571  rkS = 2;
1572
1573/* ~~~> Method coefficients */
1574  rkGamma = (KPP_REAL)1.707106781186547524400844362104849;
1575  rkA[0][0] = (KPP_REAL)1.707106781186547524400844362104849;
1576  rkA[0][1] = (KPP_REAL)(-0.707106781186547524400844362104849);
1577  rkA[1][1] = (KPP_REAL)1.707106781186547524400844362104849;
1578  rkB[0] = (KPP_REAL)(-0.707106781186547524400844362104849);
1579  rkB[1] = (KPP_REAL)1.707106781186547524400844362104849;
1580  rkBhat[0] = (KPP_REAL)0.6666666666666666666666666666666667;
1581  rkBhat[1] = (KPP_REAL)0.3333333333333333333333333333333333;
1582  rkC[0] = (KPP_REAL)1.707106781186547524400844362104849;
1583  rkC[1] = ONE;
1584
1585/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1586  rkD[0] = ZERO;
1587  rkD[1] = ONE;
1588
1589/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1590  rkE[0] = (KPP_REAL)(-0.4714045207910316829338962414032326);
1591  rkE[1] = (KPP_REAL)0.8047378541243650162672295747365659;
1592
1593/* ~~~> Local order of Err estimate */
1594  rkELO = 2;
1595
1596/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1597  rkTheta[0][1] = (KPP_REAL)(-0.414213562373095048801688724209698);
1598
1599/* ~~~> Starting value for Newton iterations */
1600  rkAlpha[0][1] = (KPP_REAL)0.5857864376269049511983112757903019;
1601
1602} /* end Sdirk2b */
1603
1604/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1605void Sdirk3a() {
1606/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1607
1608  sdMethod = S3A;
1609
1610/* ~~~> Number of stages */
1611  rkS = 3;
1612
1613/* ~~~> Method coefficients */
1614  rkGamma = (KPP_REAL)0.2113248654051871177454256097490213;
1615  rkA[0][0] = (KPP_REAL)0.2113248654051871177454256097490213;
1616  rkA[0][1] = (KPP_REAL)0.2113248654051871177454256097490213;
1617  rkA[1][1] = (KPP_REAL)0.2113248654051871177454256097490213;
1618  rkA[0][2] = (KPP_REAL)0.2113248654051871177454256097490213;
1619  rkA[1][2] = (KPP_REAL)0.5773502691896257645091487805019573;
1620  rkA[2][2] = (KPP_REAL)0.2113248654051871177454256097490213;
1621  rkB[0] = (KPP_REAL)0.2113248654051871177454256097490213;
1622  rkB[1] = (KPP_REAL)0.5773502691896257645091487805019573;
1623  rkB[2] = (KPP_REAL)0.2113248654051871177454256097490213;
1624  rkBhat[0]= (KPP_REAL)0.2113248654051871177454256097490213;
1625  rkBhat[1]= (KPP_REAL)0.6477918909913548037576239837516312;
1626  rkBhat[2]= (KPP_REAL)0.1408832436034580784969504064993475;
1627  rkC[0] = (KPP_REAL)0.2113248654051871177454256097490213;
1628  rkC[1] = (KPP_REAL)0.4226497308103742354908512194980427;
1629  rkC[2] = ONE;
1630
1631/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
1632  rkD[0] = ZERO;
1633  rkD[1] = ZERO;
1634  rkD[2] = ONE;
1635
1636/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
1637  rkE[0] = (KPP_REAL)0.9106836025229590978424821138352906;
1638  rkE[1] = (KPP_REAL)(-1.244016935856292431175815447168624);
1639  rkE[2] = (KPP_REAL)0.3333333333333333333333333333333333;
1640
1641/* ~~~> Local order of Err estimate    */
1642  rkELO    = 2;
1643
1644/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
1645  rkTheta[0][1] =  ONE;
1646  rkTheta[0][2] = (KPP_REAL)(-1.732050807568877293527446341505872);
1647  rkTheta[1][2] = (KPP_REAL)2.732050807568877293527446341505872;
1648
1649/* ~~~> Starting value for Newton iterations */
1650  rkAlpha[0][1] = (KPP_REAL)2.0;
1651  rkAlpha[0][2] = (KPP_REAL)(-12.92820323027550917410978536602349);
1652  rkAlpha[1][2] = (KPP_REAL)8.83012701892219323381861585376468;
1653
1654} /* end Sdirk3a */
1655
1656/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1657void FUN_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL P[])
1658/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1659{
1660
1661  KPP_REAL Told;
1662
1663  Told = TIME;
1664  TIME = T;
1665  Update_SUN();
1666  Update_RCONST();
1667  Fun( Y, FIX, RCONST, P );
1668  TIME = Told;
1669
1670}  /*  end FUN_CHEM  */
1671
1672/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1673void JAC_CHEM(KPP_REAL T, KPP_REAL Y[], KPP_REAL JV[]) {
1674/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1675
1676  KPP_REAL Told;
1677
1678#ifdef FULL_ALGEBRA
1679  KPP_REAL JS[LU_NONZERO];
1680  int i,j;
1681#endif
1682
1683  Told = TIME;
1684  TIME = T;
1685  Update_SUN();
1686  Update_RCONST();
1687
1688#ifdef FULL_ALGEBRA
1689  Jac_SP( Y, FIX, RCONST, JS);
1690
1691  for(j=0; j<NVAR; j++) {
1692    for(i=0; i<NVAR; i++)
1693      JV[j][i] = (KPP_REAL)0.0;
1694  } /* end for */
1695
1696  for(i=0; i<LU_NONZERO; i++)
1697    JV[LU_ICOL[i]][LU_IROW[i]] = JS[i];
1698#else
1699  Jac_SP(Y, FIX, RCONST, JV);
1700#endif
1701
1702  TIME = Told;
1703
1704} /* end JAC_CHEM  */
1705/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |