Changeset 32582 in project


Ignore:
Timestamp:
07/12/15 23:06:37 (5 years ago)
Author:
Ivan Raikov
Message:

sundials: updates to examples with events

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

Legend:

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

    r29469 r32582  
    44;; DIfferential/ALgebraic equation Solvers).
    55;;
    6 ;;  Copyright 2011-2013 Ivan Raikov.
     6;;  Copyright 2011-2015 Ivan Raikov.
    77;;
    88;;
     
    210210
    211211(define-external (ida_residual_event_cb (double t)
    212                                    (c-pointer yy)
    213                                    (c-pointer yp)
    214                                    (c-pointer rr)
    215                                    (unsigned-int data-index)) int
     212                                        (c-pointer yy)
     213                                        (c-pointer yp)
     214                                        (c-pointer rr)
     215                                        (unsigned-int data-index)) int
    216216   (let ((v (and ((dynvector-ref (ida-residual-event-global) data-index)
    217217                  t (c_nvector_data_pointer yy)
     
    652652    }
    653653
    654     if (events > 0)
     654    if (event_number > 0)
    655655    {
    656656#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
     
    841841  (let (
    842842        (n (f64vector-length variables))
    843         (nev (s32vector-length events))
     843        (nevents (s32vector-length events))
    844844        (data-index (dynvector-length (ida-data-global)))
    845845        )
     
    897897                        (move-memory! yp yp1 (fx* 8 n))
    898898                        (let ((v (residual-event t yy1 yp1 data)))
    899                           (move-memory! v rr (fx* 4 nev))
     899                          (move-memory! v rr (fx* 8 nevents))
    900900                          )))
    901901                    (lambda (t yy yp rr data) 
     
    905905                        (move-memory! yp yp1 (fx* 8 n))
    906906                        (let ((v (residual-event t yy1 yp1)))
    907                           (move-memory! v rr (* 4 nev))
     907                          (move-memory! v rr (* 8 nevents))
    908908                          )))
    909909                    )))
     
    10891089  (let (
    10901090        (n (f64vector-length variables))
    1091         (nev (s32vector-length events))
     1091        (nevents (s32vector-length events))
    10921092        (data-index (dynvector-length (cvode-data-global)))
    10931093        )
     
    11321132                        (move-memory! yy yy1 (fx* 8 n))
    11331133                        (let ((v (event-fn t yy1 data)))
    1134                           (move-memory! v gout (fx* 4 nev))
     1134                          (move-memory! v gout (fx* 8 nevents))
    11351135                          )))
    11361136                    (lambda (t yy gout data) 
     
    11381138                        (move-memory! yy yy1 (fx* 8 n))
    11391139                        (let ((v (event-fn t yy1)))
    1140                           (move-memory! v gout (fx* 4 nev))
     1140                          (move-memory! v gout (fx* 8 nevents))
    11411141                          )))
    11421142               )))
  • release/4/sundials/trunk/tests/adex.scm

    r30930 r32582  
    4040               )))
    4141
     42(define (ressc t yy yp)
     43  (let ((dd (subthreshold t yy)))
     44    (let ((v (- (f64vector-ref yp 0) (f64vector-ref dd 0) ))
     45          (w (- (f64vector-ref yp 1) (f64vector-ref dd 1) )))
     46      (f64vector v w))))
     47
     48
    4249(define (spike-detect t yy)
    4350  (let ((V (f64vector-ref yy 0)))
    44     (s32vector (if (>= V theta) 0 1))))
    45  
    46    
     51    (f64vector (- V theta))
     52    ))
     53
     54(define (residual-spike-detect t yy yp)
     55  (spike-detect t yy))
     56
    4757
    4858(define (main)
     
    8696
    8797
     98(define (ida-main)
     99  (let* ((yy (f64vector -70.0 0.0))
     100         (yp (subthreshold 0.0 yy))
     101         
     102         ;; Integration limits
     103         (t0  0.0)
     104         (dt  1e-2))
     105   
     106    ;; IDA initialization
     107    (let ((solver (ida-create-solver t0 yy yp ressc 
     108                                     tstop: tf
     109                                     abstol: 1e-6
     110                                     reltol: 1e-6
     111                                     events: (s32vector 0)
     112                                     residual-event: residual-spike-detect
     113                                     )))
     114
     115      ;; In loop, call IDASolve, print results, and test for error.
     116     
     117      (let recur ((tnext (+ t0 dt)) (iout 1))
     118
     119        (let ((flag  (ida-solve solver tnext)))
     120          (if (negative? flag) (error 'main "IDA solver error" flag))
     121
     122          (print-results/ida solver tnext)
     123
     124          (if (= 1 flag)
     125              (let ((yy (ida-yy solver)))
     126                ;; discrete event: re-initialize solver to integrate over the discontinuity
     127                (ida-reinit-solver solver
     128                                   tnext
     129                                   (f64vector Vr (+ (f64vector-ref yy 1) b))
     130                                   (f64vector 0.0 0.0)
     131                                   )))
     132
     133          (if (< tnext tf)
     134              (recur (+ tnext dt) (+ 1 iout)))
     135          ))
     136     
     137      (ida-destroy-solver solver)
     138     
     139      ))
     140  )
     141
     142
    88143(define (print-results solver)
    89144  (let ((yy (cvode-yy solver)))
     
    97152            )))
    98153
     154
     155(define (print-results/ida solver t)
     156  (let ((yy (ida-yy solver)))
     157    (printf "~A ~A ~A~%"
     158            t
     159            (f64vector-ref yy 0)
     160            (f64vector-ref yy 1)
     161            )))
     162
    99163     
    100 (main)
     164(ida-main)
  • release/4/sundials/trunk/tests/hh.scm

    r30930 r32582  
    203203    (let ((dd (rhs t yy)))
    204204
    205       (let ((dv (- (f64vector-ref dd 0) (f64vector-ref yp 0)))
    206             (dm (- (f64vector-ref dd 1) (f64vector-ref yp 1)))
    207             (dh (- (f64vector-ref dd 2) (f64vector-ref yp 2)))
    208             (dn (- (f64vector-ref dd 3) (f64vector-ref yp 3)))
     205      (let ((dv (- (f64vector-ref yp 0) (f64vector-ref dd 0) ))
     206            (dm (- (f64vector-ref yp 1) (f64vector-ref dd 1) ))
     207            (dh (- (f64vector-ref yp 2) (f64vector-ref dd 2) ))
     208            (dn (- (f64vector-ref yp 3) (f64vector-ref dd 3) ))
    209209            )
    210210        (f64vector dv dm dh dn))))
  • release/4/sundials/trunk/tests/iaf.scm

    r28012 r32582  
    3232(define (spike-detect t yy)
    3333  (let ((v (f64vector-ref yy 0)))
    34     (s32vector (if (>= v theta) 0 1))))
     34    (f64vector (- v theta))
     35    ))
    3536 
    3637   
     
    6364              (cvode-reinit-solver solver (+ tnext trefractory) (f64vector vreset)) )
    6465         
    65 ;;        (print-results solver)
     66          (print-results solver)
    6667
    6768          (if (< tnext tf)
  • release/4/sundials/trunk/tests/izh.scm

    r30965 r32582  
    5454(define (spike-detect t yy)
    5555  (let ((v (f64vector-ref yy 0)))
    56     (s32vector (floor (- v theta)))))
     56    (f64vector (floor (- v theta)))))
    5757             
    5858
     
    102102  (let ((yy (cvode-yy solver))
    103103        (t (cvode-t solver)))
    104     (let ((spike (= (s32vector-ref (spike-detect t yy) 0) 0)))
     104    (let ((spike (>= (f64vector-ref (spike-detect t yy) 0) 0)))
    105105    (printf "~A ~A ~A ~A ~A ~A ~A ~A~%"
    106106            t (or (and spike 1) 0) (or (and spike t) 0.0)
  • release/4/sundials/trunk/tests/izhfs.scm

    r30930 r32582  
    6060(define (spike-detect t yy)
    6161  (let ((v (f64vector-ref yy 0)))
    62     (s32vector (floor (- v Vpeak)))))
     62    (f64vector (floor (- v Vpeak)))))
    6363             
    6464
     
    108108  (let ((yy (cvode-yy solver))
    109109        (t (cvode-t solver)))
    110     (let ((spike (= (s32vector-ref (spike-detect t yy) 0) 0)))
     110    (let ((spike (>= (f64vector-ref (spike-detect t yy) 0) 0)))
    111111    (printf "~A ~A ~A ~A ~A ~A ~A ~A~%"
    112112            t (or (and spike 1) 0) (or (and spike t) 0.0)
  • release/4/sundials/trunk/tests/ml.scm

    r30965 r32582  
    7171    (let ((dd (rhs t yy)))
    7272
    73       (let ((v (- (f64vector-ref dd 0) (f64vector-ref yp 0)))
    74             (w (- (f64vector-ref dd 1) (f64vector-ref yp 1))))
     73      (let ((v (- (f64vector-ref yp 0) (f64vector-ref dd 0) ))
     74            (w (- (f64vector-ref yp 1) (f64vector-ref dd 1) )))
     75
    7576        (f64vector v w))))
    7677 
     
    9798          (if (negative? flag) (error 'main "IDA solver error" flag))
    9899
    99 ;;        (print-results solver tnext)
     100          (print-results/ida solver tnext)
    100101
    101102          (if (< tnext tf)
     
    132133          (if (negative? flag) (error 'main "CVODE solver error" flag))
    133134
    134           (print-results solver tnext)
     135          (print-results/cvode solver tnext)
    135136
    136137          (if (< tnext tf)
     
    171172          (if (negative? flag) (error 'main "CVODE solver error" flag))
    172173
    173 ;         (print-results solver tnext)
     174;         (print-results/cvode solver tnext)
    174175
    175176          (if (< tnext tf)
     
    188189
    189190
    190 (define (print-results solver t)
     191(define (print-results/cvode solver t)
    191192  (let ((yy (cvode-yy solver)))
    192193    (printf "~A ~A ~A ~A ~A ~A~%"
     
    198199            (cvode-get-last-step solver)
    199200            )))
     201
     202(define (print-results/ida solver t)
     203  (let ((yy (ida-yy solver)))
     204    (printf "~A ~A ~A~%"
     205            t
     206            (f64vector-ref yy 0)
     207            (f64vector-ref yy 1)
     208            )))
    200209     
    201210(cvode-main)
    202 ;(cvode-main/unsafe)
    203 ;(ida-main)     
     211(cvode-main/unsafe)
     212(ida-main)     
Note: See TracChangeset for help on using the changeset viewer.