feat: start of lift, debugging, cleanup
This commit is contained in:
parent
883e5fff01
commit
589f89c768
GPG key ID:
7065F57ED8856128
31 changed files with 1297 additions and 51 deletions
119
src/bin/lift.ml
Normal file
119
src/bin/lift.ml
Normal file
|
|
@ -0,0 +1,119 @@
|
|||
|
||||
open Hsim.Types
|
||||
open Solvers.Zls
|
||||
open Common.Ztypes
|
||||
|
||||
type ('s, 'a) state =
|
||||
{ mutable state : 's; mutable input : 'a option; mutable time : time }
|
||||
|
||||
let lift f =
|
||||
let cstate =
|
||||
{ cvec = cmake 0; dvec = cmake 0; cindex = 0; zindex = 0;
|
||||
cend = 0; zend = 0; cmax = 0; zmax = 0;
|
||||
zinvec = zmake 0; zoutvec = cmake 0;
|
||||
major = false; horizon = max_float } in
|
||||
let Node { alloc=f_alloc; step=f_step; reset=f_reset } = f cstate in
|
||||
let state = { state = f_alloc (); input = None; time = 0.0 } in
|
||||
let csize, zsize = cstate.cmax, cstate.zmax in
|
||||
let no_roots_in = zmake zsize in
|
||||
let no_roots_out = cmake zsize in
|
||||
let ignore_der = cmake csize in
|
||||
let cstates = cmake csize in
|
||||
cstate.cvec <- cstates;
|
||||
f_reset state.state;
|
||||
|
||||
let no_time = -1.0 in
|
||||
|
||||
(* the function that compute the derivatives *)
|
||||
let fder { state; _ } input y =
|
||||
cstate.major <- false;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.cvec <- y;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
ignore (f_step state (no_time, input));
|
||||
cstate.dvec in
|
||||
|
||||
(* the function that compute the zero-crossings *)
|
||||
let fzer { state; _ } input y =
|
||||
Common.Debug.print "LIFT :: Calling [fzer]";
|
||||
cstate.major <- false;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.cvec <- y;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
ignore (f_step state (no_time, input));
|
||||
cstate.zoutvec in
|
||||
|
||||
(* the function which compute the output during integration *)
|
||||
let fout { state; _ } input y =
|
||||
cstate.major <- false;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.cvec <- y;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
f_step state (no_time, input) in
|
||||
|
||||
(* the function which compute a discrete step *)
|
||||
let step ({ state; time; _ } as st) input =
|
||||
cstate.major <- true;
|
||||
cstate.horizon <- infinity;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
let o = f_step state (time, input) in
|
||||
o, { st with state; input=Some input } in
|
||||
|
||||
let reset _ ({ state; _ } as st) = f_reset state; st in
|
||||
|
||||
(* horizon *)
|
||||
let horizon _ = cstate.horizon in
|
||||
|
||||
let jump _ = true in
|
||||
|
||||
(* the function which sets the zinvector into the *)
|
||||
(* internal zero-crossing variables *)
|
||||
let zset ({ state; input; _ } as st) zinvec =
|
||||
cstate.major <- false;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.zinvec <- zinvec;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
ignore (f_step state (no_time, Option.get input));
|
||||
st in
|
||||
|
||||
let cset ({ state; input; _ } as st) _ =
|
||||
cstate.major <- false;
|
||||
cstate.horizon <- infinity;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
ignore (f_step state (no_time, Option.get input));
|
||||
st in
|
||||
|
||||
let cget { state; input; _ } =
|
||||
cstate.major <- false;
|
||||
cstate.horizon <- infinity;
|
||||
cstate.zinvec <- no_roots_in;
|
||||
cstate.zoutvec <- no_roots_out;
|
||||
cstate.dvec <- ignore_der;
|
||||
cstate.cindex <- 0;
|
||||
cstate.zindex <- 0;
|
||||
ignore (f_step state (no_time, Option.get input));
|
||||
cstate.cvec in
|
||||
|
||||
HNode
|
||||
{ state; fder; fzer; step; fout; reset;
|
||||
horizon; cset; cget; zset; zsize; csize; jump }
|
||||
|
||||
|
|
@ -4,17 +4,25 @@ open Solvers
|
|||
open Examples
|
||||
open Common
|
||||
open Types
|
||||
open Lift
|
||||
|
||||
let sample = ref 10
|
||||
let stop = ref 30.0
|
||||
let sample = ref 1
|
||||
let stop = ref 10.0
|
||||
let accel = ref false
|
||||
let inplace = ref false
|
||||
let sundials = ref false
|
||||
let speed = ref false
|
||||
let steps = ref 1
|
||||
let model = ref None
|
||||
let minstep = ref None
|
||||
let maxstep = ref None
|
||||
let mintol = ref None
|
||||
let maxtol = ref None
|
||||
let no_print = ref false
|
||||
|
||||
let gt0i v i = v := if i <= 0 then 1 else i
|
||||
let gt0f v f = v := if f <= 0.0 then 1.0 else f
|
||||
let opt r s = r := Some s
|
||||
let modelargs = ref []
|
||||
|
||||
let set_model s =
|
||||
|
|
@ -30,25 +38,43 @@ let opts = [
|
|||
"-sundials", Arg.Set sundials, "\tUse sundials (doesn't support -accelerate)";
|
||||
"-inplace", Arg.Set inplace, "\tUse imperative solvers";
|
||||
"-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
|
||||
"-speed", Arg.Set speed, "\tLog the step length";
|
||||
"-minstep", Arg.String (opt minstep), "\tSet minimum solver step length";
|
||||
"-maxstep", Arg.String (opt maxstep), "\tSet maximum solver step length";
|
||||
"-mintol", Arg.String (opt mintol), "\tSet minimum solver tolerance";
|
||||
"-maxtol", Arg.String (opt maxtol), "\tSet maximum solver tolerance";
|
||||
"-no-print", Arg.Set no_print, "\tDo not print output values";
|
||||
]
|
||||
|
||||
let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
|
||||
|
||||
let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
|
||||
|
||||
let args = List.rev !modelargs
|
||||
let () = ignore lift
|
||||
|
||||
let m =
|
||||
try match !model with
|
||||
| None -> Format.eprintf "Missing model\n"; exit 2
|
||||
| Some "ball" -> Ball.init !modelargs
|
||||
| Some "vdp" -> Vdp.init !modelargs
|
||||
| Some "sincos" -> Sincos.init !modelargs
|
||||
| Some "sqrt" -> Sqrt.init !modelargs
|
||||
| Some "ball" -> Ball.init args
|
||||
| Some "vdp" -> Vdp.init args
|
||||
| Some "sincos" -> Sincos.init args
|
||||
| Some "sqrt" -> Sqrt.init args
|
||||
| Some "sin1x" -> Sin1x.init args
|
||||
| Some "sin1xd" -> Sin1x_der.init args
|
||||
| Some "ballz" -> lift Ballz.ball
|
||||
| Some "sincosz" -> lift Sincosz.f
|
||||
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
|
||||
with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
|
||||
|
||||
let st = if !inplace then (module State.InPlaceSimState : State.SimState)
|
||||
else (module State.FunctionalSimState : State.SimState)
|
||||
|
||||
let output =
|
||||
if !no_print then Output.ignore
|
||||
else if !speed then Output.print_h
|
||||
else Output.print (* Output.ignore *)
|
||||
|
||||
let sim =
|
||||
if !sundials then
|
||||
let open StatefulSundials in
|
||||
|
|
@ -57,7 +83,7 @@ let sim =
|
|||
let z = if !inplace then InPlace.zsolve else Functional.zsolve in
|
||||
let s = Solver.solver c (d_of_dc z) in
|
||||
let open Sim.Sim(val st) in
|
||||
run_until_n (Output.print !sample (run m s))
|
||||
run_until_n (output !sample (run m s))
|
||||
else
|
||||
let open StatefulRK45 in
|
||||
let c = if !inplace then InPlace.csolve else Functional.csolve in
|
||||
|
|
@ -66,7 +92,7 @@ let sim =
|
|||
let s = Solver.solver_c c z in
|
||||
let open Sim.Sim(val st) in
|
||||
let n = if !accel then accelerate m s else run m (d_of_dc s) in
|
||||
run_until_n (Output.print !sample n)
|
||||
run_until_n (output !sample n)
|
||||
|
||||
let () = sim !stop !steps ignore
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,6 @@
|
|||
|
||||
open Hsim.Types
|
||||
open Hsim.Utils
|
||||
open Common
|
||||
|
||||
let print_entry y t =
|
||||
let n = Bigarray.Array1.dim y in
|
||||
|
|
@ -13,6 +12,16 @@ let print_entry y t =
|
|||
Format.printf "\n";
|
||||
flush stdout
|
||||
|
||||
let print_entry_h y t h =
|
||||
let n = Bigarray.Array1.dim y in
|
||||
let rec loop i =
|
||||
if i = n then ()
|
||||
else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
|
||||
Format.printf "% .10e\t% .10e" t h;
|
||||
loop 0;
|
||||
Format.printf "\n";
|
||||
flush stdout
|
||||
|
||||
let print_sample samples ({ h; u; _ }, now) =
|
||||
let step = h /. (float_of_int samples) in
|
||||
let rec loop i =
|
||||
|
|
@ -20,8 +29,16 @@ let print_sample samples ({ h; u; _ }, now) =
|
|||
else if i = samples then print_entry (u h) (now +. h)
|
||||
else let t = float_of_int i *. step in
|
||||
(print_entry (u t) (now +. t); loop (i+1)) in
|
||||
if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end
|
||||
else begin Debug.print "C: "; loop 0 end
|
||||
if h <= 0.0 then print_entry (u 0.0) now else loop 0
|
||||
|
||||
let print_sample_h samples ({ h; u; _ }, now) =
|
||||
let step = h /. (float_of_int samples) in
|
||||
let rec loop i =
|
||||
if i > samples then ()
|
||||
else if i = samples then print_entry_h (u h) (now +. h) h
|
||||
else let t = float_of_int i *. step in
|
||||
(print_entry_h (u t) (now +. t) h; loop (i+1)) in
|
||||
if h <= 0.0 then print_entry_h (u 0.0) now h else loop 0
|
||||
|
||||
let print_limits { h; _ } =
|
||||
if h <= 0.0 then Format.printf "D: % .10e\n" 0.0
|
||||
|
|
@ -30,3 +47,18 @@ let print_limits { h; _ } =
|
|||
let print samples n =
|
||||
let DNode m = compose n (compose track (map (print_sample samples))) in
|
||||
DNode { m with reset=fun p -> m.reset (p, ((), ())) }
|
||||
|
||||
let print_h samples n =
|
||||
let DNode m = compose n (compose track (map (print_sample_h samples))) in
|
||||
DNode { m with reset=fun p -> m.reset (p, ((), ())) }
|
||||
|
||||
let ignore _ n =
|
||||
let state = () in
|
||||
let step () = function
|
||||
| None -> None, ()
|
||||
| Some _ -> Some (), () in
|
||||
let reset () () = () in
|
||||
let i = DNode { state; step; reset } in
|
||||
let DNode n = compose n i in
|
||||
DNode { n with reset=fun p -> n.reset (p, ()) }
|
||||
|
||||
|
|
|
|||
|
|
@ -2,3 +2,10 @@
|
|||
let debug = ref false
|
||||
|
||||
let print s = if !debug then Format.printf "%s\n" s else ()
|
||||
|
||||
let print_entry y =
|
||||
let n = Bigarray.Array1.dim y in
|
||||
let rec loop i =
|
||||
if i = n then ()
|
||||
else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
|
||||
if !debug then (loop 0; Format.printf "\n"; flush stdout)
|
||||
|
|
|
|||
151
src/lib/common/ztypes.ml
Normal file
151
src/lib/common/ztypes.ml
Normal file
|
|
@ -0,0 +1,151 @@
|
|||
(**************************************************************************)
|
||||
(* *)
|
||||
(* Zelus *)
|
||||
(* A synchronous language for hybrid systems *)
|
||||
(* http://zelus.di.ens.fr *)
|
||||
(* *)
|
||||
(* Marc Pouzet and Timothy Bourke *)
|
||||
(* *)
|
||||
(* Copyright 2012 - 2019. All rights reserved. *)
|
||||
(* *)
|
||||
(* This file is distributed under the terms of the CeCILL-C licence *)
|
||||
(* *)
|
||||
(* Zelus is developed in the INRIA PARKAS team. *)
|
||||
(* *)
|
||||
(**************************************************************************)
|
||||
|
||||
(* Type declarations and values that must be linked with *)
|
||||
(* the generated code *)
|
||||
type 'a continuous = { mutable pos: 'a; mutable der: 'a }
|
||||
|
||||
type ('a, 'b) zerocrossing = { mutable zin: 'a; mutable zout: 'b }
|
||||
|
||||
type 'a signal = 'a * bool
|
||||
type zero = bool
|
||||
|
||||
(* a synchronous stream function with type 'a -D-> 'b *)
|
||||
(* is represented by an OCaml value of type ('a, 'b) node *)
|
||||
type ('a, 'b) node =
|
||||
Node:
|
||||
{ alloc : unit -> 's; (* allocate the state *)
|
||||
step : 's -> 'a -> 'b; (* compute a step *)
|
||||
reset : 's -> unit; (* reset/inialize the state *)
|
||||
} -> ('a, 'b) node
|
||||
|
||||
(* the same with a method copy *)
|
||||
type ('a, 'b) cnode =
|
||||
Cnode:
|
||||
{ alloc : unit -> 's; (* allocate the state *)
|
||||
copy : 's -> 's -> unit; (* copy the source into the destination *)
|
||||
step : 's -> 'a -> 'b; (* compute a step *)
|
||||
reset : 's -> unit; (* reset/inialize the state *)
|
||||
} -> ('a, 'b) cnode
|
||||
|
||||
open Bigarray
|
||||
|
||||
type time = float
|
||||
type cvec = (float, float64_elt, c_layout) Array1.t
|
||||
type dvec = (float, float64_elt, c_layout) Array1.t
|
||||
type zinvec = (int32, int32_elt, c_layout) Array1.t
|
||||
type zoutvec = (float, float64_elt, c_layout) Array1.t
|
||||
|
||||
(* The interface with the ODE solver *)
|
||||
type cstate =
|
||||
{ mutable dvec : dvec; (* the vector of derivatives *)
|
||||
mutable cvec : cvec; (* the vector of positions *)
|
||||
mutable zinvec : zinvec; (* the vector of boolean; true when the
|
||||
solver has detected a zero-crossing *)
|
||||
mutable zoutvec : zoutvec; (* the corresponding vector that define
|
||||
zero-crossings *)
|
||||
mutable cindex : int; (* the position in the vector of positions *)
|
||||
mutable zindex : int; (* the position in the vector of zero-crossings *)
|
||||
mutable cend : int; (* the end of the vector of positions *)
|
||||
mutable zend : int; (* the end of the zero-crossing vector *)
|
||||
mutable cmax : int; (* the maximum size of the vector of positions *)
|
||||
mutable zmax : int; (* the maximum number of zero-crossings *)
|
||||
mutable horizon : float; (* the next horizon *)
|
||||
mutable major : bool; (* integration iff [major = false] *)
|
||||
}
|
||||
|
||||
(* A hybrid node is a node that is parameterised by a continuous state *)
|
||||
(* all instances points to this global parameter and read/write on it *)
|
||||
type ('a, 'b) hnode = cstate -> (time * 'a, 'b) node
|
||||
|
||||
type 'b hsimu =
|
||||
Hsim:
|
||||
{ alloc : unit -> 's;
|
||||
(* allocate the initial state *)
|
||||
maxsize : 's -> int * int;
|
||||
(* returns the max length of the *)
|
||||
(* cvector and zvector *)
|
||||
csize : 's -> int;
|
||||
(* returns the current length of the continuous state vector *)
|
||||
zsize : 's -> int;
|
||||
(* returns the current length of the zero-crossing vector *)
|
||||
step : 's -> cvec -> dvec -> zinvec -> time -> 'b;
|
||||
(* computes a step *)
|
||||
derivative : 's -> cvec -> dvec -> zinvec -> zoutvec -> time -> unit;
|
||||
(* computes the derivative *)
|
||||
crossings : 's -> cvec -> zinvec -> zoutvec -> time -> unit;
|
||||
(* computes the zero-crossings *)
|
||||
reset : 's -> unit;
|
||||
(* resets the state *)
|
||||
horizon : 's -> time;
|
||||
(* gives the next time horizon *)
|
||||
} -> 'b hsimu
|
||||
|
||||
(* a function with type 'a -C-> 'b, when given to a solver, is *)
|
||||
(* represented by an OCaml value of type ('a, 'b) hsnode *)
|
||||
type ('a, 'b) hsnode =
|
||||
Hnode:
|
||||
{ state : 's;
|
||||
(* the discrete state *)
|
||||
zsize : int;
|
||||
(* the maximum size of the zero-crossing vector *)
|
||||
csize : int;
|
||||
(* the maximum size of the continuous state vector (positions) *)
|
||||
derivative : 's -> 'a -> time -> cvec -> dvec -> unit;
|
||||
(* computes the derivative *)
|
||||
crossing : 's -> 'a -> time -> cvec -> zoutvec -> unit;
|
||||
(* computes the derivative *)
|
||||
output : 's -> 'a -> cvec -> 'b;
|
||||
(* computes the zero-crossings *)
|
||||
setroots : 's -> 'a -> cvec -> zinvec -> unit;
|
||||
(* returns the zero-crossings *)
|
||||
majorstep : 's -> time -> cvec -> 'a -> 'b;
|
||||
(* computes a step *)
|
||||
reset : 's -> unit;
|
||||
(* resets the state *)
|
||||
horizon : 's -> time;
|
||||
(* gives the next time horizon *)
|
||||
} -> ('a, 'b) hsnode
|
||||
|
||||
(* An idea suggested by Adrien Guatto, 26/04/2021 *)
|
||||
(* provide a means to the type for input/outputs of nodes *)
|
||||
(* express them with GADT to ensure type safety *)
|
||||
(* type ('a, 'b) node =
|
||||
| Fun : { step : 'a -> 'b;
|
||||
typ_arg: 'a typ;
|
||||
typ_return: 'b typ
|
||||
} -> ('a, 'b) node
|
||||
| Node :
|
||||
{ state : 's; step : 's -> 'a -> 'b * 's;
|
||||
typ_arg: 'a typ;
|
||||
typ_state : 's typ;
|
||||
typ_return: 'b typ } -> ('a, 'b) node
|
||||
|
||||
and 'a typ =
|
||||
| Tunit : unit typ
|
||||
| Tarrow : 'a typ * 'b typ -> ('a * 'b) typ
|
||||
| Tint : int -> int typ
|
||||
| Ttuple : 'a typlist -> 'a typ
|
||||
| Tnode : 'a typ * 'b typ -> ('a,'b) node typ
|
||||
|
||||
and 'a typlist =
|
||||
| Tnil : unit typlist
|
||||
| Tpair : 'a typ * 'b typlist -> ('a * 'b) typlist
|
||||
|
||||
Q1: do it for records? sum types ? How?
|
||||
Q2: provide a "type_of" function for every introduced type?
|
||||
*)
|
||||
|
||||
|
|
@ -8,6 +8,7 @@ module Sim (S : SimState) =
|
|||
include S
|
||||
|
||||
let step_discrete s step hor fder fzer cget csize zsize jump reset =
|
||||
Common.Debug.print "SIMU :: DISCRETE :: start";
|
||||
let ms, ss = get_mstate s, get_sstate s in
|
||||
let i, now, stop = get_input s, get_now s, get_stop s in
|
||||
let o, ms = step ms (i.u now) in
|
||||
|
|
@ -26,9 +27,11 @@ module Sim (S : SimState) =
|
|||
let mode, stop, now = Continuous, i.h, 0.0 in
|
||||
update ms ss (set_running ~mode ~input:i ~stop ~now s)
|
||||
end else set_running ~mode:Continuous s in
|
||||
Common.Debug.print "SIMU :: DISCRETE :: end";
|
||||
Utils.dot o, s
|
||||
|
||||
let step_continuous s step cset fout zset =
|
||||
Common.Debug.print "SIMU :: CONTINUOUS :: start";
|
||||
let ms, ss = get_mstate s, get_sstate s in
|
||||
let i, now, stop = get_input s, get_now s, get_stop s in
|
||||
let (h, f, z), ss = step ss stop in
|
||||
|
|
@ -46,6 +49,7 @@ module Sim (S : SimState) =
|
|||
let s = set_running ~mode:Discrete ~now:h s in
|
||||
update (zset ms z) ss s, Discontinuous in
|
||||
let h = h -. now in
|
||||
Common.Debug.print "SIMU :: CONTINUOUS :: end";
|
||||
{ h; u=fout; c }, s, { h; c; u=fms }
|
||||
|
||||
(** Simulation of a model with any solver. *)
|
||||
|
|
@ -55,7 +59,7 @@ module Sim (S : SimState) =
|
|||
: ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
|
||||
= let state = get_init m.state s.state in
|
||||
let step_discrete st =
|
||||
let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget
|
||||
let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget
|
||||
m.csize m.zsize m.jump s.reset in
|
||||
Some o, s in
|
||||
let step_continuous st =
|
||||
|
|
@ -97,7 +101,7 @@ module Sim (S : SimState) =
|
|||
let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
|
||||
DNode { a with state }) al in
|
||||
Some o, (st, al) in
|
||||
|
||||
|
||||
let step_continuous (st, al) =
|
||||
let ({ h; _ } as o), st, u =
|
||||
step_continuous st s.step m.body.cset m.body.fout m.body.zset in
|
||||
|
|
|
|||
|
|
@ -69,7 +69,7 @@ type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
|
|||
(** Consider a node with state copying as a node without state copying. *)
|
||||
let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset }
|
||||
|
||||
(** Consider a model without assertions as a model with an empty list of
|
||||
(** Consider a model without assertions as a model with an empty list of
|
||||
assertions. *)
|
||||
let a_of_h (HNode body) = HNodeA { body; assertions=[] }
|
||||
|
||||
|
|
|
|||
|
|
@ -3,7 +3,9 @@
|
|||
(* part of the Zelus standard library. *)
|
||||
(* It is implemented with in-place modification of arrays. *)
|
||||
|
||||
let debug = ref false
|
||||
let debug () =
|
||||
(* false *)
|
||||
!Common.Debug.debug
|
||||
|
||||
let printf x = Format.printf x
|
||||
|
||||
|
|
@ -121,7 +123,7 @@ type t = {
|
|||
let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c =
|
||||
s.t1 <- t;
|
||||
g t1 c f1; (* fill f1, because it is immediately copied into f0 by next_mesh *)
|
||||
if !debug then (printf "z|---------- init(%.24e, ... ----------@." t;
|
||||
if debug () then (printf "z|---------- init(%.24e, ... ----------@." t;
|
||||
log_limit s.f1);
|
||||
s.bothf_valid <- false
|
||||
|
||||
|
|
@ -152,6 +154,7 @@ let num_roots { f0; _ } = Zls.length f0
|
|||
|
||||
(* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *)
|
||||
let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
|
||||
Common.Debug.print "ZSOL :: Calling [step]";
|
||||
(* swap f0 and f1; f0 takes the previous value of f1 *)
|
||||
s.f0 <- f1;
|
||||
s.t0 <- t1;
|
||||
|
|
@ -162,7 +165,7 @@ let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
|
|||
g t c s.f1;
|
||||
s.bothf_valid <- true;
|
||||
|
||||
if !debug then
|
||||
if debug () then
|
||||
(printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1;
|
||||
log_limits s.f0 s.f1)
|
||||
|
||||
|
|
@ -212,7 +215,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
|
|||
dky t_right 0; (* c = dky_0(t_right); update state *)
|
||||
ignore (update_roots calc_zc f_left (get_f_right f_right') roots);
|
||||
|
||||
if !debug then
|
||||
if debug () then
|
||||
(printf
|
||||
"z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@."
|
||||
t_left t_right ttol;
|
||||
|
|
@ -280,20 +283,20 @@ let find ({ g = g; bothf_valid = bothf_valid;
|
|||
|
||||
match check_interval calc_zc f_left f_mid with
|
||||
| SearchLeft ->
|
||||
if !debug then printf "z| (%.24e -- %.24e] %.24e@."
|
||||
if debug () then printf "z| (%.24e -- %.24e] %.24e@."
|
||||
t_left t_mid t_right;
|
||||
let alpha = if i >= 1 then alpha *. 0.5 else alpha in
|
||||
let n_mid = f_mid_from_f_right f_right' in
|
||||
seek (t_left, f_left, n_mid, t_mid, Some f_mid, alpha, i + 1)
|
||||
|
||||
| SearchRight ->
|
||||
if !debug then printf "z| %.24e (%.24e -- %.24e]@."
|
||||
if debug () then printf "z| %.24e (%.24e -- %.24e]@."
|
||||
t_left t_mid t_right;
|
||||
let alpha = if i >= 1 then alpha *. 2.0 else alpha in
|
||||
seek (t_mid, f_mid, f_left, t_right, f_right', alpha, i + 1)
|
||||
|
||||
| FoundMid ->
|
||||
if !debug then printf "z| %.24e [%.24e] %.24e@."
|
||||
if debug () then printf "z| %.24e [%.24e] %.24e@."
|
||||
t_left t_mid t_right;
|
||||
ignore (update_roots calc_zc f_left f_mid roots);
|
||||
let f_tmp = f_mid_from_f_right f_right' in
|
||||
|
|
@ -303,7 +306,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
|
|||
|
||||
if not bothf_valid then (clear_roots roots; assert false)
|
||||
else begin
|
||||
if !debug then
|
||||
if debug () then
|
||||
printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1;
|
||||
|
||||
match check_interval calc_zc f0 f1 with
|
||||
|
|
@ -314,7 +317,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
|
|||
end
|
||||
|
||||
| FoundMid -> begin
|
||||
if !debug then printf "z| zero-crossing at limit (%.24e)@." t1;
|
||||
if debug () then printf "z| zero-crossing at limit (%.24e)@." t1;
|
||||
ignore (update_roots calc_zc f0 f1 roots);
|
||||
s.bothf_valid <- false;
|
||||
t1
|
||||
|
|
|
|||
|
|
@ -51,7 +51,9 @@ module GenericODE (Butcher : BUTCHER_TABLEAU) : STATE_ODE_SOLVER =
|
|||
struct (* {{{1 *)
|
||||
open Bigarray
|
||||
|
||||
let debug = ref false (* !Debug.debug *)
|
||||
let debug () =
|
||||
false
|
||||
(* !Common.Debug.debug *)
|
||||
|
||||
let pow = 1.0 /. float(Butcher.order)
|
||||
|
||||
|
|
@ -274,7 +276,7 @@ struct (* {{{1 *)
|
|||
"odexx: step size < min step size (\n now=%.24e\n h=%.24e\n< min_step=%.24e)"
|
||||
t h s.min_step);
|
||||
|
||||
if !debug then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t;
|
||||
if debug () then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t;
|
||||
|
||||
let rec onestep (alreadyfailed: bool) h =
|
||||
|
||||
|
|
@ -288,11 +290,11 @@ struct (* {{{1 *)
|
|||
let tnew = if finished then max_t else t +. h *. (mA maxK) in
|
||||
mapinto ynew (make_newval y k maxK);
|
||||
f tnew ynew k.(maxK);
|
||||
if !debug then log_step t y k.(0) tnew ynew k.(maxK);
|
||||
if debug () then log_step t y k.(0) tnew ynew k.(maxK);
|
||||
|
||||
let err = h *. calculate_error (abs_tol /. rel_tol) k y ynew in
|
||||
if err > rel_tol then begin
|
||||
if !debug then Printf.printf "s| error exceeds tolerance\n";
|
||||
if debug () then Printf.printf "s| error exceeds tolerance\n";
|
||||
|
||||
if h <= hmin then failwith
|
||||
(Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e"
|
||||
|
|
|
|||
|
|
@ -22,6 +22,8 @@ module Functional =
|
|||
{ state; vec = init } in
|
||||
|
||||
let step ({ state ; vec=v } as s) h =
|
||||
Common.Debug.print "SOLVER STEP";
|
||||
Common.Debug.print_entry v;
|
||||
let y_nv = vec v in
|
||||
let h = step state h y_nv in
|
||||
let state = copy state in
|
||||
|
|
|
|||
|
|
@ -15,7 +15,10 @@ module Functional =
|
|||
vec = zmake 0 } in
|
||||
|
||||
let reset { fzer; init; size } { vec; _ } =
|
||||
let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in
|
||||
let fzer t cvec zout =
|
||||
let zout' = fzer t cvec in blit zout' zout in
|
||||
Common.Debug.print "ZSolver Reset";
|
||||
Common.Debug.print_entry init;
|
||||
{ state = initialize size fzer init;
|
||||
vec = if length vec = size then vec else zmake size } in
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue