diff --git a/pipedream_solver/nsuperlink.py b/pipedream_solver/nsuperlink.py index f094ff8..6dae001 100644 --- a/pipedream_solver/nsuperlink.py +++ b/pipedream_solver/nsuperlink.py @@ -932,12 +932,13 @@ def orifice_flow_coefficients(self, u=None): _alpha_o = self._alpha_o # Orifice flow coefficient alpha_o _beta_o = self._beta_o # Orifice flow coefficient beta_o _chi_o = self._chi_o # Orifice flow coefficient chi_o + _unidir_o = self._unidir_o # If no input signal, assume orifice is closed if u is None: u = np.zeros(self.n_o, dtype=np.float64) # Specify orifice heads at previous timestep numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_inv_j, - _z_o, _tau_o, _Co, _Ao, _y_max_o, _J_uo, _J_do) + _z_o, _tau_o, _Co, _Ao, _y_max_o, _unidir_o, _J_uo, _J_do) # Export instance variables self._alpha_o = _alpha_o self._beta_o = _beta_o @@ -1424,13 +1425,14 @@ def solve_orifice_flows(self, dt, u=None): _Co = self._Co # Discharge coefficient of orifice o _Ao = self._Ao # Maximum flow area of orifice o _V_sj = self._V_sj + _unidir_o = self._unidir_o g = 9.81 # If no input signal, assume orifice is closed if u is None: u = np.zeros(self.n_o, dtype=np.float64) # Compute orifice flows _Qo_next = numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o, _tau_o, _y_max_o, _Co, _Ao, - _J_uo, _J_do, g) + _unidir_o, _J_uo, _J_do, g) # TODO: Move this inside numba function upstream_ctrl = (H_j[_J_uo] > H_j[_J_do]) _Qo_max = np.where(upstream_ctrl, _V_sj[_J_uo], _V_sj[_J_do]) / dt @@ -2519,10 +2521,12 @@ def xi_dk(dx_dk, B_dk, theta_dk, dt): return result @njit(int64(float64[:], float64[:], float64[:], float64[:], float64[:], float64[:], float64[:], - float64[:], float64[:], float64[:], float64[:], float64[:], int64[:], int64[:]), + float64[:], float64[:], float64[:], float64[:], float64[:], boolean[:], + int64[:], int64[:]), cache=True) def numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_inv_j, - _z_o, _tau_o, _Co, _Ao, _y_max_o, _J_uo, _J_do): + _z_o, _tau_o, _Co, _Ao, _y_max_o, _unidir_o, + _J_uo, _J_do): g = 9.81 _H_uo = H_j[_J_uo] _H_do = H_j[_J_do] @@ -2539,6 +2543,7 @@ def numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_i _z_o + _z_inv_uo + (_tau_o * _y_max_o * u / 2)) cond_2 = (_omega_o * _H_uo + (1 - _omega_o) * _H_do > _z_o + _z_inv_uo) + cond_3 = (_H_do >= _H_uo) & _unidir_o # Fill coefficient arrays # Submerged on both sides a = (cond_0 & cond_1) @@ -2559,17 +2564,18 @@ def numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_i _chi_o[c] = (_gamma_o[c] * (-1)**(1 - _omega_o[c]) * (- _z_inv_uo[c] - _z_o[c])) # No flow - d = (~cond_0 & ~cond_2) + d = (~cond_0 & ~cond_2) | cond_3 _alpha_o[d] = 0.0 _beta_o[d] = 0.0 _chi_o[d] = 0.0 return 1 @njit(float64[:](float64[:], float64[:], float64[:], float64[:], - float64[:], float64[:], float64[:], float64[:], int64[:], int64[:], float64), + float64[:], float64[:], float64[:], float64[:], boolean[:], + int64[:], int64[:], float64), cache=True) def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o, - _tau_o, _y_max_o, _Co, _Ao, _J_uo, _J_do, g=9.81): + _tau_o, _y_max_o, _Co, _Ao, _unidir_o, _J_uo, _J_do, g=9.81): # Specify orifice heads at previous timestep _H_uo = H_j[_J_uo] _H_do = H_j[_J_do] @@ -2590,6 +2596,7 @@ def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o, _z_o + _z_inv_uo + (_tau_o * _y_max_o * u / 2)) cond_2 = (_omega_o * _H_uo + (1 - _omega_o) * _H_do > _z_o + _z_inv_uo) + cond_3 = (_H_do >= _H_uo) & _unidir_o # Fill coefficient arrays # Submerged on both sides a = (cond_0 & cond_1) @@ -2610,7 +2617,7 @@ def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o, _chi_o[c] = (_gamma_o[c] * (-1)**(1 - _omega_o[c]) * (- _z_inv_uo[c] - _z_o[c])) # No flow - d = (~cond_0 & ~cond_2) + d = (~cond_0 & ~cond_2) | cond_3 _alpha_o[d] = 0.0 _beta_o[d] = 0.0 _chi_o[d] = 0.0 @@ -2716,8 +2723,8 @@ def numba_pump_flow_coefficients(_alpha_p, _beta_p, _chi_p, H_j, _z_inv_j, _Qp, cond_0 = _H_up > _z_inv_up + _z_p # Condition 1: Head difference is within range of pump curve cond_1 = (_dHp > _dHp_min) & (_dHp < _dHp_max) - _dHp[_dHp > _dHp_max] = _dHp_max - _dHp[_dHp < _dHp_min] = _dHp_min + _dHp[_dHp > _dHp_max] = _dHp_max[_dHp > _dHp_max] + _dHp[_dHp < _dHp_min] = _dHp_min[_dHp < _dHp_min] # Compute universal coefficients _gamma_p = gamma_p(_Qp, _b_p, _c_p, u) # Fill coefficient arrays @@ -2748,8 +2755,8 @@ def numba_solve_pump_flows(H_j, u, _z_inv_j, _z_p, _dHp_max, _dHp_min, _a_p, _b_ _z_inv_up = _z_inv_j[_J_up] # Create conditionals _dHp = _H_dp - _H_up - _dHp[_dHp > _dHp_max] = _dHp_max - _dHp[_dHp < _dHp_min] = _dHp_min + _dHp[_dHp > _dHp_max] = _dHp_max[_dHp > _dHp_max] + _dHp[_dHp < _dHp_min] = _dHp_min[_dHp < _dHp_min] cond_0 = _H_up > _z_inv_up + _z_p # Compute universal coefficients _Qp_next = (u / _b_p * (_a_p - _dHp))**(1 / _c_p) diff --git a/pipedream_solver/superlink.py b/pipedream_solver/superlink.py index 1348bdb..2e405a7 100644 --- a/pipedream_solver/superlink.py +++ b/pipedream_solver/superlink.py @@ -173,6 +173,7 @@ class SuperLink(): | A | float | m^2 | Full area of orifice | | y_max | float | m | Full height of orifice | | z_o | float | m | Offset of bottom above upstream superjunction invert | + | oneway | bool | | Is the flow one-way only? (upstream to downstream) | |-------------+-------+------+------------------------------------------------------| weirs: pd.DataFrame (optional) @@ -409,8 +410,8 @@ def __init__(self, superlinks, superjunctions, self._ki = links['k'].values.astype(np.int64) self.start_nodes = self.superlinks['j_0'].values.astype(np.int64) self.end_nodes = self.superlinks['j_1'].values.astype(np.int64) - self._is_start = np.zeros(self._I.size, dtype=np.bool8) - self._is_end = np.zeros(self._I.size, dtype=np.bool8) + self._is_start = np.zeros(self._I.size, dtype=np.bool_) + self._is_end = np.zeros(self._I.size, dtype=np.bool_) self._is_start[self.start_nodes] = True self._is_end[self.end_nodes] = True self.middle_nodes = self._I[(~self._is_start) & (~self._is_end)] @@ -460,7 +461,7 @@ def __init__(self, superlinks, superjunctions, self._n_ik = links['roughness'].values.astype(np.float64) else: self._n_ik = links['n'].values.astype(np.float64) - self._ctrl = links['ctrl'].values.astype(np.bool8) + self._ctrl = links['ctrl'].values.astype(np.bool_) self._A_c_ik = links['A_c'].values.astype(np.float64) self._C_ik = links['C'].values.astype(np.float64) self._storage_type = superjunctions['storage'] @@ -521,6 +522,10 @@ def __init__(self, superlinks, superjunctions, self._g1_o = np.sqrt(self._Ao_max) self._g2_o = np.sqrt(self._Ao_max) self._g3_o = np.zeros(self.n_o) + if 'oneway' in self.orifices.columns: + self._unidir_o = self.orifices['oneway'].values.astype(np.bool_) + else: + self._unidir_o = np.zeros(self.n_o, dtype=np.bool_) self._Qo = np.zeros(self.n_o, dtype=np.float64) self._alpha_o = np.zeros(self.n_o, dtype=np.float64) self._beta_o = np.zeros(self.n_o, dtype=np.float64) @@ -642,7 +647,7 @@ def __init__(self, superlinks, superjunctions, self._E_Ik = np.zeros(self._I.size) self._D_Ik = np.zeros(self._I.size) # Forward recurrence relations - self._I_end = np.zeros(self._I.size, dtype=np.bool8) + self._I_end = np.zeros(self._I.size, dtype=np.bool_) self._I_end[self.end_nodes] = True self._I_1k = self.start_nodes self._I_2k = self.forward_I_I[self._I_1k] @@ -653,7 +658,7 @@ def __init__(self, superlinks, superjunctions, self._V_Ik = np.zeros(self._I.size) self._W_Ik = np.zeros(self._I.size) # Backward recurrence relations - self._I_start = np.zeros(self._I.size, dtype=np.bool8) + self._I_start = np.zeros(self._I.size, dtype=np.bool_) self._I_start[self.start_nodes] = True self._I_Np1k = self.end_nodes self._I_Nk = self.backward_I_I[self._I_Np1k] @@ -680,7 +685,7 @@ def __init__(self, superlinks, superjunctions, self.A = np.zeros((self.M, self.M)) self.b = np.zeros(self.M) self.D = np.zeros(self.M) - self.bc = self.superjunctions['bc'].values.astype(np.bool8) + self.bc = self.superjunctions['bc'].values.astype(np.bool_) if sparse: self.B = scipy.sparse.lil_matrix((self.M, self.n_o)) self.O = scipy.sparse.lil_matrix((self.M, self.M)) @@ -732,8 +737,8 @@ def __init__(self, superlinks, superjunctions, self._S_o_uk = _S_o_uk self._S_o_dk = _S_o_dk # Boundary indexers - self._link_start = np.zeros(self._ik.size, dtype=np.bool8) - self._link_end = np.zeros(self._ik.size, dtype=np.bool8) + self._link_start = np.zeros(self._ik.size, dtype=np.bool_) + self._link_end = np.zeros(self._ik.size, dtype=np.bool_) self._link_start[self._i_1k] = True self._link_end[self._i_nk] = True # End roughness @@ -1094,7 +1099,7 @@ def _configure_internals_variable(self, internal_links=4, mobile_elements=True): # xx[:, :] = np.vstack([np.linspace(j, i + j + k, njunctions) # for i, j, k in zip(dx_j, _dx_uk, _dx_dk)]) zz[:] = xx * _m.reshape(-1, 1) + _b0.reshape(-1, 1) - _fixed = np.ones(xx.shape, dtype=np.bool8) + _fixed = np.ones(xx.shape, dtype=np.bool_) _xc = None _zc = None c = None