source: project/release/4/sundials/tags/1.4/sundials.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: 30.2 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    if (time_stop > 0.0)
376    {
377      flag = IDASetStopTime(solver_handle->ida_mem, time_stop);
378    }
379    if (flag != IDA_SUCCESS) 
380         chicken_error("could not set stop time with IDASetStopTime", C_fix(flag)) ;   
381
382    return solver_handle;
383}
384
385
386void ida_reinit_solver (IDA_Solver_Handle* solver_handle, double t0, C_word y0, C_word yp0)
387{
388    int flag; N_Vector ytmp, yptmp;
389
390    ytmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yy), C_c_f64vector(y0));
391    yptmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yp), C_c_f64vector(yp0));
392
393    N_VScale(1.0, ytmp, solver_handle->yy);
394    N_VScale(1.0, yptmp, solver_handle->yp);
395
396    flag = IDAReInit(solver_handle->ida_mem, t0, solver_handle->yy, solver_handle->yp);
397    if (flag != CV_SUCCESS) 
398         chicken_error("could not set reinitialize solver time with IDAReInit", C_fix(flag)) 
399
400    N_VDestroy_Serial(ytmp);
401    N_VDestroy_Serial(yptmp);
402   
403    return;
404}
405
406
407void ida_destroy_solver (IDA_Solver_Handle* solver_handle)
408{
409    IDAFree(&(solver_handle->ida_mem));
410    N_VDestroy_Serial(solver_handle->yy);
411    N_VDestroy_Serial(solver_handle->yp);
412    N_VDestroy_Serial(solver_handle->id);
413    free(solver_handle);
414    return;
415}
416
417
418int ida_solve (IDA_Solver_Handle* solver_handle, double tout)
419{
420    int flag;
421
422    if (solver_handle->event_number > 0)
423    {
424        long int ngevals;
425        flag = IDAGetNumGEvals(solver_handle->ida_mem, &ngevals);
426
427        if (flag != IDA_SUCCESS) 
428          chicken_error("error in IDAGetNumGEvals", C_fix(flag)) ;     
429
430
431        if (ngevals == 0)
432        {
433         flag = IDASolve (solver_handle->ida_mem, 
434                          tout,
435                          &(solver_handle->tret),
436                          solver_handle->yy,
437                          solver_handle->yp, 
438                          IDA_ONE_STEP);
439
440            switch (flag)
441            {
442                case IDA_SUCCESS:      return 0;
443                case IDA_ROOT_RETURN:  return 0;
444                case IDA_TSTOP_RETURN: return 2;
445                default:               chicken_error("unknown status code returned by IDASolve", C_fix(flag)) ;
446
447            }
448        }
449    }
450
451    flag = IDASolve (solver_handle->ida_mem,
452                     tout,
453                     &(solver_handle->tret),
454                     solver_handle->yy,
455                     solver_handle->yp, 
456                     IDA_NORMAL);
457
458    switch (flag)
459    {
460     case IDA_SUCCESS:
461      return 0;
462     
463      case IDA_ROOT_RETURN:
464      flag = IDAGetRootInfo(solver_handle->ida_mem, solver_handle->events);
465      if (flag != IDA_SUCCESS) 
466         chicken_error("could not obtain even information wtih IDAGetRootInfo", C_fix(flag)) 
467      adjust_zero_crossings(solver_handle->yy, solver_handle->abstol);
468      adjust_zero_crossings(solver_handle->yp, solver_handle->abstol);
469      return 1;
470     
471      case IDA_TSTOP_RETURN:
472      return 2;
473     
474      default: 
475      chicken_error("unknown status code returned by IDASolve", C_fix(flag)) 
476
477      }
478    }
479
480
481
482CVODE_Solver_Handle *cvode_create_solver
483( int lmm,
484  int iter, 
485
486  double time_start,
487  double time_stop,
488
489  int variable_number,
490  C_word variables, 
491
492  int ewt,
493
494  int event_number, 
495  int* events, 
496
497  unsigned int data_index,
498
499  double abstol,
500  double reltol
501)
502{
503    CVODE_Solver_Handle* solver_handle;
504    assert ((solver_handle = malloc (sizeof(struct CVODE_Solver_Handle_struct))) != NULL);
505
506    solver_handle->tret = 0.0;
507    solver_handle->event_number = event_number;
508    solver_handle->events = events;
509
510    int flag = 0;
511    int i = 0;
512
513    solver_handle->syy = variables;
514    solver_handle->yy = N_VMake_Serial(variable_number, C_c_f64vector(variables));
515
516    solver_handle->abstol = abstol;
517    solver_handle->reltol = reltol;
518
519    solver_handle->cvode_mem = CVodeCreate(lmm, iter);
520    if (solver_handle->cvode_mem == NULL) 
521    {
522       chicken_error("could not allocate memory with CVodeCreate", C_SCHEME_UNDEFINED);
523    }
524
525    flag = CVodeSetUserData(solver_handle->cvode_mem, (void *)data_index);
526    if (flag != CV_SUCCESS) 
527       chicken_error("could not set user data with CVodeSetUserData", C_SCHEME_UNDEFINED);
528
529    flag = CVodeInit (solver_handle->cvode_mem, 
530                    (CVRhsFn)cvode_rhs_cb, 
531                    time_start, 
532                    solver_handle->yy);
533    if (flag != CV_SUCCESS) 
534        chicken_error("could not initialize solver with CVodeInit", C_fix(flag)) ;     
535
536
537    if (events > 0)
538    { 
539#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
540       flag = CVodeRootInit(solver_handle->cvode_mem, event_number, (CVRootFn)cvode_event_cb);
541#else
542       flag = CVodeRootInit(solver_handle->cvode_mem, event_number, (CVRootFn)cvode_event_cb, (void *)data_index);
543#endif
544       if (flag != CV_SUCCESS) 
545          chicken_error("could not initialize event variables with CVodeRootInit", C_fix(flag)) ;       
546
547    }
548
549    if (ewt > 0)
550    {
551       flag = CVodeWFtolerances(solver_handle->cvode_mem, (CVEwtFn)cvode_ewt_cb);
552       if (flag != CV_SUCCESS) 
553          chicken_error("could not set error weight function with CVodeWFtolerances", C_fix(flag)) ;   
554    } else
555    {
556       flag = CVodeSStolerances(solver_handle->cvode_mem, solver_handle->reltol, solver_handle->abstol);
557       if (flag != CV_SUCCESS) 
558          chicken_error("could not set tolerances with CVodeSStolerances", C_fix(flag)) ;       
559    }
560
561    if (time_stop > 0.0)
562    {
563     flag = CVodeSetStopTime(solver_handle->cvode_mem, time_stop);
564     if (flag != CV_SUCCESS) 
565         chicken_error("could not set stop time with CVodeSetStopTime", C_fix(flag)) 
566    }
567
568    return solver_handle;
569}
570
571
572void cvode_reinit_solver (CVODE_Solver_Handle* solver_handle, double t0, C_word y0)
573{
574    int flag; N_Vector ytmp;
575
576    ytmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yy), C_c_f64vector(y0));
577
578    N_VScale(1.0, ytmp, solver_handle->yy);
579    flag = CVodeReInit(solver_handle->cvode_mem, t0, solver_handle->yy);
580    if (flag != CV_SUCCESS) 
581         chicken_error("could not set reinitialize solver time with CVodeReInit", C_fix(flag)) ;       
582
583    N_VDestroy_Serial(ytmp);
584   
585    return;
586}
587
588
589void cvode_destroy_solver (CVODE_Solver_Handle* solver_handle)
590{
591    CVodeFree(&(solver_handle->cvode_mem));
592    N_VDestroy_Serial(solver_handle->yy);
593    free(solver_handle);
594    return;
595}
596
597
598int cvode_solve (CVODE_Solver_Handle* solver_handle, double tout)
599{
600    int flag;
601
602    if (solver_handle->event_number > 0)
603    {
604        long int nevals;
605        flag = CVodeGetNumRhsEvals(solver_handle->cvode_mem, &nevals);
606        if (flag != CV_SUCCESS) 
607           chicken_error("error in CVodeGetNumRhsEvals", C_fix(flag)) ;
608
609        if (nevals == 0)
610        {
611         flag = CVode (solver_handle->cvode_mem, 
612                       tout,
613                       solver_handle->yy,
614                       &(solver_handle->tret),
615                       CV_ONE_STEP);
616
617            switch (flag)
618            {
619                case CV_SUCCESS:      return 0;
620                case CV_ROOT_RETURN:  return 1;
621                case CV_TSTOP_RETURN: return 2;
622                default:              chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
623
624            }
625        }
626    }
627
628    flag = CVode (solver_handle->cvode_mem, 
629                  tout,
630                  solver_handle->yy,
631                  &(solver_handle->tret),
632                  CV_NORMAL);
633
634    switch (flag)
635    {
636     case CV_SUCCESS:
637      return 0;
638     
639      case CV_ROOT_RETURN:
640      flag = CVodeGetRootInfo(solver_handle->cvode_mem, solver_handle->events);
641      if (flag != CV_SUCCESS) 
642         chicken_error("could not obtain even information wtih CVodeGetRootInfo", C_fix(flag)) ;       
643      adjust_zero_crossings(solver_handle->yy, solver_handle->abstol);
644      return 1;
645     
646      case CV_TSTOP_RETURN:
647      return 2;
648     
649      default: 
650      chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
651
652      }
653    }
654
655<#
656
657(define c-cvode-create-solver
658  (foreign-safe-lambda (nonnull-c-pointer CVODESolverHandle) 
659                  "cvode_create_solver" 
660                  int     ;; linear multistep method
661                  int     ;; iter
662                  double  ;; time_start
663                  double  ;; time_stop
664                 
665                  int ;; variable_number
666                  scheme-object ;; variables
667                 
668                  int       ;; ewt
669                 
670                  int ;; event_number
671                  s32vector ;; events
672                 
673                  unsigned-int ;; user data
674                 
675                  double ;; abstol
676                  double ;; reltol
677                  ))
678
679 
680(define c-ida-create-solver
681  (foreign-safe-lambda (nonnull-c-pointer IDASolverHandle) 
682                  "ida_create_solver" 
683                  double    ;; time_start
684                  double    ;; time_stop
685                 
686                  int ;; variable_number
687                  scheme-object ;; variables
688                  scheme-object ;; derivatives
689                  bool ;; calculate initial conditions
690                  s32vector ;; alg_or_diff
691                  bool
692                 
693                  int ;; event_number
694                  s32vector ;; events
695
696                  unsigned-int ;; user data
697
698                  double ;; abstol
699                  double ;; reltol
700                  ))
701
702
703(define (ida-create-solver
704         tstart variables derivatives 
705         residual-main #!key 
706         (tstop 0.0)
707         (residual-init #f)
708         (residual-event #f)
709         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
710         (suppress #f)
711         (ic #f)
712         (user-data #f)
713         (events (make-s32vector 0))
714         (reltol  1.0e-6)
715         (abstol  1.0e-6)
716         )
717 
718  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
719 
720  (let ((n (f64vector-length variables))
721        (nev (s32vector-length events))
722        (data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
723
724    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
725
726    (let ((residual-main1
727           (if user-data
728               (lambda (t yy yp rr data) 
729                 (let ((yy1 (make-f64vector n))
730                       (yp1 (make-f64vector n)))
731                   (move-memory! yy yy1 (fx* 8 n))
732                   (move-memory! yp yp1 (fx* 8 n))
733                   (let ((v (residual-main t yy1 yp1 data)))
734                     (move-memory! v rr (fx* 8 n))
735                     )))
736               (lambda (t yy yp rr data) 
737                 (let ((yy1 (make-f64vector n))
738                       (yp1 (make-f64vector n)))
739                   (move-memory! yy yy1 (fx* 8 n))
740                   (move-memory! yp yp1 (fx* 8 n))
741                   (let ((v (residual-main t yy1 yp1)))
742                     (move-memory! v rr (* 8 n))
743                     )))
744               ))
745
746          (residual-init1
747           (and residual-init
748                (if user-data
749                    (lambda (t yy yp rr data) 
750                      (let ((yy1 (make-f64vector n))
751                            (yp1 (make-f64vector n)))
752                        (move-memory! yy yy1 (fx* 8 n))
753                        (move-memory! yp yp1 (fx* 8 n))
754                        (let ((v (residual-init t yy1 yp1 data)))
755                          (move-memory! v rr (fx* 8 n))
756                          )))
757                    (lambda (t yy yp rr data) 
758                      (let ((yy1 (make-f64vector n))
759                            (yp1 (make-f64vector n)))
760                        (move-memory! yy yy1 (fx* 8 n))
761                        (move-memory! yp yp1 (fx* 8 n))
762                        (let ((v (residual-init t yy1 yp1)))
763                          (move-memory! v rr (* 8 n))
764                          )))
765                    )))
766
767          (residual-event1 
768           (and residual-event
769                (if user-data
770                    (lambda (t yy yp rr data) 
771                      (let ((yy1 (make-f64vector n))
772                            (yp1 (make-f64vector n)))
773                        (move-memory! yy yy1 (fx* 8 n))
774                        (move-memory! yp yp1 (fx* 8 n))
775                        (let ((v (residual-event t yy1 yp1 data)))
776                          (move-memory! v rr (fx* 4 nev))
777                          )))
778                    (lambda (t yy yp rr data) 
779                      (let ((yy1 (make-f64vector n))
780                            (yp1 (make-f64vector n)))
781                        (move-memory! yy yy1 (fx* 8 n))
782                        (move-memory! yp yp1 (fx* 8 n))
783                        (let ((v (residual-event t yy1 yp1)))
784                          (move-memory! v rr (* 4 nev))
785                          )))
786                    )))
787
788          )
789
790      (ida-residual-main-global residual-main1)
791      (if residual-init (ida-residual-init-global residual-init1))
792      (if residual-event (ida-residual-event-global residual-event1)))
793   
794    (c-ida-create-solver
795     tstart tstop 
796     n (object-evict variables) (object-evict derivatives)
797     ic alg-or-diff suppress (s32vector-length events) events
798     data-index abstol reltol)
799
800    ))
801
802
803(define (ida-create-solver/unsafe
804         tstart variables derivatives 
805         residual-main #!key 
806         (tstop 0.0)
807         (residual-init #f)
808         (residual-event #f)
809         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
810         (suppress #f)
811         (ic #f)
812         (user-data #f)
813         (events (make-s32vector 0))
814         (reltol  1.0e-6)
815         (abstol  1.0e-6)
816         )
817
818  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
819 
820  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
821
822    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
823
824    (ida-residual-main-global residual-main)
825    (if residual-init (ida-residual-init-global residual-init))
826    (if residual-event (ida-residual-event-global residual-event))
827
828   
829    (c-ida-create-solver 
830     tstart tstop 
831     (f64vector-length variables) (object-evict variables) (object-evict derivatives)
832     ic alg-or-diff suppress (s32vector-length events) events
833     data-index abstol reltol)
834
835    ))
836                       
837                           
838
839(define c-ida-reinit-solver
840  (foreign-safe-lambda void "ida_reinit_solver" (nonnull-c-pointer IDASolverHandle) double scheme-object scheme-object ))
841
842(define (ida-reinit-solver solver t0 y0 yp0)
843  (assert (f64vector? y0))
844  (c-ida-reinit-solver solver t0 y0 yp0))
845
846(define c-ida-destroy-solver
847  (foreign-safe-lambda void "ida_destroy_solver" (nonnull-c-pointer IDASolverHandle) ))
848
849(define (ida-destroy-solver solver)
850  (object-release (ida-syy solver))
851  (object-release (ida-syp solver))
852  (c-ida-destroy-solver solver))
853
854(define ida-solve
855  (foreign-safe-lambda int "ida_solve" (nonnull-c-pointer IDASolverHandle) double ))
856
857(define ida-get-last-order
858  (foreign-safe-lambda* int (((nonnull-c-pointer IDASolverHandle) handle))
859#<<EOF
860    int result;
861    IDAGetLastOrder (handle->ida_mem, &result);
862    C_return (result);
863EOF
864))
865
866(define ida-get-num-steps
867  (foreign-safe-lambda* long (((nonnull-c-pointer IDASolverHandle) handle))
868#<<EOF
869    long result;
870    IDAGetNumSteps (handle->ida_mem, &result);
871    C_return (result);
872EOF
873))
874
875(define ida-get-last-step
876  (foreign-safe-lambda* double (((nonnull-c-pointer IDASolverHandle) handle))
877#<<EOF
878    double result;
879    IDAGetLastStep (handle->ida_mem, &result);
880    C_return (result);
881EOF
882))
883
884(define ida-syy 
885  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
886#<<EOF
887    C_return (handle->syy);
888EOF
889))
890
891(define ida-syp
892  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
893#<<EOF
894    C_return (handle->syp);
895EOF
896))
897
898(define c_ida_yy 
899  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
900#<<EOF
901   memcpy(result,
902          NV_DATA_S(handle->yy),
903          (NV_LENGTH_S(handle->yy))*sizeof(double));
904EOF
905))
906
907(define c_ida_yy_length
908  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
909#<<EOF
910   C_return (NV_LENGTH_S(handle->yy));
911EOF
912))
913
914(define (ida-yy handle)
915  (let ((v (make-f64vector (c_ida_yy_length handle) 0.0)))
916    (c_ida_yy handle v)
917    v))
918
919(define c_ida_yp 
920  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
921#<<EOF
922   memcpy(result,
923          NV_DATA_S(handle->yp),
924          (NV_LENGTH_S(handle->yp))*sizeof(double));
925EOF
926))
927
928(define c_ida_yp_length
929  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
930#<<EOF
931   C_return (NV_LENGTH_S(handle->yp));
932EOF
933))
934
935(define (ida-yp handle)
936  (let ((v (make-f64vector (c_ida_yp_length handle) 0.0)))
937    (c_ida_yp handle v)
938    v))
939 
940 
941
942(define (cvode-create-solver
943         tstart variables 
944         rhs-fn #!key 
945         (tstop 0.0)
946         (lmm cvode-lmm/adams)
947         (iter cvode-iter/functional)
948         (ewt-fn #f)
949         (event-fn #f)
950         (user-data #f)
951         (events (make-s32vector 0))
952         (reltol  1.0e-6)
953         (abstol  1.0e-6)
954         )
955
956  (let ((n (f64vector-length variables))
957        (nev (s32vector-length events))
958        (data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
959
960    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
961
962    (let ((rhs-fn1
963           (if user-data
964               (lambda (t yy yp data) 
965                 (let ((yy1 (make-f64vector n)))
966                   (move-memory! yy yy1 (fx* 8 n))
967                   (let ((v (rhs-fn t yy1 data)))
968                     (move-memory! v yp (fx* 8 n))
969                     )))
970               (lambda (t yy yp data) 
971                 (let ((yy1 (make-f64vector n)))
972                   (move-memory! yy yy1 (fx* 8 n))
973                   (let ((v (rhs-fn t yy1)))
974                     (move-memory! v yp (fx* 8 n))
975                     )))
976               ))
977          (ewt-fn1
978           (if user-data
979               (lambda (yy ewt data) 
980                 (let ((yy1 (make-f64vector n)))
981                   (move-memory! yy yy1 (fx* 8 n))
982                   (let ((v (ewt-fn yy1 data)))
983                     (move-memory! v ewt (fx* 8 n))
984                     )))
985               (lambda (yy ewt data) 
986                 (let ((yy1 (make-f64vector n)))
987                   (move-memory! yy yy1 (fx* 8 n))
988                   (let ((v (ewt-fn yy1)))
989                     (move-memory! v ewt (fx* 8 n))
990                     )))
991               ))
992          (event-fn1
993           (and event-fn
994                (if user-data
995                    (lambda (t yy gout data) 
996                      (let ((yy1 (make-f64vector n)))
997                        (move-memory! yy yy1 (fx* 8 n))
998                        (let ((v (event-fn t yy1 data)))
999                          (move-memory! v gout (fx* 4 nev))
1000                          )))
1001                    (lambda (t yy gout data) 
1002                      (let ((yy1 (make-f64vector n)))
1003                        (move-memory! yy yy1 (fx* 8 n))
1004                        (let ((v (event-fn t yy1)))
1005                          (move-memory! v gout (fx* 4 nev))
1006                          )))
1007               )))
1008          )
1009   
1010      (cvode-rhs-global rhs-fn1)
1011      (if ewt-fn (cvode-ewt-global ewt-fn1))
1012      (if event-fn (cvode-event-global event-fn1)))
1013   
1014    (c-cvode-create-solver
1015     lmm iter tstart tstop 
1016     n (object-evict variables )
1017     (if ewt-fn 1 0) (s32vector-length events) events
1018     data-index abstol reltol)
1019    ))
1020
1021(define (cvode-create-solver/unsafe
1022         tstart variables 
1023         rhs-fn #!key 
1024         (tstop 0.0)
1025         (lmm cvode-lmm/adams)
1026         (iter cvode-iter/functional)
1027         (ewt-fn #f)
1028         (event-fn #f)
1029         (user-data #f)
1030         (events (make-s32vector 0))
1031         (reltol  1.0e-6)
1032         (abstol  1.0e-6)
1033         )
1034 
1035  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
1036
1037    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
1038
1039    (cvode-rhs-global rhs-fn)
1040    (if ewt-fn (cvode-ewt-global ewt-fn))
1041    (if event-fn (cvode-event-global event-fn))
1042   
1043    (c-cvode-create-solver
1044     lmm iter tstart tstop 
1045     (f64vector-length variables) (object-evict variables )
1046     (if ewt-fn 1 0) (s32vector-length events) events
1047     data-index abstol reltol)
1048    ))
1049
1050
1051(define c-cvode-reinit-solver
1052  (foreign-safe-lambda void "cvode_reinit_solver" (nonnull-c-pointer CVODESolverHandle) double scheme-object ))
1053
1054(define (cvode-reinit-solver solver t0 y0)
1055  (assert (f64vector? y0))
1056  (c-cvode-reinit-solver solver t0 y0))
1057
1058                       
1059(define c-cvode-destroy-solver
1060  (foreign-safe-lambda void "cvode_destroy_solver" (nonnull-c-pointer CVODESolverHandle) ))
1061
1062(define (cvode-destroy-solver solver)
1063  (object-release (cvode-syy solver))
1064  (c-cvode-destroy-solver solver))
1065
1066
1067
1068(define cvode-solve
1069  (foreign-safe-lambda int "cvode_solve" (nonnull-c-pointer CVODESolverHandle) double ))
1070                           
1071(define cvode-syy 
1072  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer CVODESolverHandle) handle))
1073#<<EOF
1074    C_return (handle->syy);
1075EOF
1076))
1077
1078(define c_cvode_yy 
1079  (foreign-safe-lambda* void (((nonnull-c-pointer CVODESolverHandle) handle) (f64vector result))
1080#<<EOF
1081   memcpy(result,
1082          NV_DATA_S(handle->yy),
1083          (NV_LENGTH_S(handle->yy))*sizeof(double));
1084EOF
1085))
1086
1087(define c_cvode_yy_length
1088  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer CVODESolverHandle) handle))
1089#<<EOF
1090   C_return (NV_LENGTH_S(handle->yy));
1091EOF
1092))
1093
1094(define (cvode-yy handle)
1095  (let ((v (make-f64vector (c_cvode_yy_length handle) 0.0)))
1096    (c_cvode_yy handle v)
1097    v))
1098
1099(define cvode-get-last-order
1100  (foreign-safe-lambda* int (((nonnull-c-pointer CVODESolverHandle) handle))
1101#<<EOF
1102    int result;
1103    CVodeGetLastOrder (handle->cvode_mem, &result);
1104    C_return (result);
1105EOF
1106))
1107
1108(define cvode-get-num-steps
1109  (foreign-safe-lambda* long (((nonnull-c-pointer CVODESolverHandle) handle))
1110#<<EOF
1111    long result;
1112    CVodeGetNumSteps (handle->cvode_mem, &result);
1113    C_return (result);
1114EOF
1115))
1116
1117(define cvode-get-last-step
1118  (foreign-safe-lambda* double (((nonnull-c-pointer CVODESolverHandle) handle))
1119#<<EOF
1120    double result;
1121    CVodeGetLastStep (handle->cvode_mem, &result);
1122    C_return (result);
1123EOF
1124))
1125
1126
1127(define (pointer+-f64 p n)
1128  (pointer+ p (fx* 8 n)))
1129
1130
1131)
Note: See TracBrowser for help on using the repository browser.