Context Navigation

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

Last change on this file since 25539 was 25539, checked in by Ivan Raikov, 9 years ago
sundials: fixes to unit tests
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
81	static void chicken_panic (C_char *) C_noret;
82	static 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
90	static void chicken_throw_exception(C_word value) C_noret;
91	static 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
106	void 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
121	typedef 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
138	typedef 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));
184	EOF
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
250	C_externexport int ida_residual_init_cb(double,void ,void ,void *,unsigned int);
251	C_externexport int ida_residual_main_cb(double,void ,void ,void *,unsigned int);
252	C_externexport int ida_residual_event_cb(double,void ,void ,void *,unsigned int);
253
254	C_externexport int cvode_rhs_cb(double,void ,void ,unsigned int);
255	C_externexport int cvode_event_cb(double,void ,void ,unsigned int);
256	C_externexport int cvode_ewt_cb(void ,void ,unsigned int);
257
258
259	void 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
270	IDA_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
386	void 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
407	void 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
418	int 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
482	CVODE_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
572	void 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
589	void 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
598	int 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);
863	EOF
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);
872	EOF
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);
881	EOF
882	))
883
884	(define ida-syy
885	(foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
886	#<<EOF
887	C_return (handle->syy);
888	EOF
889	))
890
891	(define ida-syp
892	(foreign-safe-lambda* scheme-object (((nonnull-c-pointer IDASolverHandle) handle))
893	#<<EOF
894	C_return (handle->syp);
895	EOF
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));
904	EOF
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));
911	EOF
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));
925	EOF
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));
932	EOF
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);
1075	EOF
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));
1084	EOF
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));
1091	EOF
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);
1105	EOF
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);
1114	EOF
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);
1123	EOF
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.

Download in other formats: