Dev

.github/workflows/publish-to-pypi.yml

@@ -68,7 +68,7 @@ jobs:
         name: python-package-distributions
         path: dist/
     - name: Sign the dists with Sigstore
-      uses: sigstore/gh-action-sigstore-python@v2.1.1
+      uses: sigstore/gh-action-sigstore-python@v3.0.0
       with:
         inputs: >-
           ./dist/*.tar.gz

README.md

@@ -9,7 +9,8 @@ Solverz supports three types of abstract equation types, that are
 
 where $p$ is the parameter set of your models, $y_0$ is the previous time node value of $y$.
 
-For example, we want to know how long it takes for an apple to fall from a tree to the ground. We have the DAE 
+Say, we want to know how long it will take for an apple, launched from the ground into the air, to fall back to the
+ground. We have the differential equations 
 
 $$
 \begin{aligned}

Solverz/__init__.py

@@ -1,10 +1,9 @@
 from Solverz.equation.eqn import Eqn, Ode
 from Solverz.equation.equations import AE, FDAE, DAE
 from Solverz.equation.param import Param, IdxParam, TimeSeriesParam
-from Solverz.sym_algebra.symbols import idx, Para, iVar, iAliasVar
 from Solverz.sym_algebra.functions import (Sign, Abs, transpose, exp, Diag, Mat_Mul, sin, cos, Min, AntiWindUp,
                                            Saturation, heaviside, ln, MulVarFunc, UniVarFunc)
-from Solverz.variable.variables import Vars, TimeVars, as_Vars
+from Solverz.variable.variables import Vars, TimeVars
 from Solverz.solvers import *
 from Solverz.code_printer import made_numerical, module_printer
 from Solverz.utilities.io import save, load, save_result

Solverz/code_printer/python/module/test/test_module_generator.py

@@ -6,12 +6,13 @@
 import numpy as np
 from sympy import symbols, pycode, Integer
 
-from Solverz import as_Vars, Eqn, Ode, AE, sin, made_numerical, Model, Var, Param, TimeSeriesParam, AliasVar, Abs, exp, \
+from Solverz import Eqn, Ode, AE, sin, made_numerical, Model, Var, Param, TimeSeriesParam, AliasVar, Abs, exp, \
     cos
 from Solverz.code_printer.make_module import module_printer
 from Solverz.code_printer.python.inline.inline_printer import print_J_block
 from Solverz.code_printer.python.utilities import _print_var_parser
 from Solverz.sym_algebra.symbols import idx, iVar, Para
+from Solverz.variable.variables import as_Vars
 
 expected_dependency = r"""import os
 current_module_dir = os.path.dirname(os.path.abspath(__file__))

Solverz/equation/param.py

@@ -140,9 +140,9 @@ def get_v_t(self, t):
             # input of interp1d is zero-dimensional, we need to reshape
             # [0] is to eliminate the numpy DeprecationWarning in m3b9 test: Conversion of an array with ndim > 0 to a scalar is
             # deprecated, and will error in the future, which should be resolved.
-            vt = self.vt(t).reshape((-1,))[0]
+            vt = self.vt(t).reshape((-1,))
         else:
-            vt = self.v_series[-1]
+            vt = self.v_series[-1:]
 
         if self.index is not None:
             temp = self.v.copy()

Solverz/equation/test/test_DAE.py

@@ -1,6 +1,7 @@
 from sympy import Integer
-from Solverz import Ode, iVar, DAE, as_Vars, made_numerical
-
+from Solverz import Ode, DAE, made_numerical
+from Solverz.sym_algebra.symbols import iVar
+from Solverz.variable.variables import as_Vars
 
 def test_zero_index():
     # to test if index zero, being the sympy.Integer.instance, raises error when creating the singular mass mat

Solverz/equation/test/test_Param.py

@@ -46,7 +46,7 @@ def test_Param():
                          v_series=[0, 10, 100],
                          time_series=[0, 10, 20],
                          value=0)
-    assert Pb.get_v_t(5).__str__() == '5.0'
+    assert Pb.get_v_t(5).__str__() == '[5.]'
 
     try:
         Pb = TimeSeriesParam(name='Pb',

Solverz/model/basic.py

@@ -92,3 +92,13 @@ def create_instance(self):
             warnings.warn(f'Equation size {eqs.eqn_size} and variable size {eqs.vsize} not equal!')
 
         return eqs, y0
+
+    def add(self, m):
+        if isinstance(m, Model):
+            self.__dict__.update(m.__dict__)
+        elif isinstance(m, (Var, ParamBase, Eqn, Ode)):
+            self.__dict__[m.name] = m
+        elif isinstance(m, dict):
+            self.__dict__.update(m)
+        else:
+            raise ValueError(f"Unknown element type {type(m)}")

Solverz/solvers/daesolver/__init__.py

@@ -1,3 +1,4 @@
 from Solverz.solvers.daesolver.rodas import Rodas
+from Solverz.solvers.daesolver.ode15s import ode15s
 from Solverz.solvers.daesolver.beuler import backward_euler
 from Solverz.solvers.daesolver.trapezoidal import implicit_trapezoid

Solverz/solvers/daesolver/beuler.py

@@ -53,9 +53,12 @@ def backward_euler(dae: nDAE,
     t0 = tt
     Nt = int((T_end-T_initial)/dt) + 1
 
-    y = np.zeros((Nt, y0.shape[0]))
+    if opt.pbar:
+        pbar = tqdm(total=T_end - T_initial)
+
+    Y = np.zeros((Nt, y0.shape[0]))
     y0 = DaeIc(dae, y0, t0, opt.rtol)  # check and modify initial values
-    y[0, :] = y0
+    Y[0, :] = y0
     T = np.zeros((Nt,))
     T[0] = t0
 
@@ -73,13 +76,17 @@ def backward_euler(dae: nDAE,
 
         tt = tt + dt
         nt = nt + 1
-        y[nt] = y1
+        Y[nt] = y1
         T[nt] = tt
+        if opt.pbar:
+            pbar.update(T[nt] - T[nt - 1])
         y0 = y1
         t0 = tt
 
-    y = y[0:nt + 1]
+    Y = Y[0:nt + 1]
     T = T[0:nt + 1]
+    if opt.pbar:
+        pbar.close()
     stats.nstep = nt
 
-    return daesol(T, y, stats=stats)
+    return daesol(T, Y, stats=stats)

Solverz/solvers/daesolver/daeic.py

@@ -5,6 +5,21 @@
 
 from scipy.sparse.linalg import svds
 
+def getyp0(dae: nDAE, y0: np.ndarray, t0):
+
+    M = dae.M
+    p = dae.p
+    F0 = dae.F(t0, y0, p)
+    DiffEqn, DiffVar = M.nonzero()
+
+    yp0 = np.zeros_like(y0)
+
+    nonzero_rows = set(DiffEqn)
+    nonzero_cols = set(DiffVar)
+    Mp = M[list(nonzero_rows), :][:, list(nonzero_cols)]
+    yp0[DiffVar] = solve(Mp, F0[DiffEqn])
+    return yp0
+
 
 def DaeIc(dae: nDAE, y0: np.ndarray, t0, rtol):
     """

Solverz/solvers/daesolver/ode15s/__init__.py

@@ -0,0 +1 @@
+from .ode15s import ode15s

Solverz/solvers/daesolver/ode15s/ntrp15s.py

@@ -0,0 +1,13 @@
+from ..utilities import *
+
+
+def ntrp15s(tinterp, tnew, ynew, dt, dif, k):
+    s = (tinterp - tnew) / dt
+
+    if k == 1:
+        yinterp = ynew + dif[:, 0] * s
+    else:
+        kI = np.arange(1, k + 1).reshape((-1, 1))
+        yinterp = ynew + (dif[:, 0:k] @ np.cumprod((s + kI - 1) / kI, axis=0)).reshape(-1)
+
+    return yinterp