source: project/wiki/eggref/5/sundials @ 40066

Last change on this file since 40066 was 40066, checked in by Ivan Raikov, 7 weeks ago

added C5 sundials doc

File size: 11.7 KB
Line 
1[[tags: eggs]]
2[[toc:]]
3
4== sundials
5
6=== Description
7
8The Chicken {{sundials}} library provides bindings to the solvers from
9the SUNDIALS library. [[http://computation.llnl.gov/casc/sundials/|SUNDIALS]]
10(SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) is a
11collection of solvers for systems of ordinary differential equations
12and differential-algebraic equations.
13
14The Chicken {{sundials}} library provides interfaces to the CVODE and
15IDA solvers and has been tested with SUNDIALS version 3.1.2.
16
17
18=== Library procedures
19
20==== IDA solver interface
21
22
23<procedure>(ida-create-solver TSTART TSTOP VARIABLES DERIVATIVES RESIDUAL-MAIN [RESIDUAL-INIT] [RESIDUAL-EVENT] [EVENTS] [ALG-OR-DIFF] [SUPPRESS] [IC] [USER-DATA] [RELTOL] [ABSTOL]) => IDA-SOLVER</procedure>
24
25Creates and initializes an object representing a problem to be solved
26with the IDA solver. 
27
28Arguments {{TSTART}} and {{TSTOP}} must be real numbers that represent
29the beginning and end of the independent variable range.
30
31Arguments {{VARIABLES}} and {{DERIVATIVES}} must be SRFI-4
32{{f64vector}} objects that hold respectively the initial values and
33derivatives of the system variables.
34
35Argument {{RESIDUAL-MAIN}} is used to compute the residual function
36{{F}} and must be a procedure of the following form:
37
38 (LAMBDA T YY YP DATA)
39
40or
41
42 (LAMBDA T YY YP)
43
44depending on whether the {{USER-DATA}} optional argument is set, where
45
46; {{T}}  :  real-valued independent variable
47; {{YY}} : SRFI-4 {{f64vector}} with current variable values
48; {{YP}} : SRFI-4 {{f64vector}} with current variable derivatives
49; {{DATA}} : is a user data object (if set)
50
51This procedure must return a SRFI-4 {{f64vector}} containing the
52residual vector.
53
54Optional keyword argument {{RESIDUAL-EVENT}} must be a procedure of
55the same form as {{RESIDUAL-MAIN}}, which computes a rootfinding
56problem to be solved during the integration of the system. It is set
57only if argument {{EVENTS}} is given.
58
59Optional keyword argument {{EVENTS}} is an SRFI-4 {{s32vector}} that
60is used for storage of root finding solutions. It must be given if
61{{RESIDUAL-EVENT}} is given.
62
63Optional keyword argument {{ALG-OR-DIFF}} must be an SRFI-4
64{{s32vector}} which indicates the algebraic and differential variables
65in the system. A value of 1 indiciates differential variable, and a
66value of 0 indicates an algebraic one. This is required if the
67{{SUPPRESS}} argument is given and true.
68
69Optional keyword argument {{SUPPRESS}} is a boolean flag that
70indicates whether algebraic variables must be suppressed in the local
71error test. If it is true (suppress), then the argument
72{{ALG-OR-DIFF}} must be given.
73
74Optional keyword argument {{IC}} is a boolean flag that indicates
75whether the solver must calculate consistent initial conditions, or
76whether it must use the initial conditions given by {{VARIABLES}}.
77
78Optional keyword argument {{USER-DATA}} is an object that will be
79passed as an additional argument to the residual functions.
80
81Optional keyword arguments {{RELTOL}} and {{ABSTOL}} specify relative
82and absolute error tolerance, respectively. These both default to
831e-4.
84
85
86<procedure>(ida-reinit-solver IDA-SOLVER T0 Y0 YP0)</procedure>
87
88Re-initializes IDA for the solution of a problem.
89
90<procedure>(ida-destroy-solver IDA-SOLVER)</procedure>
91
92Deallocates the memory associated with the given solver.
93
94<procedure>(ida-solve IDA-SOLVER T)</procedure>
95
96Integrates the system over an interval in the independent
97variable. This procedure returns either when the given {{T}} is
98reached, or when a root is found.
99
100<procedure>(ida-yy IDA-SOLVER)</procedure>
101
102Returns the vector of current state values of the system.
103
104<procedure>(ida-yp IDA-SOLVER)</procedure>
105
106Returns the vector of current state derivative values of the system.
107
108<procedure>(ida-get-last-order IDA-SOLVER)</procedure>
109
110Returns the order used during the last solver step.
111
112<procedure>(ida-get-last-step IDA-SOLVER)</procedure>
113
114Returns the steps size used during the last solver step.
115
116<procedure>(ida-get-num-steps IDA-SOLVER)</procedure>
117
118Returns the cumulative number of steps taken by the solver.
119
120
121==== CVODE solver interface
122
123
124<procedure>(cvode-create-solver TSTART TSTOP VARIABLES RHS-FN [LMM] [ITER] [EWT-FN] [EVENT-FN] [EVENTS] [USER-DATA] [RELTOL] [ABSTOL]) => CVODE-SOLVER</procedure>
125
126Creates and initializes an object representing a problem to be solved
127with the CVODE solver. 
128
129Arguments {{TSTART}} and {{TSTOP}} must be real numbers that represent
130the beginning and end of the independent variable range.
131
132Arguments {{VARIABLES}} must be a SRFI-4 {{f64vector}} object that
133holds the initial values of the system variables.
134
135Argument {{RHS-FN}} is used to compute the right-hand side of the
136equations, and must be a procedure of the following form:
137
138 (LAMBDA T YY DATA)
139
140or
141
142 (LAMBDA T YY)
143
144depending on whether the {{USER-DATA}} optional argument is set, where
145
146; {{T}}  :  real-valued independent variable
147; {{YY}} : SRFI-4 {{f64vector}} with current variable values
148; {{DATA}} : is a user data object (if set)
149
150This procedure must return a SRFI-4 {{f64vector}} containing the
151residual vector.
152
153Optional keyword argument {{EWT-FN}} must be a procedure of the same
154form as {{(LAMBDA YY)}}, which computes error weights for the system
155variables, and which can be used in place of relative and absolute
156error tolerance.
157
158Optional keyword argument {{EVENT-FN}} must be a procedure of
159the same form as {{RHS-FN}}, which computes a rootfinding
160problem to be solved during the integration of the system. It is set
161only if argument {{EVENTS}} is given.
162
163Optional keyword argument {{EVENTS}} is an SRFI-4 {{s32vector}} that
164is used for storage of root finding solutions. It must be given if
165{{EVENT-FN}} is given.
166
167Optional keyword argument {{LMM}} specifies the linear multistep
168method to be used and can be one of {{cvode-lmm/adams}} (default) or
169{{cvode-lmm/bdf}}. {{cvode-lmm/bdf}} is recommended for stiff
170problems.
171
172Optional keyword argument {{ITER}} specifies the iteration type to be
173used and can be one of {{cvode-iter/functional}} (default) or
174{{cvode-iter/newton}}. {{cvode-iter/newton}} is recommended for stiff
175problems.
176
177Optional keyword argument {{USER-DATA}} is an object that will be
178passed as an additional argument to the residual functions.
179
180Optional keyword arguments {{RELTOL}} and {{ABSTOL}} specify relative
181and absolute error tolerance, respectively. These both default to
1821e-4. They are only set of {{EWT-FN}} is not specified.
183
184
185<procedure>(cvode-reinit-solver CVODE-SOLVER T0 Y0 YP0)</procedure>
186
187Re-initializes CVODE for the solution of a problem.
188
189<procedure>(cvode-destroy-solver CVODE-SOLVER)</procedure>
190
191Deallocates the memory associated with the given solver.
192
193<procedure>(cvode-solve CVODE-SOLVER T)</procedure>
194
195Integrates the system over an interval in the independent
196variable. This procedure returns either when the given {{T}} is
197reached, or when a root is found.
198
199<procedure>(cvode-yy CVODE-SOLVER)</procedure>
200
201Returns the vector of current state values of the system.
202
203
204=== Example
205
206 ;;
207 ;; Hodgkin-Huxley model
208 ;;
209 
210 (import mathh sundials srfi-4)
211 
212 (define neg -)
213 (define pow expt)
214 
215 (define TEND  500.0)
216 
217                           
218 ;; Model parameters
219 
220 (define (I_stim t) 10)
221 (define C_m       1)
222 (define E_Na      50)
223 (define E_K       -77)
224 (define E_L       -54.4)
225  (define gbar_Na   120)
226 (define gbar_K    36)
227 (define g_L       0.3)
228 
229 ;; Rate functions
230 
231 (define (amf v)   (* 0.1    (/ (+ v 40)  (- 1.0 (exp (/ (neg (+ v 40)) 10))))))
232 (define (bmf v)   (* 4.0    (exp (/ (neg (+ v 65)) 18))))
233 (define (ahf v)   (* 0.07   (exp (/ (neg (+ v 65)) 20))))
234 (define (bhf v)   (/ 1.0    (+ 1.0 (exp (/ (neg (+ v 35)) 10)))))
235 (define (anf v)   (* 0.01   (/ (+ v 55) (- 1 (exp (/ (neg (+ v 55)) 10))))))
236 (define (bnf v)   (* 0.125  (exp (/ (neg (+ v 65)) 80))))
237 
238 ;; State functions
239 
240 (define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
241 (define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
242 (define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
243                           
244 ;; Model equations
245 
246 (define (rhs t yy)
247 
248   (let ((v (f64vector-ref yy 0))
249        (m (f64vector-ref yy 1))
250        (h (f64vector-ref yy 2))
251        (n (f64vector-ref yy 3)))
252 
253     ;; transition rates at current step
254     (let ((am  (amf v))
255          (an  (anf v))
256          (ah  (ahf v))
257          (bm  (bmf v))
258          (bn  (bnf v))
259          (bh  (bhf v))
260 
261          (g_Na (* gbar_Na  (* h (pow m 3))))
262          (g_K  (* gbar_K   (pow n 4))))
263       
264       (let (
265 
266            ;; currents
267            (I_Na   (* (- v E_Na) g_Na))
268            (I_K    (* (- v E_K)  g_K))
269            (I_L    (* g_L  (- v E_L))))
270                 
271        (let (
272              ;; state equations
273              (dm (- (* am (- 1 m))  (* bm m)))
274              (dh (- (* ah (- 1 h))  (* bh h)))
275              (dn (- (* an (- 1 n))  (* bn n)))
276              (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
277              )
278   
279          (f64vector dv dm dh dn)
280         
281          )))
282     ))
283   
284  (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n
285 
286        ;; Integration limits
287        (t0  0.0)
288        (tf  TEND)
289        (dt  1e-2))
290   
291     ;; CVODE initialization
292     (let ((solver (cvode-create-solver
293                   t0 yy rhs 
294                   tstop: tf
295                   abstol: 1e-4
296                   reltol: 1e-4)))
297 
298       ;; In loop, call CVodeSolve, print results, and test for error.
299       
300       (let recur ((tnext (+ t0 dt)) (iout 1))
301 
302        (let ((flag  (cvode-solve solver tnext)))
303          (if (negative? flag) (error 'main "CVODE solver error" flag))
304 
305          (print-results solver tnext)
306 
307          (if (< tnext tf)
308              (recur (+ tnext dt) (+ 1 iout)))
309          ))
310       
311 
312      (cvode-destroy-solver solver)
313       
314 (define (print-results solver t)
315   (let ((yy (cvode-yy solver)))
316     (printf "~A ~A ~A ~A ~A ~A ~A ~A~%"
317            t
318             (f64vector-ref yy 0)
319             (f64vector-ref yy 1)
320             (f64vector-ref yy 2)
321             (f64vector-ref yy 3)
322             (cvode-get-last-order solver)
323             (cvode-get-num-steps solver)
324             (cvode-get-last-step solver)
325             )))
326
327
328=== Version History
329
330* 2.17 Ported to CHICKEN 5
331* 2.14 More robust SUNDIALS version detection procedure
332* 2.13 Compatibility with argvector CHICKEN
333* 2.12 Corrections to event interface
334* 1.11 unit test update
335* 1.5-1.9 compatibility with versions < 2.4
336* 1.4 tstop made optional; bug fixes in test cases
337* 1.3 Bug fixes in discrete event interface
338* 1.2 Added ida-reinit-solver and cvode-reinit-solver
339* 1.1 Added detection of SUNDIALS package version
340* 1.0 Initial release
341
342=== License
343
344
345  Copyright 2011-2021 Ivan Raikov.
346  All rights reserved.
347 
348  Redistribution and use in source and binary forms, with or without
349  modification, are permitted provided that the following conditions are
350  met:
351 
352  Redistributions of source code must retain the above copyright
353  notice, this list of conditions and the following disclaimer.
354 
355  Redistributions in binary form must reproduce the above copyright
356  notice, this list of conditions and the following disclaimer in the
357  documentation and/or other materials provided with the distribution.
358 
359  Neither the name of the author nor the names of its contributors may
360  be used to endorse or promote products derived from this software
361  without specific prior written permission.
362 
363  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
364  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
365  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
366  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
367  COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
368  INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
369  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
370  SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
371  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
372  STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
373  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
374  OF THE POSSIBILITY OF SUCH DAMAGE.
Note: See TracBrowser for help on using the repository browser.