Changeset 25509 in project


Ignore:
Timestamp:
11/16/11 12:10:14 (9 years ago)
Author:
Ivan Raikov
Message:

sundials: added M-L example with IDA solver

File:
1 edited

Legend:

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

    r25082 r25509  
    6666
    6767
    68 
    69 (define (main)
     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))))
    7076 
    71   (let ((yy (f64vector -60.899 0.0149)) ;; v w
     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 tf yy yp ressc 
     88                                     abstol: 1e-14
     89                                     reltol: 1e-14)))
     90
     91      ;; In loop, call IDASolve, print results, and test for error.
     92     
     93      (let recur ((tnext (+ t0 dt)) (iout 1))
     94
     95        (let ((flag  (ida-solve solver tnext)))
     96          (if (negative? flag) (error 'main "IDA solver error" flag))
     97
     98         
     99          (let ((yy (ida-yy solver)))
     100            (printf "~A ~A ~A~%"
     101                    tnext
     102                    (f64vector-ref yy 0)
     103                    (f64vector-ref yy 1)
     104                    ))
     105
     106
     107          (if (< tnext tf)
     108              (recur (+ tnext dt) (+ 1 iout)))
     109          ))
     110     
     111      (ida-destroy-solver solver)
     112     
     113      )))
     114
     115
     116
     117(define (cvode-main)
     118 
     119  (let ((yy (f64vector -60.899 0.0149));; v w
    72120
    73121        ;; Integration limits
     
    89137          (if (negative? flag) (error 'main "CVODE solver error" flag))
    90138
     139
    91140          (if (< tnext tf)
    92141              (recur (+ tnext dt) (+ 1 iout)))
     
    103152      )))
    104153
    105 (define (main/unsafe)
     154(define (cvode-main/unsafe)
    106155 
    107156  (let ((yy (f64vector -60.899 0.0149)) ;; v w
     
    151200            )))
    152201     
    153      
    154 (main)
    155 (main/unsafe)
     202;(ida-main)     
     203(cvode-main)
     204(cvode-main/unsafe)
Note: See TracChangeset for help on using the changeset viewer.