source: project/release/4/sundials/trunk/sundials.scm @ 25245

Last change on this file since 25245 was 25245, checked in by Ivan Raikov, 10 years ago

sundials: added a test for discrete events and fixed some bugs in setting up discrete event routine

File size: 30.1 KB
Line 
1;;
2;;
3;; Chicken Scheme bindings for SUNDIALS (SUite of Nonlinear and
4;; DIfferential/ALgebraic equation Solvers).
5;;
6;;  Copyright 2011 Ivan Raikov.
7;;
8;;
9;;  Redistribution and use in source and binary forms, with or without
10;;  modification, are permitted provided that the following conditions
11;;  are met:
12;;
13;;  - Redistributions of source code must retain the above copyright
14;;  notice, this list of conditions and the following disclaimer.
15;;
16;;  - Redistributions in binary form must reproduce the above
17;;  copyright notice, this list of conditions and the following
18;;  disclaimer in the documentation and/or other materials provided
19;;  with the distribution.
20;;
21;;  - Neither name of the copyright holders nor the names of its
22;;  contributors may be used to endorse or promote products derived
23;;  from this software without specific prior written permission.
24;;
25;;  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND THE
26;;  CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
27;;  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
28;;  MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
29;;  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR THE
30;;  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31;;  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32;;  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
33;;  USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
34;;  AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
35;;  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
36;;  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
37;;  POSSIBILITY OF SUCH DAMAGE.
38;;
39;; 
40;;
41
42
43(module sundials
44
45        (
46         ida-create-solver ida-create-solver/unsafe
47         ida-reinit-solver ida-destroy-solver ida-solve ida-yy ida-yp
48         ida-get-last-order ida-get-num-steps ida-get-last-step
49
50         cvode-lmm/adams cvode-lmm/bdf
51         cvode-iter/functional cvode-iter/newton
52         cvode-create-solver cvode-create-solver/unsafe
53         cvode-reinit-solver cvode-destroy-solver cvode-solve cvode-yy
54         cvode-get-last-order cvode-get-num-steps cvode-get-last-step
55
56         pointer-f64-ref pointer-f64-set! pointer+-f64
57         )
58
59        (import scheme chicken foreign)
60        (require-extension srfi-4)
61        (require-library lolevel srfi-1 data-structures)
62        (import (only lolevel move-memory! object-evict object-release pointer+ pointer-f64-ref pointer-f64-set!)
63                (only srfi-1 fold)
64                (only data-structures alist-ref))
65
66
67
68#>
69
70#include <stdlib.h>
71#include <math.h>
72#include <assert.h>
73
74#include <ida/ida.h>
75#include <ida/ida_dense.h>
76#include <cvode/cvode.h>           
77#include <cvode/cvode_band.h>     
78#include <nvector/nvector_serial.h>
79
80
81static void chicken_panic (C_char *) C_noret;
82static void chicken_panic (C_char *msg)
83{
84  C_word *a = C_alloc (C_SIZEOF_STRING (strlen (msg)));
85  C_word scmmsg = C_string2 (&a, msg);
86  C_halt (scmmsg);
87  exit (5); /* should never get here */
88}
89
90static void chicken_throw_exception(C_word value) C_noret;
91static void chicken_throw_exception(C_word value)
92{
93  char *aborthook = C_text("\003sysabort");
94
95  C_word *a = C_alloc(C_SIZEOF_STRING(strlen(aborthook)));
96  C_word abort = C_intern2(&a, aborthook);
97
98  abort = C_block_item(abort, 0);
99  if (C_immediatep(abort))
100    chicken_panic(C_text("`##sys#abort' is not defined"));
101
102  C_save(value);
103  C_do_apply(1, abort, C_SCHEME_UNDEFINED);
104}
105
106void chicken_error (char *msg, C_word obj) 
107{
108  size_t msglen;
109  C_word *a;
110  C_word scmmsg;
111  C_word list;
112
113  msglen = strlen (msg);
114  a = C_alloc (C_SIZEOF_STRING (msglen) + C_SIZEOF_LIST(2));
115  scmmsg = C_string2 (&a, (char *) msg);
116  list = C_list(&a, 2, scmmsg, obj);
117  chicken_throw_exception(list);
118}
119
120
121typedef struct IDA_Solver_Handle_struct
122{
123    void* ida_mem;
124    double tret;
125    C_word syy;
126    C_word syp;
127    N_Vector yy;
128    N_Vector yp;
129    N_Vector id;
130    int event_number;
131    int* events;
132    double abstol;
133    double reltol;
134
135} IDA_Solver_Handle;
136
137
138typedef struct CVODE_Solver_Handle_struct
139{
140    void* cvode_mem;
141    double tret;
142    C_word syy;
143    N_Vector yy;
144    N_Vector id;
145    int event_number;
146    int* events;
147    double abstol;
148    double reltol;
149
150} CVODE_Solver_Handle;
151
152<#
153
154(define-foreign-type IDASolverHandle "IDA_Solver_Handle")
155(define-foreign-type CVODESolverHandle "CVODE_Solver_Handle")
156
157
158(define cvode-lmm/adams (foreign-value "CV_ADAMS" int))
159(define cvode-lmm/bdf (foreign-value "CV_BDF" int))
160
161(define cvode-iter/functional (foreign-value "CV_FUNCTIONAL" int))
162(define cvode-iter/newton (foreign-value "CV_NEWTON" int))
163
164(define ida-residual-init-global  (make-parameter (lambda _ (begin 0))))
165(define ida-residual-main-global  (make-parameter (lambda _ (begin 0))))
166(define ida-residual-event-global (make-parameter (lambda _ (begin 0))))
167
168(define ida-data-global (make-parameter '()))
169
170
171(define cvode-rhs-global   (make-parameter (lambda _ (begin 0))))
172(define cvode-event-global  (make-parameter (lambda _ (begin 0))))
173(define cvode-ewt-global   (make-parameter (lambda _ (begin 0))))
174
175(define cvode-data-global (make-parameter '()))
176
177
178
179(define c_nvector_data_pointer 
180  (foreign-safe-lambda* nonnull-c-pointer ((nonnull-c-pointer x))
181#<<EOF
182   N_Vector v = (N_Vector)x;
183   C_return (NV_DATA_S(v));
184EOF
185))
186
187
188(define-external (ida_residual_main_cb (double t) 
189                                   (c-pointer yy) 
190                                   (c-pointer yp) 
191                                   (c-pointer rr)
192                                   (unsigned-int data-index)) int 
193   (let ((v (and ((ida-residual-main-global) t (c_nvector_data_pointer yy) 
194                  (c_nvector_data_pointer yp) (c_nvector_data_pointer rr )
195                  (alist-ref data-index (ida-data-global) )) 0)))
196     v))
197
198
199(define-external (ida_residual_init_cb (double t) 
200                                   (c-pointer yy) 
201                                   (c-pointer yp) 
202                                   (c-pointer rr)
203                                   (unsigned-int data-index)) int 
204   (let ((v ((ida-residual-init-global) t (c_nvector_data_pointer yy) 
205             (c_nvector_data_pointer yp) (c_nvector_data_pointer rr )
206             (alist-ref data-index (ida-data-global) ))))
207     v))
208
209
210
211(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 
216   (let ((v (and ((ida-residual-event-global) t (c_nvector_data_pointer yy) 
217                  (c_nvector_data_pointer yp) rr
218                  (alist-ref data-index (ida-data-global) )) 0)))
219     v))
220
221
222(define-external (cvode_rhs_cb (double t) 
223                               (c-pointer yy) 
224                               (c-pointer yp) 
225                               (unsigned-int data-index)) int 
226   (let ((v (and ((cvode-rhs-global) t (c_nvector_data_pointer yy) (c_nvector_data_pointer yp) 
227                  (alist-ref data-index (cvode-data-global) )) 0)))
228     v))
229
230
231(define-external (cvode_event_cb (double t) 
232                                 (c-pointer yy) 
233                                 (c-pointer gout) 
234                                 (unsigned-int data-index)) int 
235   (let ((v (and ((cvode-event-global) t (c_nvector_data_pointer yy) gout
236                  (alist-ref data-index (cvode-data-global) )) 0)))
237     v))
238
239
240(define-external (cvode_ewt_cb (c-pointer yy) 
241                               (c-pointer ewt) 
242                               (unsigned-int data-index)) int 
243   (let ((v (and ((cvode-ewt-global) (c_nvector_data_pointer yy) (c_nvector_data_pointer ewt) 
244                  (alist-ref data-index (cvode-data-global) )) 0)))
245     v))
246
247
248#>
249
250C_externexport  int  ida_residual_init_cb(double,void *,void *,void *,unsigned int);
251C_externexport  int  ida_residual_main_cb(double,void *,void *,void *,unsigned int);
252C_externexport  int  ida_residual_event_cb(double,void *,void *,void *,unsigned int);
253
254C_externexport  int  cvode_rhs_cb(double,void *,void *,unsigned int);
255C_externexport  int  cvode_event_cb(double,void *,void *,unsigned int);
256C_externexport  int  cvode_ewt_cb(void *,void *,unsigned int);
257
258
259void adjust_zero_crossings (N_Vector v, double abstol)
260{
261    int i;
262    for (i = 0; i < NV_LENGTH_S(v); i++)
263    {
264        if (fabs(NV_Ith_S(v,i)) < abstol) NV_Ith_S(v,i) = 0;
265    }
266    return;
267}
268
269
270IDA_Solver_Handle* ida_create_solver
271( double time_start,
272  double time_stop,
273
274  int variable_number,
275  C_word variables, 
276  C_word derivatives, 
277  int ic,
278  int* alg_or_diff,
279  int suppress,
280
281  int event_number, 
282  int* events, 
283
284  unsigned int data_index,
285
286  double abstol,
287  double reltol
288)
289{
290    IDA_Solver_Handle* solver_handle;
291    assert ((solver_handle = malloc (sizeof(struct IDA_Solver_Handle_struct))) != NULL);
292
293    solver_handle->tret = 0.0;
294    solver_handle->event_number = event_number;
295    solver_handle->events = events;
296
297    int flag = 0;
298    int i = 0;
299
300    solver_handle->syy = variables;
301    solver_handle->syp = derivatives;
302
303    solver_handle->yy = N_VMake_Serial(variable_number, C_c_f64vector(variables));
304    solver_handle->yp = N_VMake_Serial(variable_number, C_c_f64vector(derivatives));
305
306    solver_handle->id = N_VNew_Serial(variable_number);
307
308    for (i = 0; i < variable_number; i++)
309    {
310        NV_Ith_S(solver_handle->id,i) = (alg_or_diff[i] == 0) ? 0.0 : 1.0;
311    }
312
313    solver_handle->abstol = abstol;
314    solver_handle->reltol = reltol;
315
316    solver_handle->ida_mem = IDACreate();
317    if (solver_handle->ida_mem == NULL) 
318    {
319       chicken_error("could not allocate memory with IDACreate", C_SCHEME_UNDEFINED);
320    }
321
322    flag = IDASetUserData(solver_handle->ida_mem, (void *)data_index);
323    if (flag != IDA_SUCCESS) 
324    {
325       chicken_error("could not set user data with IDASetUserData", C_fix(flag));
326    }
327
328    flag = IDAInit (solver_handle->ida_mem, 
329                    (IDAResFn)ida_residual_main_cb, 
330                    time_start, 
331                    solver_handle->yy, 
332                    solver_handle->yp);
333    if (flag != IDA_SUCCESS) 
334        chicken_error("could not initialize solver with IDAInit", C_fix(flag)) ;       
335
336    flag = IDASStolerances(solver_handle->ida_mem, solver_handle->reltol, solver_handle->abstol);
337    if (flag != IDA_SUCCESS) 
338        chicken_error("could not set tolerances with IDASStolerances", C_fix(flag)) ;   
339
340    flag = IDASetId(solver_handle->ida_mem, solver_handle->id);
341    if (flag != IDA_SUCCESS) 
342        chicken_error("could not set algebraic variables with IDASetId", C_fix(flag)) ;
343       
344    flag = IDASetSuppressAlg(solver_handle->ida_mem, suppress);
345    if (flag != IDA_SUCCESS) 
346        chicken_error("could not set error suppression flag with IDASetSuppressAlg", C_fix(flag)) ;     
347
348    flag = IDADense(solver_handle->ida_mem, variable_number);
349#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
350    if (flag != IDADLS_SUCCESS) 
351#else
352    if (flag != IDADENSE_SUCCESS) 
353#endif
354       chicken_error("could not initialize linear solver with IDADense", C_fix(flag)) ;
355 
356    if (ic)
357    {
358       flag = IDACalcIC(solver_handle->ida_mem, IDA_YA_YDP_INIT, time_start + solver_handle->abstol);
359      if (flag != IDA_SUCCESS) 
360        chicken_error("could not compute initial conditions with IDACalcIC", C_fix(flag)) ;     
361    }
362
363    if (events > 0)
364    { 
365#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
366       flag = IDARootInit(solver_handle->ida_mem, event_number, (IDARootFn)ida_residual_event_cb);
367#else
368       flag = IDARootInit(solver_handle->ida_mem, event_number, (IDARootFn)ida_residual_event_cb, (void *)data_index);
369#endif
370       if (flag != IDA_SUCCESS) 
371          chicken_error("could not initialize event variables with IDARootInit", C_fix(flag)) ;
372
373    }
374
375    flag = IDASetStopTime(solver_handle->ida_mem, time_stop);
376    if (flag != IDA_SUCCESS) 
377         chicken_error("could not set stop time with IDASetStopTime", C_fix(flag)) ;   
378
379    return solver_handle;
380}
381
382
383void ida_reinit_solver (IDA_Solver_Handle* solver_handle, double t0, C_word y0, C_word yp0)
384{
385    int flag; N_Vector ytmp, yptmp;
386
387    ytmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yy), C_c_f64vector(y0));
388    yptmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yp), C_c_f64vector(yp0));
389
390    N_VScale(1.0, ytmp, solver_handle->yy);
391    N_VScale(1.0, yptmp, solver_handle->yp);
392
393    flag = IDAReInit(solver_handle->ida_mem, t0, solver_handle->yy, solver_handle->yp);
394    if (flag != CV_SUCCESS) 
395         chicken_error("could not set reinitialize solver time with IDAReInit", C_fix(flag)) 
396
397    N_VDestroy_Serial(ytmp);
398    N_VDestroy_Serial(yptmp);
399   
400    return;
401}
402
403
404void ida_destroy_solver (IDA_Solver_Handle* solver_handle)
405{
406    IDAFree(&(solver_handle->ida_mem));
407    N_VDestroy_Serial(solver_handle->yy);
408    N_VDestroy_Serial(solver_handle->yp);
409    N_VDestroy_Serial(solver_handle->id);
410    free(solver_handle);
411    return;
412}
413
414
415int ida_solve (IDA_Solver_Handle* solver_handle, double tout)
416{
417    int flag;
418
419    if (solver_handle->event_number > 0)
420    {
421        long int ngevals;
422        flag = IDAGetNumGEvals(solver_handle->ida_mem, &ngevals);
423
424        if (flag != IDA_SUCCESS) 
425          chicken_error("error in IDAGetNumGEvals", C_fix(flag)) ;     
426
427
428        if (ngevals == 0)
429        {
430         flag = IDASolve (solver_handle->ida_mem, 
431                          tout,
432                          &(solver_handle->tret),
433                          solver_handle->yy,
434                          solver_handle->yp, 
435                          IDA_ONE_STEP);
436
437            switch (flag)
438            {
439                case IDA_SUCCESS:      return 0;
440                case IDA_ROOT_RETURN:  return 0;
441                case IDA_TSTOP_RETURN: return 2;
442                default:               chicken_error("unknown status code returned by IDASolve", C_fix(flag)) ;
443
444            }
445        }
446    }
447
448    flag = IDASolve (solver_handle->ida_mem,
449                     tout,
450                     &(solver_handle->tret),
451                     solver_handle->yy,
452                     solver_handle->yp, 
453                     IDA_NORMAL);
454
455    switch (flag)
456    {
457     case IDA_SUCCESS:
458      return 0;
459     
460      case IDA_ROOT_RETURN:
461      flag = IDAGetRootInfo(solver_handle->ida_mem, solver_handle->events);
462      if (flag != IDA_SUCCESS) 
463         chicken_error("could not obtain even information wtih IDAGetRootInfo", C_fix(flag)) 
464      adjust_zero_crossings(solver_handle->yy, solver_handle->abstol);
465      adjust_zero_crossings(solver_handle->yp, solver_handle->abstol);
466      return 1;
467     
468      case IDA_TSTOP_RETURN:
469      return 2;
470     
471      default: 
472      chicken_error("unknown status code returned by IDASolve", C_fix(flag)) 
473
474      }
475    }
476
477
478
479CVODE_Solver_Handle *cvode_create_solver
480( int lmm,
481  int iter, 
482
483  double time_start,
484  double time_stop,
485
486  int variable_number,
487  C_word variables, 
488
489  int ewt,
490
491  int event_number, 
492  int* events, 
493
494  unsigned int data_index,
495
496  double abstol,
497  double reltol
498)
499{
500    CVODE_Solver_Handle* solver_handle;
501    assert ((solver_handle = malloc (sizeof(struct CVODE_Solver_Handle_struct))) != NULL);
502
503    solver_handle->tret = 0.0;
504    solver_handle->event_number = event_number;
505    solver_handle->events = events;
506
507    int flag = 0;
508    int i = 0;
509
510    solver_handle->syy = variables;
511    solver_handle->yy = N_VMake_Serial(variable_number, C_c_f64vector(variables));
512
513    solver_handle->abstol = abstol;
514    solver_handle->reltol = reltol;
515
516    solver_handle->cvode_mem = CVodeCreate(lmm, iter);
517    if (solver_handle->cvode_mem == NULL) 
518    {
519       chicken_error("could not allocate memory with CVodeCreate", C_SCHEME_UNDEFINED);
520    }
521
522    flag = CVodeSetUserData(solver_handle->cvode_mem, (void *)data_index);
523    if (flag != CV_SUCCESS) 
524       chicken_error("could not set user data with CVodeSetUserData", C_SCHEME_UNDEFINED);
525
526    flag = CVodeInit (solver_handle->cvode_mem, 
527                    (CVRhsFn)cvode_rhs_cb, 
528                    time_start, 
529                    solver_handle->yy);
530    if (flag != CV_SUCCESS) 
531        chicken_error("could not initialize solver with CVodeInit", C_fix(flag)) ;     
532
533
534    if (events > 0)
535    { 
536#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
537       flag = CVodeRootInit(solver_handle->cvode_mem, event_number, (CVRootFn)cvode_event_cb);
538#else
539       flag = CVodeRootInit(solver_handle->cvode_mem, event_number, (CVRootFn)cvode_event_cb, (void *)data_index);
540#endif
541       if (flag != CV_SUCCESS) 
542          chicken_error("could not initialize event variables with CVodeRootInit", C_fix(flag)) ;       
543
544    }
545
546    if (ewt > 0)
547    {
548       flag = CVodeWFtolerances(solver_handle->cvode_mem, (CVEwtFn)cvode_ewt_cb);
549       if (flag != CV_SUCCESS) 
550          chicken_error("could not set error weight function with CVodeWFtolerances", C_fix(flag)) ;   
551    } else
552    {
553       flag = CVodeSStolerances(solver_handle->cvode_mem, solver_handle->reltol, solver_handle->abstol);
554       if (flag != CV_SUCCESS) 
555          chicken_error("could not set tolerances with CVodeSStolerances", C_fix(flag)) ;       
556    }
557
558    flag = CVodeSetStopTime(solver_handle->cvode_mem, time_stop);
559    if (flag != CV_SUCCESS) 
560         chicken_error("could not set stop time with CVodeSetStopTime", C_fix(flag)) 
561
562    return solver_handle;
563}
564
565
566void cvode_reinit_solver (CVODE_Solver_Handle* solver_handle, double t0, C_word y0)
567{
568    int flag; N_Vector ytmp;
569
570    ytmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yy), C_c_f64vector(y0));
571
572    N_VScale(1.0, ytmp, solver_handle->yy);
573    flag = CVodeReInit(solver_handle->cvode_mem, t0, solver_handle->yy);
574    if (flag != CV_SUCCESS) 
575         chicken_error("could not set reinitialize solver time with CVodeReInit", C_fix(flag)) ;       
576
577    N_VDestroy_Serial(ytmp);
578   
579    return;
580}
581
582
583void cvode_destroy_solver (CVODE_Solver_Handle* solver_handle)
584{
585    CVodeFree(&(solver_handle->cvode_mem));
586    N_VDestroy_Serial(solver_handle->yy);
587    free(solver_handle);
588    return;
589}
590
591
592int cvode_solve (CVODE_Solver_Handle* solver_handle, double tout)
593{
594    int flag;
595
596    if (solver_handle->event_number > 0)
597    {
598        long int nevals;
599        flag = CVodeGetNumRhsEvals(solver_handle->cvode_mem, &nevals);
600        if (flag != CV_SUCCESS) 
601           chicken_error("error in CVodeGetNumRhsEvals", C_fix(flag)) ;
602
603        if (nevals == 0)
604        {
605         flag = CVode (solver_handle->cvode_mem, 
606                       tout,
607                       solver_handle->yy,
608                       &(solver_handle->tret),
609                       CV_ONE_STEP);
610
611            switch (flag)
612            {
613                case CV_SUCCESS:      return 0;
614                case CV_ROOT_RETURN:  return 1;
615                case CV_TSTOP_RETURN: return 2;
616                default:              chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
617
618            }
619        }
620    }
621
622    flag = CVode (solver_handle->cvode_mem, 
623                  tout,
624                  solver_handle->yy,
625                  &(solver_handle->tret),
626                  CV_NORMAL);
627
628    switch (flag)
629    {
630     case CV_SUCCESS:
631      return 0;
632     
633      case CV_ROOT_RETURN:
634      flag = CVodeGetRootInfo(solver_handle->cvode_mem, solver_handle->events);
635      if (flag != CV_SUCCESS) 
636         chicken_error("could not obtain even information wtih CVodeGetRootInfo", C_fix(flag)) ;       
637      adjust_zero_crossings(solver_handle->yy, solver_handle->abstol);
638      return 1;
639     
640      case CV_TSTOP_RETURN:
641      return 2;
642     
643      default: 
644      chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
645
646      }
647    }
648
649<#
650
651(define c-cvode-create-solver
652  (foreign-safe-lambda (nonnull-c-pointer CVODESolverHandle) 
653                  "cvode_create_solver" 
654                  int     ;; linear multistep method
655                  int     ;; iter
656                  double  ;; time_start
657                  double  ;; time_stop
658                 
659                  int ;; variable_number
660                  scheme-object ;; variables
661                 
662                  int       ;; ewt
663                 
664                  int ;; event_number
665                  s32vector ;; events
666                 
667                  unsigned-int ;; user data
668                 
669                  double ;; abstol
670                  double ;; reltol
671                  ))
672
673 
674(define c-ida-create-solver
675  (foreign-safe-lambda (nonnull-c-pointer IDASolverHandle) 
676                  "ida_create_solver" 
677                  double    ;; time_start
678                  double    ;; time_stop
679                 
680                  int ;; variable_number
681                  scheme-object ;; variables
682                  scheme-object ;; derivatives
683                  bool ;; calculate initial conditions
684                  s32vector ;; alg_or_diff
685                  bool
686                 
687                  int ;; event_number
688                  s32vector ;; events
689
690                  unsigned-int ;; user data
691
692                  double ;; abstol
693                  double ;; reltol
694                  ))
695
696
697(define (ida-create-solver
698         tstart tstop variables derivatives 
699         residual-main #!key 
700         (residual-init #f)
701         (residual-event #f)
702         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
703         (suppress #f)
704         (ic #f)
705         (user-data #f)
706         (events (make-s32vector 0))
707         (reltol  1.0e-6)
708         (abstol  1.0e-6)
709         )
710 
711  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
712 
713  (let ((n (f64vector-length variables))
714        (nev (s32vector-length events))
715        (data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
716
717    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
718
719    (let ((residual-main1
720           (if user-data
721               (lambda (t yy yp rr data) 
722                 (let ((yy1 (make-f64vector n))
723                       (yp1 (make-f64vector n)))
724                   (move-memory! yy yy1 (fx* 8 n))
725                   (move-memory! yp yp1 (fx* 8 n))
726                   (let ((v (residual-main t yy1 yp1 data)))
727                     (move-memory! v rr (fx* 8 n))
728                     )))
729               (lambda (t yy yp rr data) 
730                 (let ((yy1 (make-f64vector n))
731                       (yp1 (make-f64vector n)))
732                   (move-memory! yy yy1 (fx* 8 n))
733                   (move-memory! yp yp1 (fx* 8 n))
734                   (let ((v (residual-main t yy1 yp1)))
735                     (move-memory! v rr (* 8 n))
736                     )))
737               ))
738
739          (residual-init1
740           (and residual-init
741                (if user-data
742                    (lambda (t yy yp rr data) 
743                      (let ((yy1 (make-f64vector n))
744                            (yp1 (make-f64vector n)))
745                        (move-memory! yy yy1 (fx* 8 n))
746                        (move-memory! yp yp1 (fx* 8 n))
747                        (let ((v (residual-init t yy1 yp1 data)))
748                          (move-memory! v rr (fx* 8 n))
749                          )))
750                    (lambda (t yy yp rr data) 
751                      (let ((yy1 (make-f64vector n))
752                            (yp1 (make-f64vector n)))
753                        (move-memory! yy yy1 (fx* 8 n))
754                        (move-memory! yp yp1 (fx* 8 n))
755                        (let ((v (residual-init t yy1 yp1)))
756                          (move-memory! v rr (* 8 n))
757                          )))
758                    )))
759
760          (residual-event1 
761           (and residual-event
762                (if user-data
763                    (lambda (t yy yp rr data) 
764                      (let ((yy1 (make-f64vector n))
765                            (yp1 (make-f64vector n)))
766                        (move-memory! yy yy1 (fx* 8 n))
767                        (move-memory! yp yp1 (fx* 8 n))
768                        (let ((v (residual-event t yy1 yp1 data)))
769                          (move-memory! v rr (fx* 4 nev))
770                          )))
771                    (lambda (t yy yp rr data) 
772                      (let ((yy1 (make-f64vector n))
773                            (yp1 (make-f64vector n)))
774                        (move-memory! yy yy1 (fx* 8 n))
775                        (move-memory! yp yp1 (fx* 8 n))
776                        (let ((v (residual-event t yy1 yp1)))
777                          (move-memory! v rr (* 4 nev))
778                          )))
779                    )))
780
781          )
782
783      (ida-residual-main-global residual-main1)
784      (if residual-init (ida-residual-init-global residual-init1))
785      (if residual-event (ida-residual-event-global residual-event1)))
786   
787    (c-ida-create-solver
788     tstart tstop 
789     n (object-evict variables) (object-evict derivatives)
790     ic alg-or-diff suppress (s32vector-length events) events
791     data-index abstol reltol)
792
793    ))
794
795
796(define (ida-create-solver/unsafe
797         tstart tstop variables derivatives 
798         residual-main #!key 
799         (residual-init #f)
800         (residual-event #f)
801         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
802         (suppress #f)
803         (ic #f)
804         (user-data #f)
805         (events (make-s32vector 0))
806         (reltol  1.0e-6)
807         (abstol  1.0e-6)
808         )
809
810  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
811 
812  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
813
814    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
815
816    (ida-residual-main-global residual-main)
817    (if residual-init (ida-residual-init-global residual-init))
818    (if residual-event (ida-residual-event-global residual-event))
819
820   
821    (c-ida-create-solver 
822     tstart tstop 
823     (f64vector-length variables) (object-evict variables) (object-evict derivatives)
824     ic alg-or-diff suppress (s32vector-length events) events
825     data-index abstol reltol)
826
827    ))
828                       
829                           
830
831(define c-ida-reinit-solver
832  (foreign-safe-lambda void "ida_reinit_solver" (nonnull-c-pointer IDASolverHandle) double scheme-object scheme-object ))
833
834(define (ida-reinit-solver solver t0 y0 yp0)
835  (assert (f64vector? y0))
836  (c-ida-reinit-solver solver t0 y0 yp0))
837
838(define c-ida-destroy-solver
839  (foreign-safe-lambda void "ida_destroy_solver" (nonnull-c-pointer IDASolverHandle) ))
840
841(define (ida-destroy-solver solver)
842  (object-release (ida-syy solver))
843  (object-release (ida-syp solver))
844  (c-ida-destroy-solver solver))
845
846(define ida-solve
847  (foreign-safe-lambda int "ida_solve" (nonnull-c-pointer IDASolverHandle) double ))
848
849(define ida-get-last-order
850  (foreign-safe-lambda* int (((nonnull-c-pointer IDASolverHandle) handle))
851#<<EOF
852    int result;
853    IDAGetLastOrder (handle->ida_mem, &result);
854    C_return (result);
855EOF
856))
857
858(define ida-get-num-steps
859  (foreign-safe-lambda* long (((nonnull-c-pointer IDASolverHandle) handle))
860#<<EOF
861    long result;
862    IDAGetNumSteps (handle->ida_mem, &result);
863    C_return (result);
864EOF
865))
866
867(define ida-get-last-step
868  (foreign-safe-lambda* double (((nonnull-c-pointer IDASolverHandle) handle))
869#<<EOF
870    double result;
871    IDAGetLastStep (handle->ida_mem, &result);
872    C_return (result);
873EOF
874))
875
876(define ida-syy 
877  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
878#<<EOF
879    C_return (handle->syy);
880EOF
881))
882
883(define ida-syp
884  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
885#<<EOF
886    C_return (handle->syp);
887EOF
888))
889
890(define c_ida_yy 
891  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
892#<<EOF
893   memcpy(result,
894          NV_DATA_S(handle->yy),
895          (NV_LENGTH_S(handle->yy))*sizeof(double));
896EOF
897))
898
899(define c_ida_yy_length
900  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
901#<<EOF
902   C_return (NV_LENGTH_S(handle->yy));
903EOF
904))
905
906(define (ida-yy handle)
907  (let ((v (make-f64vector (c_ida_yy_length handle) 0.0)))
908    (c_ida_yy handle v)
909    v))
910
911(define c_ida_yp 
912  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
913#<<EOF
914   memcpy(result,
915          NV_DATA_S(handle->yp),
916          (NV_LENGTH_S(handle->yp))*sizeof(double));
917EOF
918))
919
920(define c_ida_yp_length
921  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
922#<<EOF
923   C_return (NV_LENGTH_S(handle->yp));
924EOF
925))
926
927(define (ida-yp handle)
928  (let ((v (make-f64vector (c_ida_yp_length handle) 0.0)))
929    (c_ida_yp handle v)
930    v))
931 
932 
933
934(define (cvode-create-solver
935         tstart tstop variables 
936         rhs-fn #!key 
937         (lmm cvode-lmm/adams)
938         (iter cvode-iter/functional)
939         (ewt-fn #f)
940         (event-fn #f)
941         (user-data #f)
942         (events (make-s32vector 0))
943         (reltol  1.0e-6)
944         (abstol  1.0e-6)
945         )
946
947  (let ((n (f64vector-length variables))
948        (nev (s32vector-length events))
949        (data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
950
951    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
952
953    (let ((rhs-fn1
954           (if user-data
955               (lambda (t yy yp data) 
956                 (let ((yy1 (make-f64vector n)))
957                   (move-memory! yy yy1 (fx* 8 n))
958                   (let ((v (rhs-fn t yy1 data)))
959                     (move-memory! v yp (fx* 8 n))
960                     )))
961               (lambda (t yy yp data) 
962                 (let ((yy1 (make-f64vector n)))
963                   (move-memory! yy yy1 (fx* 8 n))
964                   (let ((v (rhs-fn t yy1)))
965                     (move-memory! v yp (fx* 8 n))
966                     )))
967               ))
968          (ewt-fn1
969           (if user-data
970               (lambda (yy ewt data) 
971                 (let ((yy1 (make-f64vector n)))
972                   (move-memory! yy yy1 (fx* 8 n))
973                   (let ((v (ewt-fn yy1 data)))
974                     (move-memory! v ewt (fx* 8 n))
975                     )))
976               (lambda (yy ewt data) 
977                 (let ((yy1 (make-f64vector n)))
978                   (move-memory! yy yy1 (fx* 8 n))
979                   (let ((v (ewt-fn yy1)))
980                     (move-memory! v ewt (fx* 8 n))
981                     )))
982               ))
983          (event-fn1
984           (and event-fn
985                (if user-data
986                    (lambda (t yy gout data) 
987                      (let ((yy1 (make-f64vector n)))
988                        (move-memory! yy yy1 (fx* 8 n))
989                        (let ((v (event-fn t yy1 data)))
990                          (move-memory! v gout (fx* 4 nev))
991                          )))
992                    (lambda (t yy gout data) 
993                      (let ((yy1 (make-f64vector n)))
994                        (move-memory! yy yy1 (fx* 8 n))
995                        (let ((v (event-fn t yy1)))
996                          (move-memory! v gout (fx* 4 nev))
997                          )))
998               )))
999          )
1000   
1001      (cvode-rhs-global rhs-fn1)
1002      (if ewt-fn (cvode-ewt-global ewt-fn1))
1003      (if event-fn (cvode-event-global event-fn1)))
1004   
1005    (c-cvode-create-solver
1006     lmm iter tstart tstop 
1007     n (object-evict variables )
1008     (if ewt-fn 1 0) (s32vector-length events) events
1009     data-index abstol reltol)
1010    ))
1011
1012(define (cvode-create-solver/unsafe
1013         tstart tstop variables 
1014         rhs-fn #!key 
1015         (lmm cvode-lmm/adams)
1016         (iter cvode-iter/functional)
1017         (ewt-fn #f)
1018         (event-fn #f)
1019         (user-data #f)
1020         (events (make-s32vector 0))
1021         (reltol  1.0e-6)
1022         (abstol  1.0e-6)
1023         )
1024 
1025  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
1026
1027    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
1028
1029    (cvode-rhs-global rhs-fn)
1030    (if ewt-fn (cvode-ewt-global ewt-fn))
1031    (if event-fn (cvode-event-global event-fn))
1032   
1033    (c-cvode-create-solver
1034     lmm iter tstart tstop 
1035     (f64vector-length variables) (object-evict variables )
1036     (if ewt-fn 1 0) (s32vector-length events) events
1037     data-index abstol reltol)
1038    ))
1039
1040
1041(define c-cvode-reinit-solver
1042  (foreign-safe-lambda void "cvode_reinit_solver" (nonnull-c-pointer CVODESolverHandle) double scheme-object ))
1043
1044(define (cvode-reinit-solver solver t0 y0)
1045  (assert (f64vector? y0))
1046  (c-cvode-reinit-solver solver t0 y0))
1047
1048                       
1049(define c-cvode-destroy-solver
1050  (foreign-safe-lambda void "cvode_destroy_solver" (nonnull-c-pointer CVODESolverHandle) ))
1051
1052(define (cvode-destroy-solver solver)
1053  (object-release (cvode-syy solver))
1054  (c-cvode-destroy-solver solver))
1055
1056
1057
1058(define cvode-solve
1059  (foreign-safe-lambda int "cvode_solve" (nonnull-c-pointer CVODESolverHandle) double ))
1060                           
1061(define cvode-syy 
1062  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer CVODESolverHandle) handle))
1063#<<EOF
1064    C_return (handle->syy);
1065EOF
1066))
1067
1068(define c_cvode_yy 
1069  (foreign-safe-lambda* void (((nonnull-c-pointer CVODESolverHandle) handle) (f64vector result))
1070#<<EOF
1071   memcpy(result,
1072          NV_DATA_S(handle->yy),
1073          (NV_LENGTH_S(handle->yy))*sizeof(double));
1074EOF
1075))
1076
1077(define c_cvode_yy_length
1078  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer CVODESolverHandle) handle))
1079#<<EOF
1080   C_return (NV_LENGTH_S(handle->yy));
1081EOF
1082))
1083
1084(define (cvode-yy handle)
1085  (let ((v (make-f64vector (c_cvode_yy_length handle) 0.0)))
1086    (c_cvode_yy handle v)
1087    v))
1088
1089(define cvode-get-last-order
1090  (foreign-safe-lambda* int (((nonnull-c-pointer CVODESolverHandle) handle))
1091#<<EOF
1092    int result;
1093    CVodeGetLastOrder (handle->cvode_mem, &result);
1094    C_return (result);
1095EOF
1096))
1097
1098(define cvode-get-num-steps
1099  (foreign-safe-lambda* long (((nonnull-c-pointer CVODESolverHandle) handle))
1100#<<EOF
1101    long result;
1102    CVodeGetNumSteps (handle->cvode_mem, &result);
1103    C_return (result);
1104EOF
1105))
1106
1107(define cvode-get-last-step
1108  (foreign-safe-lambda* double (((nonnull-c-pointer CVODESolverHandle) handle))
1109#<<EOF
1110    double result;
1111    CVodeGetLastStep (handle->cvode_mem, &result);
1112    C_return (result);
1113EOF
1114))
1115
1116
1117(define (pointer+-f64 p n)
1118  (pointer+ p (fx* 8 n)))
1119
1120
1121)
Note: See TracBrowser for help on using the repository browser.