source: branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/sundials/include/cvodes/cvodes_spils.h @ 16222

Last change on this file since 16222 was 16222, checked in by gkronber, 4 years ago

#2925:

  • added comments about parameter identification for differential equation models
  • added source code of cvodes library (part of sundials) which provides functionality to calculate gradients for the parameters of partial differential equation models efficiently using the 'adjoint state method'.
  • added compiled version of cvodes
File size: 24.1 KB
Line 
1/*-----------------------------------------------------------------
2 * Programmer(s): Daniel R. Reynolds @ SMU
3 *                Radu Serban @ LLNL
4 * -----------------------------------------------------------------
5 * LLNS/SMU Copyright Start
6 * Copyright (c) 2017, Southern Methodist University and
7 * Lawrence Livermore National Security
8 *
9 * This work was performed under the auspices of the U.S. Department
10 * of Energy by Southern Methodist University and Lawrence Livermore
11 * National Laboratory under Contract DE-AC52-07NA27344.
12 * Produced at Southern Methodist University and the Lawrence
13 * Livermore National Laboratory.
14 *
15 * All rights reserved.
16 * For details, see the LICENSE file.
17 * LLNS/SMU Copyright End
18 * -----------------------------------------------------------------
19 * This is the header file for the Scaled, Preconditioned Iterative
20 * Linear Solver interface in CVODES.
21 *
22 * Part I contains type definitions and functions for using the
23 * iterative linear solvers on forward problems
24 * (IVP integration and/or FSA)
25 *
26 * Part II contains type definitions and functions for using the
27 * iterative linear solvers on adjoint (backward) problems
28 * -----------------------------------------------------------------*/
29
30#ifndef _CVSSPILS_H
31#define _CVSSPILS_H
32
33#include <sundials/sundials_iterative.h>
34#include <sundials/sundials_nvector.h>
35#include <sundials/sundials_linearsolver.h>
36
37#ifdef __cplusplus  /* wrapper to enable C++ usage */
38extern "C" {
39#endif
40
41/*-----------------------------------------------------------------
42  CVSSPILS return values
43  -----------------------------------------------------------------*/
44
45#define CVSPILS_SUCCESS          0
46#define CVSPILS_MEM_NULL        -1
47#define CVSPILS_LMEM_NULL       -2
48#define CVSPILS_ILL_INPUT       -3
49#define CVSPILS_MEM_FAIL        -4
50#define CVSPILS_PMEM_NULL       -5
51#define CVSPILS_SUNLS_FAIL      -6
52
53/* Return values for the adjoint module */
54
55#define CVSPILS_NO_ADJ          -101
56#define CVSPILS_LMEMB_NULL      -102
57
58/*-----------------------------------------------------------------
59  CVSSPILS solver constants
60  -----------------------------------------------------------------
61  CVSPILS_MSBPRE : maximum number of steps between
62                   preconditioner evaluations
63 
64  CVSPILS_DGMAX  : maximum change in gamma between
65                   preconditioner evaluations
66 
67  CVSPILS_EPLIN  : default value for factor by which the
68                   tolerance on the nonlinear iteration is
69                   multiplied to get a tolerance on the linear
70                   iteration
71  -----------------------------------------------------------------*/
72
73#define CVSPILS_MSBPRE 50
74#define CVSPILS_DGMAX  RCONST(0.2)
75#define CVSPILS_EPLIN  RCONST(0.05)
76
77/*-----------------------------------------------------------------
78  PART I - forward problems
79  -----------------------------------------------------------------*/
80
81/*-----------------------------------------------------------------
82  Type : CVSpilsPrecSetupFn
83  -----------------------------------------------------------------
84  The user-supplied preconditioner setup function PrecSetup and
85  the user-supplied preconditioner solve function PrecSolve
86  together must define left and right preconditoner matrices
87  P1 and P2 (either of which may be trivial), such that the
88  product P1*P2 is an approximation to the Newton matrix
89  M = I - gamma*J.  Here J is the system Jacobian J = df/dy,
90  and gamma is a scalar proportional to the integration step
91  size h.  The solution of systems P z = r, with P = P1 or P2,
92  is to be carried out by the PrecSolve function, and PrecSetup
93  is to do any necessary setup operations.
94 
95  The user-supplied preconditioner setup function PrecSetup
96  is to evaluate and preprocess any Jacobian-related data
97  needed by the preconditioner solve function PrecSolve.
98  This might include forming a crude approximate Jacobian,
99  and performing an LU factorization on the resulting
100  approximation to M.  This function will not be called in
101  advance of every call to PrecSolve, but instead will be called
102  only as often as necessary to achieve convergence within the
103  Inexact Newton iteration.  If the PrecSolve function needs no
104  preparation, the PrecSetup function can be NULL.
105 
106  For greater efficiency, the PrecSetup function may save
107  Jacobian-related data and reuse it, rather than generating it
108  from scratch.  In this case, it should use the input flag jok
109  to decide whether to recompute the data, and set the output
110  flag *jcurPtr accordingly.
111 
112  Each call to the PrecSetup function is preceded by a call to
113  the RhsFn f with the same (t,y) arguments.  Thus the PrecSetup
114  function can use any auxiliary data that is computed and
115  saved by the f function and made accessible to PrecSetup.
116 
117  A function PrecSetup must have the prototype given below.
118  Its parameters are as follows:
119 
120  t       is the current value of the independent variable.
121 
122  y       is the current value of the dependent variable vector,
123           namely the predicted value of y(t).
124 
125  fy      is the vector f(t,y).
126 
127  jok     is an input flag indicating whether Jacobian-related
128          data needs to be recomputed, as follows:
129            jok == SUNFALSE means recompute Jacobian-related data
130                   from scratch.
131            jok == SUNTRUE  means that Jacobian data, if saved from
132                   the previous PrecSetup call, can be reused
133                   (with the current value of gamma).
134          A Precset call with jok == SUNTRUE can only occur after
135          a call with jok == SUNFALSE.
136 
137  jcurPtr is a pointer to an output integer flag which is
138          to be set by PrecSetup as follows:
139          Set *jcurPtr = SUNTRUE if Jacobian data was recomputed.
140          Set *jcurPtr = SUNFALSE if Jacobian data was not recomputed,
141                         but saved data was reused.
142 
143  gamma   is the scalar appearing in the Newton matrix.
144 
145  user_data  is a pointer to user data - the same as the user_data
146          parameter passed to the CVodeSetUserData function.
147 
148  NOTE: If the user's preconditioner needs other quantities,
149        they are accessible as follows: hcur (the current stepsize)
150        and ewt (the error weight vector) are accessible through
151        CVodeGetCurrentStep and CVodeGetErrWeights, respectively).
152        The unit roundoff is available as UNIT_ROUNDOFF defined in
153        sundials_types.h.
154 
155  Returned value:
156  The value to be returned by the PrecSetup function is a flag
157  indicating whether it was successful.  This value should be
158    0   if successful,
159    > 0 for a recoverable error (step will be retried),
160    < 0 for an unrecoverable error (integration is halted).
161  -----------------------------------------------------------------*/
162typedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy,
163          booleantype jok, booleantype *jcurPtr,
164          realtype gamma, void *user_data);
165
166/*-----------------------------------------------------------------
167  Type : CVSpilsPrecSolveFn
168  -----------------------------------------------------------------
169  The user-supplied preconditioner solve function PrecSolve
170  is to solve a linear system P z = r in which the matrix P is
171  one of the preconditioner matrices P1 or P2, depending on the
172  type of preconditioning chosen.
173 
174  A function PrecSolve must have the prototype given below.
175  Its parameters are as follows:
176 
177  t      is the current value of the independent variable.
178 
179  y      is the current value of the dependent variable vector.
180 
181  fy     is the vector f(t,y).
182 
183  r      is the right-hand side vector of the linear system.
184 
185  z      is the output vector computed by PrecSolve.
186 
187  gamma  is the scalar appearing in the Newton matrix.
188 
189  delta  is an input tolerance for use by PSolve if it uses
190         an iterative method in its solution.  In that case,
191         the residual vector Res = r - P z of the system
192         should be made less than delta in weighted L2 norm,
193         i.e., sqrt [ Sum (Res[i]*ewt[i])^2 ] < delta.
194         Note: the error weight vector ewt can be obtained
195         through a call to the routine CVodeGetErrWeights.
196 
197  lr     is an input flag indicating whether PrecSolve is to use
198         the left preconditioner P1 or right preconditioner
199         P2: lr = 1 means use P1, and lr = 2 means use P2.
200 
201  user_data is a pointer to user data - the same as the user_data
202         parameter passed to the CVodeSetUserData function.
203 
204  Returned value:
205  The value to be returned by the PrecSolve function is a flag
206  indicating whether it was successful.  This value should be
207    0 if successful,
208    positive for a recoverable error (step will be retried),
209    negative for an unrecoverable error (integration is halted).
210  -----------------------------------------------------------------*/
211typedef int (*CVSpilsPrecSolveFn)(realtype t, N_Vector y, N_Vector fy,
212          N_Vector r, N_Vector z, realtype gamma,
213          realtype delta, int lr, void *user_data);
214
215/*---------------------------------------------------------------
216  Type: CVSpilsJacTimesSetupFn
217 
218  The user-supplied Jacobian-times-vector product setup function
219  JacTimesSetup and the user-supplied Jacobian-times-vector
220  product function JTimes together must generate the product
221  J*v for v, where J is the Jacobian df/dy, or an approximation
222  to it, and v is a given vector. It should return 0 if
223  successful a positive value for a recoverable error or a
224  negative value for an unrecoverable failure.
225 
226  Each call to the JacTimesSetup function is preceded by a call
227  to the RhsFn fi with the same (t,y) arguments.  Thus the
228  JacTimesSetup function can use any auxiliary data that is
229  computed and saved by the f function and made accessible to
230  JacTimesSetup.
231 
232  A function JacTimesSetup must have the prototype given below.
233  Its parameters are as follows:
234 
235  t       is the current value of the independent variable.
236 
237  y       is the current value of the dependent variable vector,
238           namely the predicted value of y(t).
239 
240  fy      is the vector f(t,y).
241 
242  user_data  is a pointer to user data - the same as the user_data
243          parameter passed to the CVodeSetUserData function.
244 
245  Returned value:
246  The value to be returned by the JacTimesSetup function is a flag
247  indicating whether it was successful.  This value should be
248    0   if successful,
249    > 0 for a recoverable error (step will be retried),
250    < 0 for an unrecoverable error (integration is halted).
251  ---------------------------------------------------------------*/
252typedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y,
253                                      N_Vector fy, void *user_data);
254
255/*-----------------------------------------------------------------
256  Type : CVSpilsJacTimesVecFn
257  -----------------------------------------------------------------
258  The user-supplied function jtimes is to generate the product
259  J*v for given v, where J is the Jacobian df/dy, or an
260  approximation to it, and v is a given vector. It should return
261  0 if successful a positive value for a recoverable error or
262  a negative value for an unrecoverable failure.
263 
264  A function jtimes must have the prototype given below. Its
265  parameters are as follows:
266 
267    v        is the N_Vector to be multiplied by J.
268 
269    Jv       is the output N_Vector containing J*v.
270 
271    t        is the current value of the independent variable.
272 
273    y        is the current value of the dependent variable
274             vector.
275 
276    fy       is the vector f(t,y).
277 
278    user_data   is a pointer to user data, the same as the user_data
279             parameter passed to the CVodeSetUserData function.
280 
281    tmp      is a pointer to memory allocated for an N_Vector
282             which can be used by Jtimes for work space.
283  -----------------------------------------------------------------*/
284typedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
285            N_Vector y, N_Vector fy,
286            void *user_data, N_Vector tmp);
287
288/*=================================================================
289  CVSSPILS Exported functions
290  =================================================================*/
291
292/*-----------------------------------------------------------------
293  Required inputs to the CVSSPILS linear solver interface
294  -----------------------------------------------------------------
295 
296  CVSpilsSetLinearSolver specifies the iterative SUNLinearSolver
297  object that CVode should use.  This is required if CVode is
298  solving a problem with the Newton nonlinear solver (i.e. not the
299  functional iteration).
300 
301  The return value is one of:
302     CVSPILS_SUCCESS   if successful
303     CVSPILS_MEM_NULL  if the CVODE memory was NULL
304     CVSPILS_ILL_INPUT if the linear solver memory was NULL
305 ---------------------------------------------------------------*/
306SUNDIALS_EXPORT int CVSpilsSetLinearSolver(void *cvode_mem,
307                                           SUNLinearSolver LS);
308
309
310/*-----------------------------------------------------------------
311  Optional inputs to the CVSSPILS linear solver
312  -----------------------------------------------------------------
313 
314  CVSpilsSetEpsLin specifies the factor by which the tolerance on
315                 the nonlinear iteration is multiplied to get a
316                 tolerance on the linear iteration.
317                 Default value is 0.05.
318 
319  CVSpilsSetPreconditioner specifies the PrecSetup and PrecSolve
320                 functions.  Default is NULL for both arguments
321                 (no preconditioning)
322 
323  CVSpilsSetJacTimes specifies the jtsetup and jtimes functions.
324                 Default is to use an internal finite difference
325                 approximation routine with no extra jtsetup.
326 
327  The return value of CVSpilsSet* is one of:
328     CVSPILS_SUCCESS   if successful
329     CVSPILS_MEM_NULL  if the cvode memory was NULL
330     CVSPILS_LMEM_NULL if the linear solver memory was NULL
331     CVSPILS_ILL_INPUT if an input has an illegal value
332  -----------------------------------------------------------------*/
333
334SUNDIALS_EXPORT int CVSpilsSetEpsLin(void *cvode_mem, realtype eplifac);
335SUNDIALS_EXPORT int CVSpilsSetPreconditioner(void *cvode_mem,
336                                             CVSpilsPrecSetupFn pset,
337               CVSpilsPrecSolveFn psolve);
338SUNDIALS_EXPORT int CVSpilsSetJacTimes(void *cvode_mem,
339                                       CVSpilsJacTimesSetupFn jtsetup,
340                                       CVSpilsJacTimesVecFn jtimes);
341
342/*-----------------------------------------------------------------
343  Optional outputs from the CVSSPILS linear solver
344  -----------------------------------------------------------------
345  CVSpilsGetWorkSpace returns the real and integer workspace used
346                 by the SPILS module.
347 
348  CVSpilsGetNumPrecEvals returns the number of preconditioner
349                  evaluations, i.e. the number of calls made
350                  to PrecSetup with jok==SUNFALSE.
351 
352  CVSpilsGetNumPrecSolves returns the number of calls made to
353                  PrecSolve.
354 
355  CVSpilsGetNumLinIters returns the number of linear iterations.
356 
357  CVSpilsGetNumConvFails returns the number of linear
358                  convergence failures.
359 
360  CVSpilsGetNumJTSetupEvals returns the number of calls to jtsetup.
361 
362  CVSpilsGetNumJtimesEvals returns the number of calls to jtimes.
363 
364  CVSpilsGetNumRhsEvals returns the number of calls to the user
365                  f routine due to finite difference Jacobian
366                  times vector evaluation.
367 
368  CVSpilsGetLastFlag returns the last error flag set by any of
369                  the CVSPILS interface functions.
370 
371  The return value of CVSpilsGet* is one of:
372     CVSPILS_SUCCESS   if successful
373     CVSPILS_MEM_NULL  if the cvode memory was NULL
374     CVSPILS_LMEM_NULL if the linear solver memory was NULL
375  -----------------------------------------------------------------*/
376
377SUNDIALS_EXPORT int CVSpilsGetWorkSpace(void *cvode_mem,
378                                        long int *lenrwLS,
379                                        long int *leniwLS);
380SUNDIALS_EXPORT int CVSpilsGetNumPrecEvals(void *cvode_mem,
381                                           long int *npevals);
382SUNDIALS_EXPORT int CVSpilsGetNumPrecSolves(void *cvode_mem,
383                                            long int *npsolves);
384SUNDIALS_EXPORT int CVSpilsGetNumLinIters(void *cvode_mem,
385                                          long int *nliters);
386SUNDIALS_EXPORT int CVSpilsGetNumConvFails(void *cvode_mem,
387                                           long int *nlcfails);
388SUNDIALS_EXPORT int CVSpilsGetNumJTSetupEvals(void *cvode_mem,
389                                              long int *njtsetups);
390SUNDIALS_EXPORT int CVSpilsGetNumJtimesEvals(void *cvode_mem,
391                                             long int *njvevals);
392SUNDIALS_EXPORT int CVSpilsGetNumRhsEvals(void *cvode_mem,
393                                          long int *nfevalsLS);
394SUNDIALS_EXPORT int CVSpilsGetLastFlag(void *cvode_mem,
395                                       long int *flag);
396
397/*-----------------------------------------------------------------
398  The following function returns the name of the constant
399  associated with a CVSSPILS return flag
400  -----------------------------------------------------------------*/
401SUNDIALS_EXPORT char *CVSpilsGetReturnFlagName(long int flag);
402
403
404/*-----------------------------------------------------------------
405  PART II - backward problems
406  -----------------------------------------------------------------*/
407
408/*-----------------------------------------------------------------
409  Type : CVSpilsPrecSetupFnB
410  -----------------------------------------------------------------
411  A function PrecSetupB for the adjoint (backward) problem must
412  have the prototype given below.
413  -----------------------------------------------------------------*/
414typedef int (*CVSpilsPrecSetupFnB)(realtype t, N_Vector y, N_Vector yB,
415                                   N_Vector fyB, booleantype jokB,
416                                   booleantype *jcurPtrB,
417                                   realtype gammaB, void *user_dataB);
418
419
420/*----------------------------------------------------------------
421  Type : CVSpilsPrecSetupFnBS
422  -----------------------------------------------------------------
423  A function PrecSetupBS for the adjoint (backward) problem must
424  have the prototype given below.
425  -----------------------------------------------------------------*/
426typedef int (*CVSpilsPrecSetupFnBS)(realtype t, N_Vector y,
427                                    N_Vector *yS, N_Vector yB,
428                                    N_Vector fyB, booleantype jokB,
429                                    booleantype *jcurPtrB,
430                                    realtype gammaB, void *user_dataB);
431
432
433/*-----------------------------------------------------------------
434  Type : CVSpilsPrecSolveFnB
435  -----------------------------------------------------------------
436  A function PrecSolveB for the adjoint (backward) problem  must
437  have the prototype given below.
438  -----------------------------------------------------------------*/
439typedef int (*CVSpilsPrecSolveFnB)(realtype t, N_Vector y, N_Vector yB,
440                                   N_Vector fyB, N_Vector rB,
441                                   N_Vector zB, realtype gammaB,
442                                   realtype deltaB, int lrB,
443                                   void *user_dataB);
444
445/*-----------------------------------------------------------------
446  Type : CVSpilsPrecSolveFnBS
447  -----------------------------------------------------------------
448  A function PrecSolveBS for the adjoint (backward) problem  must
449  have the prototype given below.
450  -----------------------------------------------------------------*/
451typedef int (*CVSpilsPrecSolveFnBS)(realtype t, N_Vector y, N_Vector *yS,
452                                    N_Vector yB, N_Vector fyB,
453                                    N_Vector rB, N_Vector zB,
454                                    realtype gammaB, realtype deltaB,
455                                    int lrB, void *user_dataB);
456
457/*-----------------------------------------------------------------
458  Type : CVSpilsJacTimesSetupFnB
459  -----------------------------------------------------------------
460  A function jtsetupB for the adjoint (backward) problem must have
461  the prototype given below.
462  -----------------------------------------------------------------*/
463typedef int (*CVSpilsJacTimesSetupFnB)(realtype t, N_Vector y, N_Vector yB,
464                                       N_Vector fyB, void *jac_dataB);
465
466/*-----------------------------------------------------------------
467  Type : CVSpilsJacTimesSetupFnBS
468  -----------------------------------------------------------------
469  A function jtsetupBS for the adjoint (backward) problem must have
470  the prototype given below.
471  -----------------------------------------------------------------*/
472typedef int (*CVSpilsJacTimesSetupFnBS)(realtype t, N_Vector y,
473                                        N_Vector *yS, N_Vector yB,
474                                        N_Vector fyB, void *jac_dataB);
475
476/*-----------------------------------------------------------------
477  Type : CVSpilsJacTimesVecFnB
478  -----------------------------------------------------------------
479  A function jtimesB for the adjoint (backward) problem must have
480  the prototype given below.
481  -----------------------------------------------------------------*/
482typedef int (*CVSpilsJacTimesVecFnB)(N_Vector vB, N_Vector JvB, realtype t,
483                                     N_Vector y, N_Vector yB, N_Vector fyB,
484                                     void *jac_dataB, N_Vector tmpB);
485
486/*-----------------------------------------------------------------
487  Type : CVSpilsJacTimesVecFnBS
488  -----------------------------------------------------------------
489  A function jtimesBS for the adjoint (backward) problem must have
490  the prototype given below.
491  -----------------------------------------------------------------*/
492typedef int (*CVSpilsJacTimesVecFnBS)(N_Vector vB, N_Vector JvB,
493                                      realtype t, N_Vector y, N_Vector *yS,
494                                      N_Vector yB, N_Vector fyB,
495                                      void *jac_dataB, N_Vector tmpB);
496
497/*-----------------------------------------------------------------
498  Functions
499  -----------------------------------------------------------------*/
500
501/*---------------------------------------------------------------
502  Required input for the CVSSPILS linear solver interface:
503
504  CVSpilsSetLinearSolverB specifies the iterative SUNLinearSolver
505  object that should be used for the backwards integration.  The
506  'which' argument is the int returned by CVodeCreateB.
507
508  The return value is one of:
509     CVSPILS_SUCCESS   if successful
510     CVSPILS_MEM_NULL  if the cvode memory was NULL
511     CVSPILS_ILL_INPUT if the linear solver memory was NULL
512  ---------------------------------------------------------------*/
513SUNDIALS_EXPORT int CVSpilsSetLinearSolverB(void *cvode_mem,
514                                            int which,
515                                            SUNLinearSolver LS);
516
517/*-----------------------------------------------------------------
518  Each CVSpilsSet***B or CVSpilsSet***BS function below links the
519  main CVODES integrator with the corresponding CVSpilsSet***
520  optional input function for the backward integration.
521  The 'which' argument is the int returned by CVodeCreateB.
522  -----------------------------------------------------------------*/
523
524SUNDIALS_EXPORT int CVSpilsSetEpsLinB(void *cvode_mem, int which, realtype eplifacB);
525
526SUNDIALS_EXPORT int CVSpilsSetPreconditionerB(void *cvode_mem, int which,
527                                              CVSpilsPrecSetupFnB psetB,
528                CVSpilsPrecSolveFnB psolveB);
529SUNDIALS_EXPORT int CVSpilsSetPreconditionerBS(void *cvode_mem, int which,
530                                               CVSpilsPrecSetupFnBS psetBS,
531                 CVSpilsPrecSolveFnBS psolveBS);
532
533SUNDIALS_EXPORT int CVSpilsSetJacTimesB(void *cvode_mem, int which,
534                                        CVSpilsJacTimesSetupFnB jtsetupB,
535                                        CVSpilsJacTimesVecFnB jtimesB);
536SUNDIALS_EXPORT int CVSpilsSetJacTimesBS(void *cvode_mem, int which,
537                                         CVSpilsJacTimesSetupFnBS jtsetupBS,
538                                         CVSpilsJacTimesVecFnBS jtimesBS);
539
540#ifdef __cplusplus
541}
542#endif
543
544#endif
Note: See TracBrowser for help on using the repository browser.