source: project/release/4/sundials/tags/1.4/tests/e.scm @ 25541

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

sundials release 1.4

File size: 2.1 KB
Line 
1;;
2;; Compute Euler's number
3;;
4
5(use sundials srfi-4)
6
7;; Problem Constants
8
9(define NEQ 1)
10
11(define TEND  1.0)
12
13
14(define (ressc t yy yp)
15  (let ((v (- (f64vector-ref yp 0) (f64vector-ref yy 0))))
16    (f64vector v)))
17
18
19(define (ressc/unsafe t yy yp rr data)
20  (let ((v (- (pointer-f64-ref yp) (pointer-f64-ref yy))))
21    (pointer-f64-set! rr v)))
22
23
24(define (main)
25 
26  (let ((yy (f64vector 1.0))
27        (yp (f64vector 1.0))
28
29         ;; Integration limits
30         (t0  0.0)
31         (tf  TEND)
32         (dt  1e-2))
33   
34    ;; IDA initialization
35    (let ((solver (ida-create-solver t0 yy yp ressc 
36                                     tstop: tf
37                                     abstol: 1e-14
38                                     reltol: 1e-14)))
39
40      ;; In loop, call IDASolve, print results, and test for error.
41     
42      (let recur ((tnext (+ t0 dt)) (iout 1))
43
44        (let ((flag  (ida-solve solver tnext)))
45          (if (negative? flag) (error 'main "IDA solver error" flag))
46
47          (if (< tnext tf)
48              (recur (+ tnext dt) (+ 1 iout)))
49          ))
50
51      (let ((yy (ida-yy solver)))
52        (assert (< (abs (- 2.71828182846 (f64vector-ref yy 0) )) 1e-12)) )
53     
54      (ida-destroy-solver solver)
55     
56      )))
57
58
59(define (main/unsafe)
60 
61  (let ((yy (f64vector 1.0))
62        (yp (f64vector 1.0))
63
64         ;; Integration limits
65         (t0  0.0)
66         (tf  TEND)
67         (dt  1e-2))
68   
69    ;; IDA initialization
70    (let ((solver (ida-create-solver/unsafe t0 yy yp ressc/unsafe
71                                            tstop: tf
72                                            abstol: 1e-14
73                                            reltol: 1e-14)))
74
75      ;; In loop, call IDASolve, print results, and test for error.
76     
77      (let recur ((tnext (+ t0 dt)) (iout 1))
78
79        (let ((flag  (ida-solve solver tnext)))
80          (if (negative? flag) (error 'main "IDA solver error" flag))
81
82          (if (< tnext tf)
83              (recur (+ tnext dt) (+ 1 iout)))
84          ))
85
86      (let ((yy (ida-yy solver)))
87        (assert (< (abs (- 2.71828182846 (f64vector-ref yy 0) )) 1e-12)) )
88     
89      (ida-destroy-solver solver)
90     
91      )))
92
93
94(define (print-results solver t)
95  (let ((yy (ida-yy solver)))
96    (printf "~A ~A ~A ~A~%" 
97            t 
98            (f64vector-ref yy 0)
99            (ida-get-last-order solver)
100            (ida-get-num-steps solver)
101            (ida-get-last-step solver)
102            )))
103     
104     
105(main)
106(main/unsafe)
Note: See TracBrowser for help on using the repository browser.