feat: sundials

exm/vdp.ml

@@ -31,7 +31,7 @@ let horizon _ = max_float
 let van_der_pol ()
 : (state, _, _, carray, carray, carray, zarray, carray) hnode
 = let yd         = cmake csize in
-  let zout       = cmake zsize in 
+  let zout       = cmake zsize in
   let fder _ _ y = fder y yd in
   let fzer _ _ _ = zout in
   let init _     = { lx=of_array [| x0; y0 |]; i=true } in

src/bin/main.ml

@@ -52,26 +52,28 @@ let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
                      else (module State.FunctionalSimState : State.SimState)
 
 let () =
-  (* if !sundials then *)
-  (*   if !greedy then *)
-  (*     (Format.eprintf "Sundials does not support greedy simulation\n"; exit 2) *)
-  (*   else *)
-  (*     let open StatefulSundials in *)
-  (*     let c = if !inplace then (module InPlace    : Csolver.Csolver) *)
-  (*                         else (module Functional : Csolver.Csolver) in *)
-  (*     let open (val c) in let open (val z) in *)
-  (*     let s = Solver.solver csolve (d_of_dc zsolve) in *)
-  (*     let open Sim.LazySim(val st) in run_until_n m s *)
-  (* else *)
   let open (val m) in
   let m = init !modelargs in
-  let open StatefulRK45 in
+  let sim =
-  let c = if !inplace then (module InPlace    : Csolver.CsolverC)
+    if !sundials then
-          else (module Functional : Csolver.CsolverC) in
+      if !greedy then
-  let open (val c) in let open (val z) in
+        (Format.eprintf "Sundials does not support greedy simulation\n"; exit 2)
-  let s = Solver.solver_c csolve zsolve in
+      else
-  let sim = if !greedy then let open Sim.GreedySim(val st) in run_until_n m s
+        let open StatefulSundials in
-            else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s) in
+        let c = if !inplace then (module InPlace    : Csolver.Csolver)
+                            else (module Functional : Csolver.Csolver) in
+        let open (val c) in let open (val z) in
+        let s = Solver.solver csolve (d_of_dc zsolve) in
+        let open Sim.LazySim(val st) in run_until_n m s
+    else
+      let open StatefulRK45 in
+      let c = if !inplace then (module InPlace    : Csolver.CsolverC)
+                          else (module Functional : Csolver.CsolverC) in
+      let open (val c) in let open (val z) in
+      let s = Solver.solver_c csolve zsolve in
+      if !greedy then let open Sim.GreedySim(val st) in run_until_n m s
+                 else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s)
+  in
   let open Solver in
   let HNode { init; csize; zsize; fder; fzer; cget; _ } = m in
   let state = init () in

src/lib/hsim/sim.ml

@@ -89,7 +89,7 @@ module LazySim (S : SimState) =
           match o with None -> state | Some o -> use o; loop state in
         loop state) (sim.init p) inputs
 
-    (** Run the model autonomously until [length], or until the model stops 
+    (** Run the model autonomously until [length], or until the model stops
         answering. *)
     let run_until model solver length =
       run_on model solver { start = 0.0; length; u = fun _ -> () }
@@ -135,7 +135,7 @@ module GreedySim (S : SimState) =
                 let ivp    = { fder; stop = stop -. now; init; size = model.csize } in
                 let zc     = { init; fzer; size = model.zsize } in
                 let ss     = solver.reset (ivp, zc) ss in
-                let i      = { start=i.start +. now; length=i.length -. now; 
+                let i      = { start=i.start +. now; length=i.length -. now;
                                u=Utils.offset i now } in
                 let mode, stop, now = Continuous, i.length, 0.0 in
                 let s = S.set_running ~mode ~input:i ~stop ~now s in
@@ -188,7 +188,7 @@ module GreedySim (S : SimState) =
         o::acc, state) ([], sim.init p) inputs in
       List.iter use (List.concat (List.rev o))
 
-    (** Run the model autonomously until [length], or until the model stops 
+    (** Run the model autonomously until [length], or until the model stops
         answering. *)
     let run_until model solver length =
       run_on model solver { start = 0.0; length; u = fun _ -> () }

src/lib/solvers/dune

@@ -2,7 +2,4 @@
 
 (library
   (name solvers)
-  (libraries
+  (libraries hsim sundialsml))
-    hsim
-    ;sundialsml
-    ))

src/lib/solvers/statefulRK45.ml

@@ -17,7 +17,7 @@ module Functional : Csolver.CsolverC =
         let state = initialize (fun _ _ _ -> ()) (vec v) in
         set_stop_time state 1.0; { state; vec=v } in
 
-      let reset { fder; init; stop; _ } _ = 
+      let reset { fder; init; stop; _ } _ =
         let fder t cvec dvec = Zls.blit (fder t cvec) dvec in
         let state = initialize fder (vec init) in
         set_stop_time state stop;

src/lib/solvers/statefulSundials.ml

@@ -1,4 +1,4 @@
-(*
+
 open Hsim.Types
 open Hsim.Solver
 open Zls
@@ -64,8 +64,8 @@ module InPlace : Csolver.Csolver =
         let y = Nvector_serial.wrap s.vec in
         let h, _ = solve_one_step s.state h y in
         let f t = get_dky s.state y t 0; Nvector_serial.unwrap y in
-        (h, f), s in 
+        (h, f), s in
 
       DNode { init; reset; step }
   end
-*)
+

src/lib/solvers/statefulZ.ml

@@ -64,6 +64,6 @@ module InPlace : Zsolver.ZsolverC =
         else (h, None), s in
 
       let copy _ = raise Common.Errors.TODO in
-      
+
       DNodeC { init; step; reset; copy }
   end

src/lib/solvers/zls.ml

@@ -42,7 +42,7 @@ let copy c1 = let c2 = cmake (length c1) in blit c1 c2; c2
 
 let blit_matrix m1 m2 = Array.iter2 blit m1 m2
 let copy_matrix m = Array.map copy m
- 
+
 type 's f_alloc = unit -> 's
 type 's f_maxsize = 's -> int * int
 type 's f_csize = 's -> int
@@ -126,7 +126,7 @@ module type STATE_ODE_SOLVER =
        time [t] since the last mesh-point or the initial instant. *)
     val get_dky : t -> dkyfn
 
-    
+
     (* generic solver parameters *)
     val set_stop_time  : t -> float -> unit
     val set_min_step   : t -> float -> unit
@@ -196,7 +196,7 @@ module type STATE_ZEROC_SOLVER =
        - the returned values is the simulation time of the earliest
          zero-crossing that was found. *)
     val find         : t -> ((float -> int -> unit) * carray) -> zarray -> float
-  
+
     (* locate the fields for which there is a takeoff *)
     val find_takeoff : t -> zarray -> float
   end