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

Last change on this file since 25538 was 25538, checked in by Ivan Raikov, 9 years ago

sundials: end time in cvode constructor is optional

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    if (time_stop > 0.0)
559    {
560     flag = CVodeSetStopTime(solver_handle->cvode_mem, time_stop);
561     if (flag != CV_SUCCESS) 
562         chicken_error("could not set stop time with CVodeSetStopTime", C_fix(flag)) 
563    }
564
565    return solver_handle;
566}
567
568
569void cvode_reinit_solver (CVODE_Solver_Handle* solver_handle, double t0, C_word y0)
570{
571    int flag; N_Vector ytmp;
572
573    ytmp = N_VMake_Serial(NV_LENGTH_S(solver_handle->yy), C_c_f64vector(y0));
574
575    N_VScale(1.0, ytmp, solver_handle->yy);
576    flag = CVodeReInit(solver_handle->cvode_mem, t0, solver_handle->yy);
577    if (flag != CV_SUCCESS) 
578         chicken_error("could not set reinitialize solver time with CVodeReInit", C_fix(flag)) ;       
579
580    N_VDestroy_Serial(ytmp);
581   
582    return;
583}
584
585
586void cvode_destroy_solver (CVODE_Solver_Handle* solver_handle)
587{
588    CVodeFree(&(solver_handle->cvode_mem));
589    N_VDestroy_Serial(solver_handle->yy);
590    free(solver_handle);
591    return;
592}
593
594
595int cvode_solve (CVODE_Solver_Handle* solver_handle, double tout)
596{
597    int flag;
598
599    if (solver_handle->event_number > 0)
600    {
601        long int nevals;
602        flag = CVodeGetNumRhsEvals(solver_handle->cvode_mem, &nevals);
603        if (flag != CV_SUCCESS) 
604           chicken_error("error in CVodeGetNumRhsEvals", C_fix(flag)) ;
605
606        if (nevals == 0)
607        {
608         flag = CVode (solver_handle->cvode_mem, 
609                       tout,
610                       solver_handle->yy,
611                       &(solver_handle->tret),
612                       CV_ONE_STEP);
613
614            switch (flag)
615            {
616                case CV_SUCCESS:      return 0;
617                case CV_ROOT_RETURN:  return 1;
618                case CV_TSTOP_RETURN: return 2;
619                default:              chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
620
621            }
622        }
623    }
624
625    flag = CVode (solver_handle->cvode_mem, 
626                  tout,
627                  solver_handle->yy,
628                  &(solver_handle->tret),
629                  CV_NORMAL);
630
631    switch (flag)
632    {
633     case CV_SUCCESS:
634      return 0;
635     
636      case CV_ROOT_RETURN:
637      flag = CVodeGetRootInfo(solver_handle->cvode_mem, solver_handle->events);
638      if (flag != CV_SUCCESS) 
639         chicken_error("could not obtain even information wtih CVodeGetRootInfo", C_fix(flag)) ;       
640      adjust_zero_crossings(solver_handle->yy, solver_handle->abstol);
641      return 1;
642     
643      case CV_TSTOP_RETURN:
644      return 2;
645     
646      default: 
647      chicken_error("unknown status code returned by CVode", C_fix(flag)) ;     
648
649      }
650    }
651
652<#
653
654(define c-cvode-create-solver
655  (foreign-safe-lambda (nonnull-c-pointer CVODESolverHandle) 
656                  "cvode_create_solver" 
657                  int     ;; linear multistep method
658                  int     ;; iter
659                  double  ;; time_start
660                  double  ;; time_stop
661                 
662                  int ;; variable_number
663                  scheme-object ;; variables
664                 
665                  int       ;; ewt
666                 
667                  int ;; event_number
668                  s32vector ;; events
669                 
670                  unsigned-int ;; user data
671                 
672                  double ;; abstol
673                  double ;; reltol
674                  ))
675
676 
677(define c-ida-create-solver
678  (foreign-safe-lambda (nonnull-c-pointer IDASolverHandle) 
679                  "ida_create_solver" 
680                  double    ;; time_start
681                  double    ;; time_stop
682                 
683                  int ;; variable_number
684                  scheme-object ;; variables
685                  scheme-object ;; derivatives
686                  bool ;; calculate initial conditions
687                  s32vector ;; alg_or_diff
688                  bool
689                 
690                  int ;; event_number
691                  s32vector ;; events
692
693                  unsigned-int ;; user data
694
695                  double ;; abstol
696                  double ;; reltol
697                  ))
698
699
700(define (ida-create-solver
701         tstart tstop variables derivatives 
702         residual-main #!key 
703         (residual-init #f)
704         (residual-event #f)
705         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
706         (suppress #f)
707         (ic #f)
708         (user-data #f)
709         (events (make-s32vector 0))
710         (reltol  1.0e-6)
711         (abstol  1.0e-6)
712         )
713 
714  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
715 
716  (let ((n (f64vector-length variables))
717        (nev (s32vector-length events))
718        (data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
719
720    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
721
722    (let ((residual-main1
723           (if user-data
724               (lambda (t yy yp rr data) 
725                 (let ((yy1 (make-f64vector n))
726                       (yp1 (make-f64vector n)))
727                   (move-memory! yy yy1 (fx* 8 n))
728                   (move-memory! yp yp1 (fx* 8 n))
729                   (let ((v (residual-main t yy1 yp1 data)))
730                     (move-memory! v rr (fx* 8 n))
731                     )))
732               (lambda (t yy yp rr data) 
733                 (let ((yy1 (make-f64vector n))
734                       (yp1 (make-f64vector n)))
735                   (move-memory! yy yy1 (fx* 8 n))
736                   (move-memory! yp yp1 (fx* 8 n))
737                   (let ((v (residual-main t yy1 yp1)))
738                     (move-memory! v rr (* 8 n))
739                     )))
740               ))
741
742          (residual-init1
743           (and residual-init
744                (if user-data
745                    (lambda (t yy yp rr data) 
746                      (let ((yy1 (make-f64vector n))
747                            (yp1 (make-f64vector n)))
748                        (move-memory! yy yy1 (fx* 8 n))
749                        (move-memory! yp yp1 (fx* 8 n))
750                        (let ((v (residual-init t yy1 yp1 data)))
751                          (move-memory! v rr (fx* 8 n))
752                          )))
753                    (lambda (t yy yp rr data) 
754                      (let ((yy1 (make-f64vector n))
755                            (yp1 (make-f64vector n)))
756                        (move-memory! yy yy1 (fx* 8 n))
757                        (move-memory! yp yp1 (fx* 8 n))
758                        (let ((v (residual-init t yy1 yp1)))
759                          (move-memory! v rr (* 8 n))
760                          )))
761                    )))
762
763          (residual-event1 
764           (and residual-event
765                (if user-data
766                    (lambda (t yy yp rr data) 
767                      (let ((yy1 (make-f64vector n))
768                            (yp1 (make-f64vector n)))
769                        (move-memory! yy yy1 (fx* 8 n))
770                        (move-memory! yp yp1 (fx* 8 n))
771                        (let ((v (residual-event t yy1 yp1 data)))
772                          (move-memory! v rr (fx* 4 nev))
773                          )))
774                    (lambda (t yy yp rr data) 
775                      (let ((yy1 (make-f64vector n))
776                            (yp1 (make-f64vector n)))
777                        (move-memory! yy yy1 (fx* 8 n))
778                        (move-memory! yp yp1 (fx* 8 n))
779                        (let ((v (residual-event t yy1 yp1)))
780                          (move-memory! v rr (* 4 nev))
781                          )))
782                    )))
783
784          )
785
786      (ida-residual-main-global residual-main1)
787      (if residual-init (ida-residual-init-global residual-init1))
788      (if residual-event (ida-residual-event-global residual-event1)))
789   
790    (c-ida-create-solver
791     tstart tstop 
792     n (object-evict variables) (object-evict derivatives)
793     ic alg-or-diff suppress (s32vector-length events) events
794     data-index abstol reltol)
795
796    ))
797
798
799(define (ida-create-solver/unsafe
800         tstart tstop variables derivatives 
801         residual-main #!key 
802         (residual-init #f)
803         (residual-event #f)
804         (alg-or-diff (make-s32vector (f64vector-length variables) 1))
805         (suppress #f)
806         (ic #f)
807         (user-data #f)
808         (events (make-s32vector 0))
809         (reltol  1.0e-6)
810         (abstol  1.0e-6)
811         )
812
813  (assert (= (f64vector-length variables) (f64vector-length derivatives)))
814 
815  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (ida-data-global)))) 0)))
816
817    (if user-data (ida-data-global (cons (cons data-index user-data) (ida-data-global))))
818
819    (ida-residual-main-global residual-main)
820    (if residual-init (ida-residual-init-global residual-init))
821    (if residual-event (ida-residual-event-global residual-event))
822
823   
824    (c-ida-create-solver 
825     tstart tstop 
826     (f64vector-length variables) (object-evict variables) (object-evict derivatives)
827     ic alg-or-diff suppress (s32vector-length events) events
828     data-index abstol reltol)
829
830    ))
831                       
832                           
833
834(define c-ida-reinit-solver
835  (foreign-safe-lambda void "ida_reinit_solver" (nonnull-c-pointer IDASolverHandle) double scheme-object scheme-object ))
836
837(define (ida-reinit-solver solver t0 y0 yp0)
838  (assert (f64vector? y0))
839  (c-ida-reinit-solver solver t0 y0 yp0))
840
841(define c-ida-destroy-solver
842  (foreign-safe-lambda void "ida_destroy_solver" (nonnull-c-pointer IDASolverHandle) ))
843
844(define (ida-destroy-solver solver)
845  (object-release (ida-syy solver))
846  (object-release (ida-syp solver))
847  (c-ida-destroy-solver solver))
848
849(define ida-solve
850  (foreign-safe-lambda int "ida_solve" (nonnull-c-pointer IDASolverHandle) double ))
851
852(define ida-get-last-order
853  (foreign-safe-lambda* int (((nonnull-c-pointer IDASolverHandle) handle))
854#<<EOF
855    int result;
856    IDAGetLastOrder (handle->ida_mem, &result);
857    C_return (result);
858EOF
859))
860
861(define ida-get-num-steps
862  (foreign-safe-lambda* long (((nonnull-c-pointer IDASolverHandle) handle))
863#<<EOF
864    long result;
865    IDAGetNumSteps (handle->ida_mem, &result);
866    C_return (result);
867EOF
868))
869
870(define ida-get-last-step
871  (foreign-safe-lambda* double (((nonnull-c-pointer IDASolverHandle) handle))
872#<<EOF
873    double result;
874    IDAGetLastStep (handle->ida_mem, &result);
875    C_return (result);
876EOF
877))
878
879(define ida-syy 
880  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
881#<<EOF
882    C_return (handle->syy);
883EOF
884))
885
886(define ida-syp
887  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
888#<<EOF
889    C_return (handle->syp);
890EOF
891))
892
893(define c_ida_yy 
894  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
895#<<EOF
896   memcpy(result,
897          NV_DATA_S(handle->yy),
898          (NV_LENGTH_S(handle->yy))*sizeof(double));
899EOF
900))
901
902(define c_ida_yy_length
903  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
904#<<EOF
905   C_return (NV_LENGTH_S(handle->yy));
906EOF
907))
908
909(define (ida-yy handle)
910  (let ((v (make-f64vector (c_ida_yy_length handle) 0.0)))
911    (c_ida_yy handle v)
912    v))
913
914(define c_ida_yp 
915  (foreign-safe-lambda* void (((nonnull-c-pointer IDASolverHandle) handle) (f64vector result))
916#<<EOF
917   memcpy(result,
918          NV_DATA_S(handle->yp),
919          (NV_LENGTH_S(handle->yp))*sizeof(double));
920EOF
921))
922
923(define c_ida_yp_length
924  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer IDASolverHandle) handle))
925#<<EOF
926   C_return (NV_LENGTH_S(handle->yp));
927EOF
928))
929
930(define (ida-yp handle)
931  (let ((v (make-f64vector (c_ida_yp_length handle) 0.0)))
932    (c_ida_yp handle v)
933    v))
934 
935 
936
937(define (cvode-create-solver
938         tstart tstop variables 
939         rhs-fn #!key 
940         (lmm cvode-lmm/adams)
941         (iter cvode-iter/functional)
942         (ewt-fn #f)
943         (event-fn #f)
944         (user-data #f)
945         (events (make-s32vector 0))
946         (reltol  1.0e-6)
947         (abstol  1.0e-6)
948         )
949
950  (let ((n (f64vector-length variables))
951        (nev (s32vector-length events))
952        (data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
953
954    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
955
956    (let ((rhs-fn1
957           (if user-data
958               (lambda (t yy yp data) 
959                 (let ((yy1 (make-f64vector n)))
960                   (move-memory! yy yy1 (fx* 8 n))
961                   (let ((v (rhs-fn t yy1 data)))
962                     (move-memory! v yp (fx* 8 n))
963                     )))
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)))
968                     (move-memory! v yp (fx* 8 n))
969                     )))
970               ))
971          (ewt-fn1
972           (if user-data
973               (lambda (yy ewt data) 
974                 (let ((yy1 (make-f64vector n)))
975                   (move-memory! yy yy1 (fx* 8 n))
976                   (let ((v (ewt-fn yy1 data)))
977                     (move-memory! v ewt (fx* 8 n))
978                     )))
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)))
983                     (move-memory! v ewt (fx* 8 n))
984                     )))
985               ))
986          (event-fn1
987           (and event-fn
988                (if user-data
989                    (lambda (t yy gout data) 
990                      (let ((yy1 (make-f64vector n)))
991                        (move-memory! yy yy1 (fx* 8 n))
992                        (let ((v (event-fn t yy1 data)))
993                          (move-memory! v gout (fx* 4 nev))
994                          )))
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)))
999                          (move-memory! v gout (fx* 4 nev))
1000                          )))
1001               )))
1002          )
1003   
1004      (cvode-rhs-global rhs-fn1)
1005      (if ewt-fn (cvode-ewt-global ewt-fn1))
1006      (if event-fn (cvode-event-global event-fn1)))
1007   
1008    (c-cvode-create-solver
1009     lmm iter tstart tstop 
1010     n (object-evict variables )
1011     (if ewt-fn 1 0) (s32vector-length events) events
1012     data-index abstol reltol)
1013    ))
1014
1015(define (cvode-create-solver/unsafe
1016         tstart tstop variables 
1017         rhs-fn #!key 
1018         (lmm cvode-lmm/adams)
1019         (iter cvode-iter/functional)
1020         (ewt-fn #f)
1021         (event-fn #f)
1022         (user-data #f)
1023         (events (make-s32vector 0))
1024         (reltol  1.0e-6)
1025         (abstol  1.0e-6)
1026         )
1027 
1028  (let ((data-index (if user-data (+ 1 (fold max 0 (map car (cvode-data-global)))) 0)))
1029
1030    (if user-data (cvode-data-global (cons (cons data-index user-data) (cvode-data-global))))
1031
1032    (cvode-rhs-global rhs-fn)
1033    (if ewt-fn (cvode-ewt-global ewt-fn))
1034    (if event-fn (cvode-event-global event-fn))
1035   
1036    (c-cvode-create-solver
1037     lmm iter tstart tstop 
1038     (f64vector-length variables) (object-evict variables )
1039     (if ewt-fn 1 0) (s32vector-length events) events
1040     data-index abstol reltol)
1041    ))
1042
1043
1044(define c-cvode-reinit-solver
1045  (foreign-safe-lambda void "cvode_reinit_solver" (nonnull-c-pointer CVODESolverHandle) double scheme-object ))
1046
1047(define (cvode-reinit-solver solver t0 y0)
1048  (assert (f64vector? y0))
1049  (c-cvode-reinit-solver solver t0 y0))
1050
1051                       
1052(define c-cvode-destroy-solver
1053  (foreign-safe-lambda void "cvode_destroy_solver" (nonnull-c-pointer CVODESolverHandle) ))
1054
1055(define (cvode-destroy-solver solver)
1056  (object-release (cvode-syy solver))
1057  (c-cvode-destroy-solver solver))
1058
1059
1060
1061(define cvode-solve
1062  (foreign-safe-lambda int "cvode_solve" (nonnull-c-pointer CVODESolverHandle) double ))
1063                           
1064(define cvode-syy 
1065  (foreign-safe-lambda* scheme-object (((nonnull-c-pointer CVODESolverHandle) handle))
1066#<<EOF
1067    C_return (handle->syy);
1068EOF
1069))
1070
1071(define c_cvode_yy 
1072  (foreign-safe-lambda* void (((nonnull-c-pointer CVODESolverHandle) handle) (f64vector result))
1073#<<EOF
1074   memcpy(result,
1075          NV_DATA_S(handle->yy),
1076          (NV_LENGTH_S(handle->yy))*sizeof(double));
1077EOF
1078))
1079
1080(define c_cvode_yy_length
1081  (foreign-safe-lambda* unsigned-int (((nonnull-c-pointer CVODESolverHandle) handle))
1082#<<EOF
1083   C_return (NV_LENGTH_S(handle->yy));
1084EOF
1085))
1086
1087(define (cvode-yy handle)
1088  (let ((v (make-f64vector (c_cvode_yy_length handle) 0.0)))
1089    (c_cvode_yy handle v)
1090    v))
1091
1092(define cvode-get-last-order
1093  (foreign-safe-lambda* int (((nonnull-c-pointer CVODESolverHandle) handle))
1094#<<EOF
1095    int result;
1096    CVodeGetLastOrder (handle->cvode_mem, &result);
1097    C_return (result);
1098EOF
1099))
1100
1101(define cvode-get-num-steps
1102  (foreign-safe-lambda* long (((nonnull-c-pointer CVODESolverHandle) handle))
1103#<<EOF
1104    long result;
1105    CVodeGetNumSteps (handle->cvode_mem, &result);
1106    C_return (result);
1107EOF
1108))
1109
1110(define cvode-get-last-step
1111  (foreign-safe-lambda* double (((nonnull-c-pointer CVODESolverHandle) handle))
1112#<<EOF
1113    double result;
1114    CVodeGetLastStep (handle->cvode_mem, &result);
1115    C_return (result);
1116EOF
1117))
1118
1119
1120(define (pointer+-f64 p n)
1121  (pointer+ p (fx* 8 n)))
1122
1123
1124)
Note: See TracBrowser for help on using the repository browser.