Changeset 25082 in project


Ignore:
Timestamp:
09/13/11 09:31:05 (10 years ago)
Author:
Ivan Raikov
Message:

sundials: introduced safe and unsafe variants of the solver creation procedures

Location:
release/4/sundials/trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • release/4/sundials/trunk/sundials.meta

    r25057 r25082  
    2020 (needs)
    2121
     22 (test-depends mathh)
     23
    2224 (author "Ivan Raikov")
    2325
  • release/4/sundials/trunk/sundials.scm

    r25079 r25082  
    4444
    4545        (
    46          ida-create-solver ida-destroy-solver ida-solve ida-yy ida-yp
     46         ida-create-solver ida-create-solver/unsafe
     47         ida-destroy-solver ida-solve ida-yy ida-yp
    4748         ida-get-last-order ida-get-num-steps ida-get-last-step
    4849
    4950         cvode-lmm/adams cvode-lmm/bdf
    5051         cvode-iter/functional cvode-iter/newton
    51          cvode-create-solver cvode-destroy-solver cvode-solve cvode-yy
     52         cvode-create-solver cvode-create-solver/unsafe
     53         cvode-destroy-solver cvode-solve cvode-yy
    5254         cvode-get-last-order cvode-get-num-steps cvode-get-last-step
     55
     56         pointer-f64-ref pointer-f64-set! pointer+-f64
    5357         )
    5458
     
    5660        (require-extension srfi-4)
    5761        (require-library lolevel srfi-1 data-structures)
    58         (import (only lolevel object-evict object-release)
     62        (import (only lolevel move-memory! object-evict object-release pointer+ pointer-f64-ref pointer-f64-set!)
    5963                (only srfi-1 fold)
    6064                (only data-structures alist-ref))
     
    643647
    644648
    645 (define (ida-create-solver tstart tstop variables derivatives
    646                            residual-main #!key
    647                            (residual-init #f)
    648                            (residual-event #f)
    649                            (alg-or-diff (make-s32vector (f64vector-length variables) 1))
    650                            (suppress #f)
    651                            (ic #f)
    652                            (user-data #f)
    653                            (events (make-s32vector 0))
    654                            (reltol  1.0e-6)
    655                            (abstol  1.0e-6)
    656                            )
     649(define (ida-create-solver
     650         tstart tstop variables derivatives
     651         residual-main #!key
     652         (residual-init #f)
     653         (residual-event #f)
     654         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
     655         (suppress #f)
     656         (ic #f)
     657         (user-data #f)
     658         (events (make-s32vector 0))
     659         (reltol  1.0e-6)
     660         (abstol  1.0e-6)
     661         )
     662 
     663  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
     664 
     665  (let ((n (f64vector-length variables))
     666        (data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
     667
     668    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
     669
     670    (let ((residual-main1
     671           (if user-data
     672               (lambda (t yy yp rr data) 
     673                 (let ((yy1 (make-f64vector n))
     674                       (yp1 (make-f64vector n)))
     675                   (move-memory! yy yy1 (fx* 8 n))
     676                   (move-memory! yp yp1 (fx* 8 n))
     677                   (let ((v (residual-main t yy1 yp1 data)))
     678                     (move-memory! v rr (fx* 8 n))
     679                     )))
     680               (lambda (t yy yp rr data) 
     681                 (let ((yy1 (make-f64vector n))
     682                       (yp1 (make-f64vector n)))
     683                   (move-memory! yy yy1 (fx* 8 n))
     684                   (move-memory! yp yp1 (fx* 8 n))
     685                   (let ((v (residual-main t yy1 yp1)))
     686                     (move-memory! v rr (* 8 n))
     687                     )))
     688               ))
     689
     690          (residual-init1
     691           (and residual-init
     692                (if user-data
     693                    (lambda (t yy yp rr data) 
     694                      (let ((yy1 (make-f64vector n))
     695                            (yp1 (make-f64vector n)))
     696                        (move-memory! yy yy1 (fx* 8 n))
     697                        (move-memory! yp yp1 (fx* 8 n))
     698                        (let ((v (residual-init t yy1 yp1 data)))
     699                          (move-memory! v rr (fx* 8 n))
     700                          )))
     701                    (lambda (t yy yp rr data) 
     702                      (let ((yy1 (make-f64vector n))
     703                            (yp1 (make-f64vector n)))
     704                        (move-memory! yy yy1 (fx* 8 n))
     705                        (move-memory! yp yp1 (fx* 8 n))
     706                        (let ((v (residual-init t yy1 yp1)))
     707                          (move-memory! v rr (* 8 n))
     708                          )))
     709                    )))
     710
     711          (residual-event1
     712           (and residual-event
     713                (if user-data
     714                    (lambda (t yy yp rr data) 
     715                      (let ((yy1 (make-f64vector n))
     716                            (yp1 (make-f64vector n)))
     717                        (move-memory! yy yy1 (fx* 8 n))
     718                        (move-memory! yp yp1 (fx* 8 n))
     719                        (let ((v (residual-event t yy1 yp1 data)))
     720                          (move-memory! v rr (fx* 8 n))
     721                          )))
     722                    (lambda (t yy yp rr data) 
     723                      (let ((yy1 (make-f64vector n))
     724                            (yp1 (make-f64vector n)))
     725                        (move-memory! yy yy1 (fx* 8 n))
     726                        (move-memory! yp yp1 (fx* 8 n))
     727                        (let ((v (residual-event t yy1 yp1)))
     728                          (move-memory! v rr (* 8 n))
     729                          )))
     730                    )))
     731
     732          )
     733
     734      (ida-residual-main-global residual-main1)
     735      (if residual-init (ida-residual-init-global residual-init1))
     736      (if residual-event (ida-residual-event-global residual-event1)))
     737   
     738    (c-ida-create-solver
     739     tstart tstop 
     740     n (object-evict variables) (object-evict derivatives)
     741     ic alg-or-diff suppress (s32vector-length events) events
     742     data-index abstol reltol)
     743
     744    ))
     745
     746
     747(define (ida-create-solver/unsafe
     748         tstart tstop variables derivatives
     749         residual-main #!key
     750         (residual-init #f)
     751         (residual-event #f)
     752         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
     753         (suppress #f)
     754         (ic #f)
     755         (user-data #f)
     756         (events (make-s32vector 0))
     757         (reltol  1.0e-6)
     758         (abstol  1.0e-6)
     759         )
     760
     761  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
    657762 
    658763  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
     
    664769    (if residual-event (ida-residual-event-global residual-event))
    665770
    666     (assert (= (f64vector-length variables) (f64vector-length derivatives)))
    667771   
    668     (c-ida-create-solver tstart tstop 
    669                          (f64vector-length variables) (object-evict variables) (object-evict derivatives)
    670                          ic alg-or-diff suppress (s32vector-length events) events
    671                          data-index abstol reltol)
     772    (c-ida-create-solver
     773     tstart tstop 
     774     (f64vector-length variables) (object-evict variables) (object-evict derivatives)
     775     ic alg-or-diff suppress (s32vector-length events) events
     776     data-index abstol reltol)
     777
    672778    ))
    673779                       
     
    770876 
    771877
    772 (define (cvode-create-solver tstart tstop variables 
    773                              rhs-fn #!key
    774                              (lmm cvode-lmm/adams)
    775                              (iter cvode-iter/functional)
    776                              (ewt-fn #f)
    777                              (event-fn #f)
    778                              (user-data #f)
    779                              (events (make-s32vector 0))
    780                              (reltol  1.0e-6)
    781                              (abstol  1.0e-6)
    782                              )
     878(define (cvode-create-solver
     879         tstart tstop variables 
     880         rhs-fn #!key
     881         (lmm cvode-lmm/adams)
     882         (iter cvode-iter/functional)
     883         (ewt-fn #f)
     884         (event-fn #f)
     885         (user-data #f)
     886         (events (make-s32vector 0))
     887         (reltol  1.0e-6)
     888         (abstol  1.0e-6)
     889         )
     890
     891  (let ((n (f64vector-length variables))
     892        (data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
     893
     894    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
     895
     896    (let ((rhs-fn1
     897           (if user-data
     898               (lambda (t yy yp data) 
     899                 (let ((yy1 (make-f64vector n)))
     900                   (move-memory! yy yy1 (fx* 8 n))
     901                   (let ((v (rhs-fn t yy1 data)))
     902                     (move-memory! v yp (fx* 8 n))
     903                     )))
     904               (lambda (t yy yp data) 
     905                 (let ((yy1 (make-f64vector n)))
     906                   (move-memory! yy yy1 (fx* 8 n))
     907                   (let ((v (rhs-fn t yy1)))
     908                     (move-memory! v yp (fx* 8 n))
     909                     )))
     910               ))
     911          (ewt-fn1
     912           (if user-data
     913               (lambda (yy ewt data) 
     914                 (let ((yy1 (make-f64vector n)))
     915                   (move-memory! yy yy1 (fx* 8 n))
     916                   (let ((v (ewt-fn yy1 data)))
     917                     (move-memory! v ewt (fx* 8 n))
     918                     )))
     919               (lambda (yy ewt data) 
     920                 (let ((yy1 (make-f64vector n)))
     921                   (move-memory! yy yy1 (fx* 8 n))
     922                   (let ((v (ewt-fn yy1)))
     923                     (move-memory! v ewt (fx* 8 n))
     924                     )))
     925               ))
     926          (event-fn1
     927           (if user-data
     928               (lambda (t yy gout data) 
     929                 (let ((yy1 (make-f64vector n)))
     930                   (move-memory! yy yy1 (fx* 8 n))
     931                   (let ((v (event-fn t yy1 data)))
     932                     (move-memory! v gout (fx* 8 n))
     933                     )))
     934               (lambda (t yy gout data) 
     935                 (let ((yy1 (make-f64vector n)))
     936                   (move-memory! yy yy1 (fx* 8 n))
     937                   (let ((v (event-fn t yy1)))
     938                     (move-memory! v gout (fx* 8 n))
     939                     )))
     940               ))
     941          )
     942   
     943      (cvode-rhs-global rhs-fn1)
     944      (if ewt-fn (cvode-ewt-global ewt-fn1))
     945      (if event-fn (cvode-event-global event-fn1)))
     946   
     947    (c-cvode-create-solver
     948     lmm iter tstart tstop 
     949     n (object-evict variables )
     950     (if ewt-fn 1 0) (s32vector-length events) events
     951     data-index abstol reltol)
     952    ))
     953
     954(define (cvode-create-solver/unsafe
     955         tstart tstop variables 
     956         rhs-fn #!key
     957         (lmm cvode-lmm/adams)
     958         (iter cvode-iter/functional)
     959         (ewt-fn #f)
     960         (event-fn #f)
     961         (user-data #f)
     962         (events (make-s32vector 0))
     963         (reltol  1.0e-6)
     964         (abstol  1.0e-6)
     965         )
    783966 
    784967  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
     
    790973    (if event-fn (cvode-event-global event-fn))
    791974   
    792     (c-cvode-create-solver lmm iter tstart tstop 
    793                            (f64vector-length variables) (object-evict variables )
    794                            (if ewt-fn 1 0) (s32vector-length events) events
    795                            data-index abstol reltol)
     975    (c-cvode-create-solver
     976     lmm iter tstart tstop 
     977     (f64vector-length variables) (object-evict variables )
     978     (if ewt-fn 1 0) (s32vector-length events) events
     979     data-index abstol reltol)
    796980    ))
    797981                       
     
    8611045))
    8621046
     1047
     1048(define (pointer+-f64 p n)
     1049  (pointer+ p (fx* 8 n)))
     1050
     1051
    8631052)
  • release/4/sundials/trunk/tests/e.scm

    r25054 r25082  
    33;;
    44
    5 (use sundials lolevel)
    6 
    7 (define (f64-pointer+ p n)
    8   (pointer+ p (fx* 8 n)))
     5(use sundials srfi-4)
    96
    107;; Problem Constants
     
    1512
    1613
    17 (define (ressc tres yy yp rr _)
     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)
    1820  (let ((v (- (pointer-f64-ref yp) (pointer-f64-ref yy))))
    1921    (pointer-f64-set! rr v)))
    20 
    2122
    2223
     
    5556
    5657
     58(define (main/unsafe)
     59 
     60  (let ((yy (f64vector 1.0))
     61        (yp (f64vector 1.0))
     62
     63         ;; Integration limits
     64         (t0  0.0)
     65         (tf  TEND)
     66         (dt  1e-2))
     67   
     68    ;; IDA initialization
     69    (let ((solver (ida-create-solver/unsafe t0 tf yy yp ressc/unsafe
     70                                            abstol: 1e-14
     71                                            reltol: 1e-14)))
     72
     73      ;; In loop, call IDASolve, print results, and test for error.
     74     
     75      (let recur ((tnext (+ t0 dt)) (iout 1))
     76
     77        (let ((flag  (ida-solve solver tnext)))
     78          (if (negative? flag) (error 'main "IDA solver error" flag))
     79
     80          (if (< tnext tf)
     81              (recur (+ tnext dt) (+ 1 iout)))
     82          ))
     83
     84      (let ((yy (ida-yy solver)))
     85        (assert (< (abs (- 2.71828182846 (f64vector-ref yy 0) )) 1e-12)) )
     86     
     87      (ida-destroy-solver solver)
     88     
     89      )))
     90
     91
    5792(define (print-results solver t)
    5893  (let ((yy (ida-yy solver)))
     
    67102     
    68103(main)
     104(main/unsafe)
  • release/4/sundials/trunk/tests/hh.scm

    r25054 r25082  
    33;;
    44
    5 (use mathh sundials lolevel)
    6 
    7 
    8 
    9 (define (f64-pointer+ p n)
    10   (pointer+ p (fx* 8 n)))
     5(use mathh sundials srfi-4)
    116
    127(define neg -)
     
    4540;; Model equations
    4641
    47 (define (rhs tres yy yp _)
    48   (let ((v (pointer-f64-ref yy))
    49         (m (pointer-f64-ref (f64-pointer+ yy 1)))
    50         (h (pointer-f64-ref (f64-pointer+ yy 2)))
    51         (n (pointer-f64-ref (f64-pointer+ yy 3))))
    52 
     42(define (rhs t yy)
     43
     44  (let ((v (f64vector-ref yy 0))
     45        (m (f64vector-ref yy 1))
     46        (h (f64vector-ref yy 2))
     47        (n (f64vector-ref yy 3)))
    5348
    5449    ;; transition rates at current step
     
    7570              (dh (- (* ah (- 1 h))  (* bh h)))
    7671              (dn (- (* an (- 1 n))  (* bn n)))
    77               (dv (/ (- (I_stim tres) I_L I_Na I_K) C_m))
     72              (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
    7873              )
    7974   
     75          (f64vector dv dm dh dn)
     76         
     77          )))
     78    ))
     79
     80
     81(define (rhs/unsafe t yy yp _)
     82  (let ((v (pointer-f64-ref yy))
     83        (m (pointer-f64-ref (pointer+-f64 yy 1)))
     84        (h (pointer-f64-ref (pointer+-f64 yy 2)))
     85        (n (pointer-f64-ref (pointer+-f64 yy 3))))
     86
     87
     88    ;; transition rates at current step
     89    (let ((am  (amf v))
     90          (an  (anf v))
     91          (ah  (ahf v))
     92          (bm  (bmf v))
     93          (bn  (bnf v))
     94          (bh  (bhf v))
     95
     96          (g_Na (* gbar_Na  (* h (pow m 3))))
     97          (g_K  (* gbar_K   (pow n 4))))
     98     
     99      (let (
     100
     101            ;; currents
     102            (I_Na   (* (- v E_Na) g_Na))
     103            (I_K    (* (- v E_K)  g_K))
     104            (I_L    (* g_L  (- v E_L))))
     105                 
     106        (let (
     107              ;; state equations
     108              (dm (- (* am (- 1 m))  (* bm m)))
     109              (dh (- (* ah (- 1 h))  (* bh h)))
     110              (dn (- (* an (- 1 n))  (* bn n)))
     111              (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
     112              )
     113   
    80114          (pointer-f64-set! yp dv)
    81           (pointer-f64-set! (f64-pointer+ yp 1) dm)
    82           (pointer-f64-set! (f64-pointer+ yp 2) dh)
    83           (pointer-f64-set! (f64-pointer+ yp 3) dn)
     115          (pointer-f64-set! (pointer+-f64 yp 1) dm)
     116          (pointer-f64-set! (pointer+-f64 yp 2) dh)
     117          (pointer-f64-set! (pointer+-f64 yp 3) dn)
    84118         
    85119          )))
     
    98132   
    99133    ;; CVODE initialization
    100     (let ((solver (cvode-create-solver t0 tf yy rhs 
    101                                      abstol: 1e-4
    102                                      reltol: 1e-4)))
     134    (let ((solver (cvode-create-solver
     135                   t0 tf yy rhs 
     136                   abstol: 1e-4
     137                   reltol: 1e-4)))
     138
     139      ;; In loop, call CVodeSolve, print results, and test for error.
     140     
     141      (let recur ((tnext (+ t0 dt)) (iout 1))
     142
     143        (let ((flag  (cvode-solve solver tnext)))
     144          (if (negative? flag) (error 'main "CVODE solver error" flag))
     145
     146          (if (< tnext tf)
     147              (recur (+ tnext dt) (+ 1 iout)))
     148          ))
     149     
     150
     151      (let ((yy (cvode-yy solver)))
     152        (let ((v (f64vector-ref yy 0))
     153              (m (f64vector-ref yy 1))
     154              (h (f64vector-ref yy 2))
     155              (n (f64vector-ref yy 3)))
     156          (assert (< (abs (- v 12.8334007736321)) 1e-12))
     157          (assert (< (abs (- m 0.988277249273334)) 1e-12))
     158          (assert (< (abs (- h 0.151908759594835)) 1e-12))
     159          (assert (< (abs (- n 0.681407916363)) 1e-12))
     160          ))
     161     
     162      (cvode-destroy-solver solver)
     163     
     164      )))
     165
     166(define (main/unsafe)
     167 
     168  (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n
     169
     170        ;; Integration limits
     171        (t0  0.0)
     172        (tf  TEND)
     173        (dt  1e-2))
     174   
     175    ;; CVODE initialization
     176    (let ((solver (cvode-create-solver/unsafe
     177                   t0 tf yy rhs/unsafe 
     178                   abstol: 1e-4
     179                   reltol: 1e-4)))
    103180
    104181      ;; In loop, call CVodeSolve, print results, and test for error.
     
    145222     
    146223(main)
     224(main/unsafe)
  • release/4/sundials/trunk/tests/ml.scm

    r25054 r25082  
    33;;
    44
    5 (use mathh sundials lolevel)
    6 
    7 (define (f64-pointer+ p n)
    8   (pointer+ p (fx* 8 n)))
     5(use mathh sundials srfi-4)
    96
    107
     
    3734;; Model equations
    3835
    39 (define (rhs tres yy yp _)
     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 _)
    4052  (let ((v (pointer-f64-ref yy))
    41         (w (pointer-f64-ref (f64-pointer+ yy 1))))
     53        (w (pointer-f64-ref (pointer+-f64 yy 1))))
    4254
    4355  (let ((ica (* gca (* (minf v)  (- vca v))))
     
    4860   
    4961      (pointer-f64-set! yp dv)
    50       (pointer-f64-set! (f64-pointer+ yp 1) dw)
     62      (pointer-f64-set! (pointer+-f64 yp 1) dw)
    5163
    5264      ))
     
    6577   
    6678    ;; CVODE initialization
    67     (let ((solver (cvode-create-solver t0 tf yy rhs 
    68                                      abstol: 1e-4
    69                                      reltol: 1e-4)))
     79    (let ((solver (cvode-create-solver
     80                   t0 tf yy rhs
     81                   abstol: 1e-4
     82                   reltol: 1e-4)))
     83
     84      ;; In loop, call CVodeSolve, print results, and test for error.
     85
     86      (let recur ((tnext (+ t0 dt)) (iout 1))
     87
     88        (let ((flag  (cvode-solve solver tnext)))
     89          (if (negative? flag) (error 'main "CVODE solver error" flag))
     90
     91          (if (< tnext tf)
     92              (recur (+ tnext dt) (+ 1 iout)))
     93          ))
     94     
     95      (let ((yy (cvode-yy solver)))
     96        (let ((v (f64vector-ref yy 0))
     97              (w (f64vector-ref yy 1)))
     98          (assert (< (abs (- v -30.0536210523504)) 1e-12))
     99          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))
     100
     101      (cvode-destroy-solver solver)
     102     
     103      )))
     104
     105(define (main/unsafe)
     106 
     107  (let ((yy (f64vector -60.899 0.0149)) ;; v w
     108
     109        ;; Integration limits
     110        (t0  0.0)
     111        (tf  TEND)
     112        (dt  1e-2))
     113   
     114    ;; CVODE initialization
     115    (let ((solver (cvode-create-solver/unsafe
     116                   t0 tf yy rhs/unsafe
     117                   abstol: 1e-4
     118                   reltol: 1e-4)))
    70119
    71120      ;; In loop, call CVodeSolve, print results, and test for error.
     
    104153     
    105154(main)
     155(main/unsafe)
Note: See TracChangeset for help on using the changeset viewer.