source: project/release/4/sundials/trunk/tests/ml.scm @ 25539

Last change on this file since 25539 was 25539, checked in by Ivan Raikov, 9 years ago

sundials: fixes to unit tests

File size: 4.4 KB
Line 
1;;
2;; Morris-Lecar model
3;;
4
5(use mathh sundials srfi-4)
6
7
8(define TEND  500.0)
9
10                           
11;; Model parameters
12
13(define Istim  20)
14(define vk     -84)
15(define vl     -60)
16(define vca    120)
17(define gk     2.5)
18(define gl     0.5)
19(define gca    1)
20(define c      2)
21(define v1     -1.2)
22(define v2     18)
23(define v3     12)
24(define v4     17)
25(define phi    0.0667)
26
27
28;; State functions
29
30(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
31(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
32(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
33                           
34;; Model equations
35
36(define (rhs t yy)
37  (let ((v (f64vector-ref yy 0))
38        (w (f64vector-ref yy 1)))
39
40  (let ((ica (* gca (* (minf v)  (- vca v))))
41        (ik  (* gk  (* w (- vk v )))))
42   
43    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
44          (dw (* (lamw v) (- (winf v) w))))
45   
46      (f64vector dv dw)
47
48      ))
49  ))
50
51(define (rhs/unsafe t yy yp _)
52  (let ((v (pointer-f64-ref yy))
53        (w (pointer-f64-ref (pointer+-f64 yy 1))))
54
55  (let ((ica (* gca (* (minf v)  (- vca v))))
56        (ik  (* gk  (* w (- vk v )))))
57   
58    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
59          (dw (* (lamw v) (- (winf v) w))))
60   
61      (pointer-f64-set! yp dv)
62      (pointer-f64-set! (pointer+-f64 yp 1) dw)
63
64      ))
65  ))
66
67
68(define (ida-main)
69
70  (define (ressc t yy yp)
71    (let ((dd (rhs t yy)))
72
73      (let ((v (- (f64vector-ref dd 0) (f64vector-ref yp 0)))
74            (w (- (f64vector-ref dd 1) (f64vector-ref yp 1))))
75        (f64vector v w))))
76 
77  (let* ((yy (f64vector -60.899 0.0149)) ;; v w
78         (yp (rhs 0.0 yy))
79         
80         ;; Integration limits
81         (t0  0.0)
82         (tf  TEND)
83         (dt  1e-2))
84   
85
86    ;; IDA initialization
87    (let ((solver (ida-create-solver t0 yy yp ressc 
88                                     tstop: tf
89                                     abstol: 1e-14
90                                     reltol: 1e-14)))
91
92      ;; In loop, call IDASolve, print results, and test for error.
93     
94      (let recur ((tnext (+ t0 dt)) (iout 1))
95
96        (let ((flag  (ida-solve solver tnext)))
97          (if (negative? flag) (error 'main "IDA solver error" flag))
98
99          ;;(print-results solver tnext)
100         
101          (let ((yy (ida-yy solver)))
102            (printf "~A ~A ~A~%" 
103                    tnext
104                    (f64vector-ref yy 0)
105                    (f64vector-ref yy 1)
106                    ))
107
108
109          (if (< tnext tf)
110              (recur (+ tnext dt) (+ 1 iout)))
111          ))
112     
113      (ida-destroy-solver solver)
114     
115      )))
116
117
118
119(define (cvode-main)
120 
121  (let ((yy (f64vector -60.899 0.0149));; v w
122
123        ;; Integration limits
124        (t0  0.0)
125        (tf  TEND)
126        (dt  1e-2))
127   
128    ;; CVODE initialization
129    (let ((solver (cvode-create-solver
130                   t0 yy rhs
131                   tstop: tf
132                   abstol: 1e-4
133                   reltol: 1e-4)))
134
135      ;; In loop, call CVodeSolve, print results, and test for error.
136
137      (let recur ((tnext (+ t0 dt)) (iout 1))
138
139        (let ((flag  (cvode-solve solver tnext)))
140          (if (negative? flag) (error 'main "CVODE solver error" flag))
141
142          ;;(print-results solver tnext)
143
144          (if (< tnext tf)
145              (recur (+ tnext dt) (+ 1 iout)))
146          ))
147     
148      (let ((yy (cvode-yy solver)))
149        (let ((v (f64vector-ref yy 0))
150              (w (f64vector-ref yy 1)))
151          (assert (< (abs (- v -30.0333833471119 )) 1e-12))
152          (assert (< (abs (- w 0.04023343585461 )) 1e-12))
153          ))
154
155      (cvode-destroy-solver solver)
156     
157      )))
158
159(define (cvode-main/unsafe)
160 
161  (let ((yy (f64vector -60.899 0.0149)) ;; v w
162
163        ;; Integration limits
164        (t0  0.0)
165        (tf  TEND)
166        (dt  1e-2))
167   
168    ;; CVODE initialization
169    (let ((solver (cvode-create-solver/unsafe
170                   t0 yy rhs/unsafe
171                   tstop: tf
172                   abstol: 1e-4
173                   reltol: 1e-4)))
174
175      ;; In loop, call CVodeSolve, print results, and test for error.
176
177      (let recur ((tnext (+ t0 dt)) (iout 1))
178
179        (let ((flag  (cvode-solve solver tnext)))
180          (if (negative? flag) (error 'main "CVODE solver error" flag))
181
182          ;;(print-results solver tnext)
183
184          (if (< tnext tf)
185              (recur (+ tnext dt) (+ 1 iout)))
186          ))
187     
188      (let ((yy (cvode-yy solver)))
189        (let ((v (f64vector-ref yy 0))
190              (w (f64vector-ref yy 1)))
191          (assert (< (abs (- v -30.0333833471119 )) 1e-12))
192          (assert (< (abs (- w 0.04023343585461 )) 1e-12))
193          ))
194
195      (cvode-destroy-solver solver)
196     
197      )))
198
199
200(define (print-results solver t)
201  (let ((yy (cvode-yy solver)))
202    (printf "~A ~A ~A ~A ~A ~A~%" 
203            t 
204            (f64vector-ref yy 0)
205            (f64vector-ref yy 1)
206            (cvode-get-last-order solver)
207            (cvode-get-num-steps solver)
208            (cvode-get-last-step solver)
209            )))
210     
211;(ida-main)     
212(cvode-main)
213(cvode-main/unsafe)
Note: See TracBrowser for help on using the repository browser.