New features

docs/examples/ode/class_based.py

@@ -0,0 +1,113 @@
+# Authors: B. Malengier
+"""
+This example shows the most simple way of using a solver.
+We solve free vibration of a simple oscillator::
+        m \ddot{u} + k u = 0, u(0) = u_0, \dot{u}(0) = \dot{u}_0
+using the CVODE solver, which means we use a rhs function of \dot{u}.
+Solution::
+        u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m)
+
+"""
+from __future__ import print_function
+from numpy import asarray, cos, sin, sqrt
+import numpy as np
+from scikits.odes.sundials.cvode import CVODE, StatusEnum, CV_WrapJacRhsFunction
+from collections import namedtuple
+
+#data
+k = 4.0
+m = 1.0
+t1 = 10.
+#initial data on t=0, x[0] = u, x[1] = \dot{u}, xp = \dot{x}
+initx = [1, 0.1]
+
+def rhseqn(t, x):
+    """ we create rhs equations for the problem"""
+    return [
+        x[1],
+        - k/m * x[0]
+    ]
+
+def jaceqn(t, x,):
+    jac = np.zeros((2,2))
+    jac[0,1] = 1
+    jac[1,0] = -k/m
+    return jac
+
+
+def rootfn(t, x):
+    return (
+        x[0],
+        x[1],
+        t - t1,
+        np.sin(t),
+    )
+
+Root = namedtuple("Root", ["index", "rootsfound"])
+Results = namedtuple("Results", ["t", "x", "e"],)
+
+class System:
+    def __init__(self, dots, jac, events, num_events):
+        self._dots = dots
+        self._jac = jac
+        self._events = events
+        self.num_events = num_events
+        self.solver = CVODE(
+            self.dots,
+            jacfn=self.jac,
+            old_api=False,
+            one_step_compute=True,
+            rootfn=self.events,
+            nr_rootfns=self.num_events,
+        )
+
+    def dots(self, t, x, xdot,):# userdata=None,):
+        xdot[:] = self._dots(t, x)
+
+    def jac(self, t, x, xdot, jac, userdata=None,):
+        jac[...] = self._jac(t, x)
+
+    def events(self, t, x, g):
+        g[:] = self._events(t, x)
+
+    def simulate(self, t0, x0, tf, results=None):
+
+        if results is None:
+            results = Results([],[],[])
+
+        dense_t = results.t
+        dense_y = results.x
+        roots = results.e
+
+        solver = self.solver
+
+        solver.init_step(t0, x0)
+        solver.set_options(tstop=tf)
+
+        dense_t.append(np.copy(t0))
+        dense_y.append(np.copy(x0))
+
+        for cnt in range(1000):
+            res = solver.step(t1)
+            print(cnt, res.flag, res.values.t)
+            dense_t.append(np.copy(res.values.t))
+            dense_y.append(np.copy(res.values.y))
+            match res.flag:
+                case StatusEnum.ROOT_RETURN:
+                    rootsfound = solver.rootinfo()
+                    roots.append(Root(cnt, rootsfound))
+
+                    if res.values.t == tf:
+                        break
+
+                case StatusEnum.TSTOP_RETURN:
+                    #continue
+                    break
+
+        return results
+
+
+
+sys = System(rhseqn, jaceqn, rootfn, 4)
+res1 = sys.simulate(0., initx, 11.)
+res2 = sys.simulate(0., initx, 10.)

docs/examples/ode/cvodes_test.py

@@ -0,0 +1,81 @@
+# Authors: B. Malengier
+"""
+This example shows the most simple way of using a solver.
+We solve free vibration of a simple oscillator::
+        m \ddot{u} + k u = 0, u(0) = u_0, \dot{u}(0) = \dot{u}_0
+using the CVODE solver, which means we use a rhs function of \dot{u}.
+Solution::
+        u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m)
+
+"""
+from __future__ import print_function
+from numpy import asarray, cos, sin, sqrt
+import numpy as np
+
+#data
+userdata = dict(
+k = 4.0,
+m = 1.0,
+t1 = 10.,
+rhs_calls = 0,
+jac_calls = 0,
+)
+#initial data on t=0, x[0] = u, x[1] = \dot{u}, xp = \dot{x}
+initx = [1, 0.1]
+
+#define function for the right-hand-side equations which has specific signature
+def rhseqn(t, x, xdot, my_user_data):
+    """ we create rhs equations for the problem"""
+
+    k = my_user_data['k']
+    m = my_user_data['m']
+    my_user_data['rhs_calls'] += 1
+    xdot[0] = x[1]
+    xdot[1] = - k/m * x[0]
+
+def jaceqn(t, x, fx, jac, my_user_data):#=None):
+    my_user_data['jac_calls'] += 1
+    if my_user_data is None:
+        print("ERROR")
+        return 0
+
+    k = my_user_data['k']
+    m = my_user_data['m']
+
+    jac[0,1] = 1
+    jac[1,0] = -k/m
+
+#instantiate the solver
+if False:
+    from scikits.odes.sundials.cvode import CVODE, StatusEnum, CV_WrapJacRhsFunction
+    SolverClass = CVODE
+else:
+    from scikits.odes.sundials.cvode import CV_WrapJacRhsFunction
+    from scikits.odes.sundials.cvodes import CVODES, StatusEnum
+    SolverClass = CVODES
+
+from collections import namedtuple
+
+def rootfn(t, x, g, my_user_data):
+    t1 = my_user_data['t1']
+    g[0] = x[0]
+    g[1] = x[1]
+    g[2] = t - t1
+
+solver = SolverClass(
+    rhseqn, user_data=userdata,# jacfn=jaceqn,
+    old_api=False, one_step_compute=True,
+               rootfn=rootfn, nr_rootfns=3, )
+
+next_tstop = 10.
+#solver.init_step(0., initx)
+solver.set_options(tstop = next_tstop)
+dense_t = []
+dense_y = []
+roots = []
+Root = namedtuple("Root", ["index", "rootsfound"])
+print("starting loop")
+#print(solver.get_info(),)
+res = solver.solve([0, 10.0], initx)
+print("completed:",res)
+#print(solver.get_info(), solver.num_chk_pts, solver.options["rfn"])

docs/examples/ode/feature_test.py

@@ -0,0 +1,106 @@
+# Authors: B. Malengier
+"""
+This example shows the most simple way of using a solver.
+We solve free vibration of a simple oscillator::
+        m \ddot{u} + k u = 0, u(0) = u_0, \dot{u}(0) = \dot{u}_0
+using the CVODE solver, which means we use a rhs function of \dot{u}.
+Solution::
+        u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m)
+
+"""
+from __future__ import print_function
+from numpy import asarray, cos, sin, sqrt
+import numpy as np
+
+#data
+userdata = dict(
+k = 4.0,
+m = 1.0,
+t1 = 10.,
+rhs_calls = 0,
+jac_calls = 0,
+)
+#initial data on t=0, x[0] = u, x[1] = \dot{u}, xp = \dot{x}
+initx = [1, 0.1]
+
+#define function for the right-hand-side equations which has specific signature
+def rhseqn(t, x, xdot, my_user_data):
+    """ we create rhs equations for the problem"""
+
+    k = my_user_data['k']
+    m = my_user_data['m']
+    my_user_data['rhs_calls'] += 1
+    xdot[0] = x[1]
+    xdot[1] = - k/m * x[0]
+
+def jaceqn(t, x, fx, jac, my_user_data):#=None):
+    my_user_data['jac_calls'] += 1
+    if my_user_data is None:
+        print("ERROR")
+        return 0
+
+    k = my_user_data['k']
+    m = my_user_data['m']
+
+    jac[0,1] = 1
+    jac[1,0] = -k/m
+
+#instantiate the solver
+#from scikits.odes.sundials.cvode import CVODE, StatusEnum, CV_WrapJacRhsFunction
+#SolverClass = CVODE
+from scikits.odes.sundials.cvode import CV_WrapJacRhsFunction
+from scikits.odes.sundials.cvodes import CVODES, StatusEnum
+SolverClass = CVODES
+from collections import namedtuple
+
+def rootfn(t, x, g, my_user_data):
+    t1 = my_user_data['t1']
+    g[0] = x[0]
+    g[1] = x[1]
+    g[2] = t - t1
+
+solver = SolverClass(
+    rhseqn, user_data=userdata,# jacfn=jaceqn,
+    old_api=False, one_step_compute=True,
+               rootfn=rootfn, nr_rootfns=3, )
+
+next_tstop = 10.
+solver.init_step(0., initx)
+solver.set_options(tstop = next_tstop)
+dense_t = []
+dense_y = []
+roots = []
+Root = namedtuple("Root", ["index", "rootsfound"])
+print("starting loop")
+#print(solver.get_info(), solver.num_chk_pts, solver.options["rfn"])
+res = solver.solve([0, 10.0], initx)
+
+for cnt in range(1000):
+    #res = solver.step(1.)
+    s
+    print(cnt, res.flag, res.values.t)
+    dense_t.append(np.copy(res.values.t))
+    dense_y.append(np.copy(res.values.y))
+    match res.flag:
+        case StatusEnum.ROOT_RETURN:
+            rootsfound = solver.rootinfo()
+            roots.append(Root(cnt, rootsfound))
+
+        case StatusEnum.TSTOP_RETURN:
+            #continue
+            break
+
+
+    if res.values.t > 10.01:
+        print("broke from t")
+        break
+    #print(cnt, res)
+    #if res.values.y[0] <= 0.:
+    #    break
+t = np.array(dense_t)
+y = np.array(dense_y)
+print(solver.get_info())
+print(userdata)
+print(cnt)
+print(solver.num_chk_pts)
+

docs/examples/ode/simpleoscillator.py

@@ -25,7 +25,7 @@ def rhseqn(t, x, xdot):
 
 #instantiate the solver
 from scikits.odes import ode
-solver = ode('cvode', rhseqn, old_api=True)
+solver = ode('cvode', rhseqn, old_api=True, rtol=1E-12, atol=1E-15)
 #obtain solution at a required time
 result = solver.solve([0., 10., 20.], initx)
 
@@ -84,4 +84,4 @@ def scrhseqn(t, x):
     print('%4.2f %15.6g %15.6g' % (solver.t, solver.y[0], initx[0]*cos(sqrt(k/m)*solver.t)+initx[1]*sin(sqrt(k/m)*solver.t)/sqrt(k/m)))
 
 solver.integrate(solver.t+100)
-print('%4.2f %15.6g %15.6g' % (solver.t, solver.y[0], initx[0]*cos(sqrt(k/m)*solver.t)+initx[1]*sin(sqrt(k/m)*solver.t)/sqrt(k/m)))
+print('%4.2f %15.6g %15.6g' % (solver.t, solver.y[0], initx[0]*cos(sqrt(k/m)*solver.t)+initx[1]*sin(sqrt(k/m)*solver.t)/sqrt(k/m)))

docs/examples/ode/simpleoscillator_jac.py

@@ -10,6 +10,7 @@
 """
 from __future__ import print_function
 from numpy import asarray, cos, sin, sqrt
+import numpy as np
 
 #data
 k = 4.0
@@ -29,20 +30,63 @@ def jaceqn(t, x, fx, jac):
 
 #instantiate the solver
 from scikits.odes import ode
-solver = ode('cvode', rhseqn, jacfn=jaceqn)
+solver = ode('cvode', rhseqn, jacfn=jaceqn, )
 #obtain solution at a required time
-result = solver.solve([0., 1., 2.], initx)
+result = solver.solve([0., 10., 20.], initx)
 
 print('\n   t        Solution          Exact')
 print('------------------------------------')
-for t, u in zip(result[1], result[2]):
+#for t, u in zip(result[1], result[2]):
+for t, u in zip(result.values.t, result.values.y):
     print('%4.2f %15.6g %15.6g' % (t, u[0], initx[0]*cos(sqrt(k/m)*t)+initx[1]*sin(sqrt(k/m)*t)/sqrt(k/m)))
 
 #continue the solver
-result = solver.solve([result[1][-1], result[1][-1]+1], result[2][-1])
+#result = solver.solve([result[1][-1], result[1][-1]+1], result[2][-1])
+result = solver.solve([result.values.t[-1], result.values.t[-1]+1, result.values.t[-1]+110], result.values.y[-1])
 print('------------------------------------')
 print('  ...continuation of the solution')
 print('------------------------------------')
 
-for t, u in zip(result[1], result[2]):
+#for t, u in zip(result[1], result[2]):
+for t, u in zip(result.values.t, result.values.y):
     print ('%4.2f %15.6g %15.6g' % (t, u[0], initx[0]*cos(sqrt(k/m)*t)+initx[1]*sin(sqrt(k/m)*t)/sqrt(k/m)))
+
+
+from scikits.odes.sundials.cvode import CVODE, StatusEnum
+from collections import namedtuple
+
+def rootfn(t, x, g):
+    g[0] = x[0]
+    g[1] = x[1]
+    g[2] = t - 10.
+solver = CVODE(rhseqn, jacfn=jaceqn, old_api=False, one_step_compute=True,
+               rootfn=rootfn, nr_rootfns=3)
+
+next_tstop = 10.
+solver.init_step(0., initx)
+solver.set_options(tstop = next_tstop)
+dense_t = []
+dense_y = []
+roots = []
+Root = namedtuple("Root", ["index", "rootsfound"])
+for cnt in range(1000):
+    res = solver.step(1.)
+    print(cnt, res.flag, res.values.t)
+    dense_t.append(np.copy(res.values.t))
+    dense_y.append(np.copy(res.values.y))
+    match res.flag:
+        case StatusEnum.ROOT_RETURN:
+            rootsfound = solver.rootinfo()
+            roots.append(Root(cnt, rootsfound))
+
+        case StatusEnum.TSTOP_RETURN:
+            break
+t = np.array(dense_t)
+y = np.array(dense_y)
+
+
+    #if res.values.t > 9.99:
+    #    break
+    #print(cnt, res)
+    #if res.values.y[0] <= 0.:
+    #    break