Skip to content

Commit

Permalink
Merge pull request #67 from mdbartos/wds_ctrl_structures
Browse files Browse the repository at this point in the history 
Fix control structures for water distribution networks
  • Loading branch information
@mdbartos
mdbartos authored Nov 15, 2023
2 parents ee07af3 + b4c598c commit d381c66
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 21 deletions.
31 changes: 19 additions & 12 deletions pipedream_solver/nsuperlink.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand All @@ -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)
Expand All @@ -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]
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
23 changes: 14 additions & 9 deletions pipedream_solver/superlink.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)]
Expand Down Expand Up @@ -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']
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit d381c66

Please sign in to comment.