Merge pull request #218 from cuihantao/develop

Develop

.github/workflows/pythonapp.yml

@@ -18,7 +18,6 @@ jobs:
         python -m pip install --upgrade pip
         python -m pip install -r requirements.txt
         python -m pip install -r requirements-dev.txt
-        python -m pip install -r requirements-opt.txt
         python -m pip install .
         python -m pip install nbmake==0.5 pytest-xdist line_profiler  # add'l packages for notebook tests.
     - name: Lint with flake8

andes/cli.py

@@ -139,6 +139,8 @@ def create_parser():
     prep_mode.add_argument('-f', '--full', action='store_true', help='full codegen')
     prep_mode.add_argument('-i', '--incremental', action='store_true',
                            help='rapid incrementally generate for updated models')
+    prep.add_argument('-c', '--compile', help='compile the code with numba after codegen',
+                      action='store_true', dest='precompile')
     prep.add_argument('--pycode-path', help='Save path for generated pycode')
     prep.add_argument('-m', '--models', nargs='*', help='model names to be individually prepared',
                       )

andes/core/model.py

@@ -58,6 +58,9 @@ def __getattr__(self, item):
 
         return self.__getattribute__(item)
 
+    def __getstate__(self):
+        return self.__dict__
+
     def add_callback(self, name: str, callback):
         """
         Add a cache attribute and a callback function for updating the attribute.
@@ -80,8 +83,6 @@ def refresh(self, name=None):
         name : str, list, optional
             name or list of cached to refresh, by default None for refreshing all
 
-        TODO: bug found in Example notebook 2. Time domain initialization fails
-        after refreshing.
         """
         if name is None:
             for name in self._callbacks.keys():
@@ -626,8 +627,6 @@ def __init__(self, system=None, config=None):
         self.tex_names = OrderedDict((('dae_t', 't_{dae}'),
                                       ('sys_f', 'f_{sys}'),
                                       ('sys_mva', 'S_{b,sys}'),
-                                      ('__ones', 'O_{nes}'),
-                                      ('__zeros', 'Z_{eros}'),
                                       ))
 
         # Model behavior flags
@@ -663,7 +662,7 @@ def __init__(self, system=None, config=None):
         self.cache.add_callback('e_setters', self._e_setters)
 
         self._input = OrderedDict()    # cached dictionary of inputs
-        self._input_z = OrderedDict()  # discrete flags in an OrderedDict
+        self._input_z = OrderedDict()  # discrete flags, storage only.
         self._rhs_f = OrderedDict()    # RHS of external f
         self._rhs_g = OrderedDict()    # RHS of external g
 
@@ -962,20 +961,43 @@ def refresh_inputs(self):
         for instance in self.cache.all_vars.values():
             self._input[instance.name] = instance.v
 
-        # append config variables
+        # append config variables as arrays
         for key, val in self.config.as_dict(refresh=True).items():
             self._input[key] = np.array(val)
 
-        # update`dae_t` and `sys_f`
+        # --- below are numpy scalars ---
+        # update`dae_t` and `sys_f`, and `sys_mva`
+        self._input['sys_f'] = np.array(self.system.config.freq, dtype=float)
+        self._input['sys_mva'] = np.array(self.system.config.mva, dtype=float)
         self._input['dae_t'] = self.system.dae.t
-        self._input['sys_f'] = self.system.config.freq
-        self._input['sys_mva'] = self.system.config.mva
 
-        # two vectors with the length of the number of devices.
-        #   Useful in the choices of `PieceWise`, which need to be vectors
-        #   for numba to compile.
-        self._input['__ones'] = np.ones(self.n)
-        self._input['__zeros'] = np.zeros(self.n)
+    def mock_refresh_inputs(self):
+        """
+        Use mock data to fill the inputs.
+
+        This function is used to generate input data of the desired type
+        to trigget JIT compilation.
+        """
+
+        self.get_inputs()
+        mock_arr = np.array([1.])
+
+        for key in self._input.keys():
+            try:
+                key_ndim = self._input[key].ndim
+            except AttributeError:
+                logger.error(key)
+                logger.error(self.class_name)
+
+            key_type = self._input[key].dtype
+
+            if key_ndim == 0:
+                self._input[key] = mock_arr.reshape(()).astype(key_type)
+            elif key_ndim == 1:
+                self._input[key] = mock_arr.astype(key_type)
+
+            else:
+                raise NotImplementedError("Unkonwn input data dimension %s" % key_ndim)
 
     def refresh_inputs_arg(self):
         """
@@ -1621,7 +1643,7 @@ def post_init_check(self):
             for item in self.services_icheck.values():
                 item.check()
 
-    def numba_jitify(self, parallel=False, cache=False, nopython=False):
+    def numba_jitify(self, parallel=False, cache=True, nopython=False):
         """
         Optionally convert `self.calls.f` and `self.calls.g` to
         JIT compiled functions.
@@ -1640,23 +1662,49 @@ def numba_jitify(self, parallel=False, cache=False, nopython=False):
         if self.flags.jited is True:
             return
 
-        self.calls.f = to_jit(self.calls.f,
-                              parallel=parallel,
-                              cache=cache,
-                              nopython=nopython,
-                              )
-        self.calls.g = to_jit(self.calls.g,
-                              parallel=parallel,
-                              cache=cache,
-                              nopython=nopython,
-                              )
+        kwargs = {'parallel': parallel,
+                  'cache': cache,
+                  'nopython': nopython,
+                  }
+
+        self.calls.f = to_jit(self.calls.f, **kwargs)
+        self.calls.g = to_jit(self.calls.g, **kwargs)
 
         for jname in self.calls.j:
-            self.calls.j[jname] = to_jit(self.calls.j[jname],
-                                         parallel=parallel, cache=cache)
+            self.calls.j[jname] = to_jit(self.calls.j[jname], **kwargs)
+
+        for name, instance in self.services_var.items():
+            if instance.v_str is not None:
+                self.calls.s[name] = to_jit(self.calls.s[name], **kwargs)
 
         self.flags.jited = True
 
+    def precompile(self):
+        """
+        Trigger numba compilation for this model.
+
+        This function requires the system to be setup, i.e.,
+        memory allocated for storage.
+        """
+
+        self.get_inputs()
+        if self.n == 0:
+            self.mock_refresh_inputs()
+        self.refresh_inputs_arg()
+
+        if callable(self.calls.f):
+            self.calls.f(*self.f_args)
+
+        if callable(self.calls.g):
+            self.calls.g(*self.g_args)
+
+        for jname, jfunc in self.calls.j.items():
+            jfunc(*self.j_args[jname])
+
+        for name, instance in self.services_var.items():
+            if instance.v_str is not None:
+                self.calls.s[name](*self.s_args[name])
+
     def __repr__(self):
         dev_text = 'device' if self.n == 1 else 'devices'
 

andes/core/symprocessor.py

@@ -138,8 +138,6 @@ def generate_symbols(self):
         self.inputs_dict['dae_t'] = Symbol('dae_t')
         self.inputs_dict['sys_f'] = Symbol('sys_f')
         self.inputs_dict['sys_mva'] = Symbol('sys_mva')
-        self.inputs_dict['__ones'] = Symbol('__ones')
-        self.inputs_dict['__zeros'] = Symbol('__zeros')
 
         # custom functions
         self.lambdify_func[0]['Indicator'] = lambda x: x
@@ -399,11 +397,11 @@ def generate_pycode(self, pycode_path, yapf_pycode):
 import numpy
 
 
-from numpy import ones_like, zeros_like, full, array                 # NOQA
+from numpy import ones_like, zeros_like, full, array                # NOQA
 from numpy import nan, pi, sin, cos, tan, sqrt, exp, select         # NOQA
 from numpy import greater_equal, less_equal, greater, less, equal   # NOQA
 from numpy import logical_and, logical_or, logical_not              # NOQA
-from numpy import real, imag, conj, angle, radians                  # NOQA
+from numpy import real, imag, conj, angle, radians, abs             # NOQA
 from numpy import arcsin, arccos, arctan, arctan2                   # NOQA
 from numpy import log                                               # NOQA
 

andes/main.py

@@ -685,7 +685,7 @@ def misc(edit_config='', save_config='', show_license=False, clean=True, recursi
 
 
 def prepare(quick=False, incremental=False, models=None,
-            nomp=False, **kwargs):
+            precompile=False, nomp=False, **kwargs):
     """
     Run code generation.
 
@@ -699,6 +699,8 @@ def prepare(quick=False, incremental=False, models=None,
     cli : bool
         True to indicate running from CLI.
         It will set `quick` to True if not `full`.
+    precompile : bool
+        True to compile model function calls after code generation.
 
     Warnings
     --------
@@ -728,6 +730,9 @@ def prepare(quick=False, incremental=False, models=None,
     system.prepare(quick=quick, incremental=incremental, models=models,
                    nomp=nomp, ncpu=ncpu)
 
+    if precompile:
+        system.precompile(models, nomp=nomp, ncpu=ncpu)
+
     if cli is True:
         return 0
     else:

andes/models/exciter/esst3a.py

@@ -196,11 +196,11 @@ def __init__(self, system, config):
 
         self.FEX = Piecewise(u=self.IN,
                              points=(0, 0.433, 0.75, 1),
-                             funs=('__ones + (__zeros * IN)',
-                                   '__ones * (1 - 0.577*IN)',
-                                   '__ones * sqrt(0.75 - IN ** 2)',
-                                   '__ones * 1.732*(1 - IN)',
-                                   '__zeros * IN'),
+                             funs=('1',
+                                   '1 - 0.577*IN',
+                                   'sqrt(0.75 - IN ** 2)',
+                                   '1.732 * (1 - IN)',
+                                   '0'),
                              info='Piecewise function FEX',
                              )
 

andes/models/renewable/reeca1.py

@@ -532,11 +532,11 @@ def __init__(self, system, config):
 
         self.VDL1 = Piecewise(u=self.s0_y,
                               points=('Vq1', 'Vq2', 'Vq3', 'Vq4'),
-                              funs=(f'({self.s0_y.name} * __zeros) + Iq1',
+                              funs=('Iq1',
                                     f'({self.s0_y.name} - Vq1) * kVq12 + Iq1',
                                     f'({self.s0_y.name} - Vq2) * kVq23 + Iq2',
                                     f'({self.s0_y.name} - Vq3) * kVq34 + Iq3',
-                                    f'({self.s0_y.name} * __zeros) + Iq4'),
+                                    'Iq4'),
                               tex_name='V_{DL1}',
                               info='Piecewise linear characteristics of Vq-Iq',
                               )
@@ -559,11 +559,11 @@ def __init__(self, system, config):
 
         self.VDL2 = Piecewise(u=self.s0_y,
                               points=('Vp1', 'Vp2', 'Vp3', 'Vp4'),
-                              funs=(f'({self.s0_y.name} * __zeros) + Ip1',
+                              funs=('Ip1',
                                     f'({self.s0_y.name} - Vp1) * kVp12 + Ip1',
                                     f'({self.s0_y.name} - Vp2) * kVp23 + Ip2',
                                     f'({self.s0_y.name} - Vp3) * kVp34 + Ip3',
-                                    f'({self.s0_y.name} * __zeros) + Ip4'),
+                                    'Ip4'),
                               tex_name='V_{DL2}',
                               info='Piecewise linear characteristics of Vp-Ip',
                               )
@@ -588,12 +588,12 @@ def __init__(self, system, config):
         Ipmax2sq = '(Imax**2 - IqHL_y**2)'
 
         # `Ipmax20`-squared (non-negative)
-        self.Ipmax2sq0 = ConstService(v_str=f'Piecewise((__zeros, Le({Ipmax2sq0}, 0.0)), ({Ipmax2sq0}, True), \
+        self.Ipmax2sq0 = ConstService(v_str=f'Piecewise((0, Le({Ipmax2sq0}, 0.0)), ({Ipmax2sq0}, True), \
                                               evaluate=False)',
                                       tex_name='I_{pmax20,nn}^2',
                                       )
 
-        self.Ipmax2sq = VarService(v_str=f'Piecewise((__zeros, Le({Ipmax2sq}, 0.0)), ({Ipmax2sq}, True), \
+        self.Ipmax2sq = VarService(v_str=f'Piecewise((0, Le({Ipmax2sq}, 0.0)), ({Ipmax2sq}, True), \
                                            evaluate=False)',
                                    tex_name='I_{pmax2}^2',
                                    )
@@ -615,12 +615,12 @@ def __init__(self, system, config):
 
         Iqmax2sq0 = '(Imax**2 - Ipcmd0**2)'  # initialization equation by using `Ipcmd0`
 
-        self.Iqmax2sq0 = ConstService(v_str=f'Piecewise((__zeros, Le({Iqmax2sq0}, 0.0)), ({Iqmax2sq0}, True), \
+        self.Iqmax2sq0 = ConstService(v_str=f'Piecewise((0, Le({Iqmax2sq0}, 0.0)), ({Iqmax2sq0}, True), \
                                               evaluate=False)',
                                       tex_name='I_{qmax,nn}^2',
                                       )
 
-        self.Iqmax2sq = VarService(v_str=f'Piecewise((__zeros, Le({Iqmax2sq}, 0.0)), ({Iqmax2sq}, True), \
+        self.Iqmax2sq = VarService(v_str=f'Piecewise((0, Le({Iqmax2sq}, 0.0)), ({Iqmax2sq}, True), \
                                            evaluate=False)',
                                    tex_name='I_{qmax2}^2')
 

andes/models/renewable/regca1.py

@@ -211,7 +211,7 @@ def __init__(self, system, config):
                                  )
 
         self.LVG = Piecewise(u=self.v, points=('Lvpnt0', 'Lvpnt1'),
-                             funs=('__zeros', '(v - Lvpnt0) * kLVG', '__ones'),
+                             funs=('0', '(v - Lvpnt0) * kLVG', '1'),
                              info='Ip gain during low voltage',
                              tex_name='L_{VG}',
                              )
@@ -229,7 +229,7 @@ def __init__(self, system, config):
                               points=('Zerox', 'Brkpt'),
                               funs=('0 + 9999*(1-Lvplsw)',
                                     '(S2_y - Zerox) * kLVPL + 9999 * (1-Lvplsw)',
-                                    '9999 * __ones'),
+                                    '9999'),
                               info='Low voltage Ipcmd upper limit',
                               tex_name='L_{VPL}',
                               )

andes/routines/pflow.py

@@ -8,7 +8,7 @@
 from andes.utils.misc import elapsed
 from andes.routines.base import BaseRoutine
 from andes.variables.report import Report
-from andes.shared import np, matrix, sparse, newton_krylov, IP_ADD
+from andes.shared import np, matrix, sparse, newton_krylov
 
 logger = logging.getLogger(__name__)
 
@@ -76,7 +76,7 @@ def init(self):
         self.system.init(self.models, routine='pflow')
         logger.info('Power flow initialized.')
 
-        # force precompile if numba is on - improves timing accuracy
+        # force compile if numba is on - improves timing accuracy
         if system.config.numba:
             system.f_update(self.models)
             system.g_update(self.models)
@@ -149,25 +149,17 @@ def summary(self):
         """
         Output a summary for the PFlow routine.
         """
-        ipadd_status = 'Standard (ipadd not available)'
 
         # extract package name, `kvxopt` or `kvxopt`
         sp_module = sparse.__module__
         if '.' in sp_module:
             sp_module = sp_module.split('.')[0]
 
-        if IP_ADD:
-            if self.system.config.ipadd:
-                ipadd_status = f'Fast in-place ({sp_module})'
-            else:
-                ipadd_status = 'Standard (ipadd disabled in config)'
-
         out = list()
         out.append('')
         out.append('-> Power flow calculation')
         out.append(f'{"Sparse solver":>16s}: {self.solver.sparselib.upper()}')
         out.append(f'{"Solution method":>16s}: {self.config.method} method')
-        out.append(f'{"Sparse addition":>16s}: {ipadd_status}')
         out_str = '\n'.join(out)
         logger.info(out_str)
 

andes/routines/tds.py

@@ -295,6 +295,7 @@ def run(self, no_pbar=False, no_summary=False, **kwargs):
             resume = True
             logger.debug("Resuming simulation from t=%.4fs.", system.dae.t)
             self._calc_h_first()
+            logger.debug("Initial step size for resumed simulation is h=%.4fs.", self.h)
 
         self.pbar = tqdm(total=100, ncols=70, unit='%', file=sys.stdout, disable=no_pbar)
 
@@ -406,6 +407,10 @@ def itm_step(self):
         system = self.system
         dae = self.system.dae
 
+        if self.h == 0:
+            logger.error("Current step size is zero. Integration is not permitted.")
+            return False
+
         self.mis = [1, 1]
         self.niter = 0
         self.converged = False