Skip to content

Commit d381c66

Browse files
mdbartosmdbartos
authored
Merge pull request #67 from mdbartos/wds_ctrl_structures
Fix control structures for water distribution networks
2 parents ee07af3 + b4c598c commit d381c66

2 files changed

Lines changed: 33 additions & 21 deletions

File tree

‎pipedream_solver/nsuperlink.py‎

Lines changed: 19 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -932,12 +932,13 @@ def orifice_flow_coefficients(self, u=None):
932932
_alpha_o = self._alpha_o # Orifice flow coefficient alpha_o
933933
_beta_o = self._beta_o # Orifice flow coefficient beta_o
934934
_chi_o = self._chi_o # Orifice flow coefficient chi_o
935+
_unidir_o = self._unidir_o
935936
# If no input signal, assume orifice is closed
936937
if u is None:
937938
u = np.zeros(self.n_o, dtype=np.float64)
938939
# Specify orifice heads at previous timestep
939940
numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_inv_j,
940-
_z_o, _tau_o, _Co, _Ao, _y_max_o, _J_uo, _J_do)
941+
_z_o, _tau_o, _Co, _Ao, _y_max_o, _unidir_o, _J_uo, _J_do)
941942
# Export instance variables
942943
self._alpha_o = _alpha_o
943944
self._beta_o = _beta_o
@@ -1424,13 +1425,14 @@ def solve_orifice_flows(self, dt, u=None):
14241425
_Co = self._Co # Discharge coefficient of orifice o
14251426
_Ao = self._Ao # Maximum flow area of orifice o
14261427
_V_sj = self._V_sj
1428+
_unidir_o = self._unidir_o
14271429
g = 9.81
14281430
# If no input signal, assume orifice is closed
14291431
if u is None:
14301432
u = np.zeros(self.n_o, dtype=np.float64)
14311433
# Compute orifice flows
14321434
_Qo_next = numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o, _tau_o, _y_max_o, _Co, _Ao,
1433-
_J_uo, _J_do, g)
1435+
_unidir_o, _J_uo, _J_do, g)
14341436
# TODO: Move this inside numba function
14351437
upstream_ctrl = (H_j[_J_uo] > H_j[_J_do])
14361438
_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):
25192521
return result
25202522

25212523
@njit(int64(float64[:], float64[:], float64[:], float64[:], float64[:], float64[:], float64[:],
2522-
float64[:], float64[:], float64[:], float64[:], float64[:], int64[:], int64[:]),
2524+
float64[:], float64[:], float64[:], float64[:], float64[:], boolean[:],
2525+
int64[:], int64[:]),
25232526
cache=True)
25242527
def numba_orifice_flow_coefficients(_alpha_o, _beta_o, _chi_o, H_j, _Qo, u, _z_inv_j,
2525-
_z_o, _tau_o, _Co, _Ao, _y_max_o, _J_uo, _J_do):
2528+
_z_o, _tau_o, _Co, _Ao, _y_max_o, _unidir_o,
2529+
_J_uo, _J_do):
25262530
g = 9.81
25272531
_H_uo = H_j[_J_uo]
25282532
_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
25392543
_z_o + _z_inv_uo + (_tau_o * _y_max_o * u / 2))
25402544
cond_2 = (_omega_o * _H_uo + (1 - _omega_o) * _H_do >
25412545
_z_o + _z_inv_uo)
2546+
cond_3 = (_H_do >= _H_uo) & _unidir_o
25422547
# Fill coefficient arrays
25432548
# Submerged on both sides
25442549
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
25592564
_chi_o[c] = (_gamma_o[c] * (-1)**(1 - _omega_o[c])
25602565
* (- _z_inv_uo[c] - _z_o[c]))
25612566
# No flow
2562-
d = (~cond_0 & ~cond_2)
2567+
d = (~cond_0 & ~cond_2) | cond_3
25632568
_alpha_o[d] = 0.0
25642569
_beta_o[d] = 0.0
25652570
_chi_o[d] = 0.0
25662571
return 1
25672572

25682573
@njit(float64[:](float64[:], float64[:], float64[:], float64[:],
2569-
float64[:], float64[:], float64[:], float64[:], int64[:], int64[:], float64),
2574+
float64[:], float64[:], float64[:], float64[:], boolean[:],
2575+
int64[:], int64[:], float64),
25702576
cache=True)
25712577
def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o,
2572-
_tau_o, _y_max_o, _Co, _Ao, _J_uo, _J_do, g=9.81):
2578+
_tau_o, _y_max_o, _Co, _Ao, _unidir_o, _J_uo, _J_do, g=9.81):
25732579
# Specify orifice heads at previous timestep
25742580
_H_uo = H_j[_J_uo]
25752581
_H_do = H_j[_J_do]
@@ -2590,6 +2596,7 @@ def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o,
25902596
_z_o + _z_inv_uo + (_tau_o * _y_max_o * u / 2))
25912597
cond_2 = (_omega_o * _H_uo + (1 - _omega_o) * _H_do >
25922598
_z_o + _z_inv_uo)
2599+
cond_3 = (_H_do >= _H_uo) & _unidir_o
25932600
# Fill coefficient arrays
25942601
# Submerged on both sides
25952602
a = (cond_0 & cond_1)
@@ -2610,7 +2617,7 @@ def numba_solve_orifice_flows(H_j, u, _z_inv_j, _z_o,
26102617
_chi_o[c] = (_gamma_o[c] * (-1)**(1 - _omega_o[c])
26112618
* (- _z_inv_uo[c] - _z_o[c]))
26122619
# No flow
2613-
d = (~cond_0 & ~cond_2)
2620+
d = (~cond_0 & ~cond_2) | cond_3
26142621
_alpha_o[d] = 0.0
26152622
_beta_o[d] = 0.0
26162623
_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,
27162723
cond_0 = _H_up > _z_inv_up + _z_p
27172724
# Condition 1: Head difference is within range of pump curve
27182725
cond_1 = (_dHp > _dHp_min) & (_dHp < _dHp_max)
2719-
_dHp[_dHp > _dHp_max] = _dHp_max
2720-
_dHp[_dHp < _dHp_min] = _dHp_min
2726+
_dHp[_dHp > _dHp_max] = _dHp_max[_dHp > _dHp_max]
2727+
_dHp[_dHp < _dHp_min] = _dHp_min[_dHp < _dHp_min]
27212728
# Compute universal coefficients
27222729
_gamma_p = gamma_p(_Qp, _b_p, _c_p, u)
27232730
# 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_
27482755
_z_inv_up = _z_inv_j[_J_up]
27492756
# Create conditionals
27502757
_dHp = _H_dp - _H_up
2751-
_dHp[_dHp > _dHp_max] = _dHp_max
2752-
_dHp[_dHp < _dHp_min] = _dHp_min
2758+
_dHp[_dHp > _dHp_max] = _dHp_max[_dHp > _dHp_max]
2759+
_dHp[_dHp < _dHp_min] = _dHp_min[_dHp < _dHp_min]
27532760
cond_0 = _H_up > _z_inv_up + _z_p
27542761
# Compute universal coefficients
27552762
_Qp_next = (u / _b_p * (_a_p - _dHp))**(1 / _c_p)

‎pipedream_solver/superlink.py‎

Lines changed: 14 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -173,6 +173,7 @@ class SuperLink():
173173
| A | float | m^2 | Full area of orifice |
174174
| y_max | float | m | Full height of orifice |
175175
| z_o | float | m | Offset of bottom above upstream superjunction invert |
176+
| oneway | bool | | Is the flow one-way only? (upstream to downstream) |
176177
|-------------+-------+------+------------------------------------------------------|
177178
178179
weirs: pd.DataFrame (optional)
@@ -409,8 +410,8 @@ def __init__(self, superlinks, superjunctions,
409410
self._ki = links['k'].values.astype(np.int64)
410411
self.start_nodes = self.superlinks['j_0'].values.astype(np.int64)
411412
self.end_nodes = self.superlinks['j_1'].values.astype(np.int64)
412-
self._is_start = np.zeros(self._I.size, dtype=np.bool8)
413-
self._is_end = np.zeros(self._I.size, dtype=np.bool8)
413+
self._is_start = np.zeros(self._I.size, dtype=np.bool_)
414+
self._is_end = np.zeros(self._I.size, dtype=np.bool_)
414415
self._is_start[self.start_nodes] = True
415416
self._is_end[self.end_nodes] = True
416417
self.middle_nodes = self._I[(~self._is_start) & (~self._is_end)]
@@ -460,7 +461,7 @@ def __init__(self, superlinks, superjunctions,
460461
self._n_ik = links['roughness'].values.astype(np.float64)
461462
else:
462463
self._n_ik = links['n'].values.astype(np.float64)
463-
self._ctrl = links['ctrl'].values.astype(np.bool8)
464+
self._ctrl = links['ctrl'].values.astype(np.bool_)
464465
self._A_c_ik = links['A_c'].values.astype(np.float64)
465466
self._C_ik = links['C'].values.astype(np.float64)
466467
self._storage_type = superjunctions['storage']
@@ -521,6 +522,10 @@ def __init__(self, superlinks, superjunctions,
521522
self._g1_o = np.sqrt(self._Ao_max)
522523
self._g2_o = np.sqrt(self._Ao_max)
523524
self._g3_o = np.zeros(self.n_o)
525+
if 'oneway' in self.orifices.columns:
526+
self._unidir_o = self.orifices['oneway'].values.astype(np.bool_)
527+
else:
528+
self._unidir_o = np.zeros(self.n_o, dtype=np.bool_)
524529
self._Qo = np.zeros(self.n_o, dtype=np.float64)
525530
self._alpha_o = np.zeros(self.n_o, dtype=np.float64)
526531
self._beta_o = np.zeros(self.n_o, dtype=np.float64)
@@ -642,7 +647,7 @@ def __init__(self, superlinks, superjunctions,
642647
self._E_Ik = np.zeros(self._I.size)
643648
self._D_Ik = np.zeros(self._I.size)
644649
# Forward recurrence relations
645-
self._I_end = np.zeros(self._I.size, dtype=np.bool8)
650+
self._I_end = np.zeros(self._I.size, dtype=np.bool_)
646651
self._I_end[self.end_nodes] = True
647652
self._I_1k = self.start_nodes
648653
self._I_2k = self.forward_I_I[self._I_1k]
@@ -653,7 +658,7 @@ def __init__(self, superlinks, superjunctions,
653658
self._V_Ik = np.zeros(self._I.size)
654659
self._W_Ik = np.zeros(self._I.size)
655660
# Backward recurrence relations
656-
self._I_start = np.zeros(self._I.size, dtype=np.bool8)
661+
self._I_start = np.zeros(self._I.size, dtype=np.bool_)
657662
self._I_start[self.start_nodes] = True
658663
self._I_Np1k = self.end_nodes
659664
self._I_Nk = self.backward_I_I[self._I_Np1k]
@@ -680,7 +685,7 @@ def __init__(self, superlinks, superjunctions,
680685
self.A = np.zeros((self.M, self.M))
681686
self.b = np.zeros(self.M)
682687
self.D = np.zeros(self.M)
683-
self.bc = self.superjunctions['bc'].values.astype(np.bool8)
688+
self.bc = self.superjunctions['bc'].values.astype(np.bool_)
684689
if sparse:
685690
self.B = scipy.sparse.lil_matrix((self.M, self.n_o))
686691
self.O = scipy.sparse.lil_matrix((self.M, self.M))
@@ -732,8 +737,8 @@ def __init__(self, superlinks, superjunctions,
732737
self._S_o_uk = _S_o_uk
733738
self._S_o_dk = _S_o_dk
734739
# Boundary indexers
735-
self._link_start = np.zeros(self._ik.size, dtype=np.bool8)
736-
self._link_end = np.zeros(self._ik.size, dtype=np.bool8)
740+
self._link_start = np.zeros(self._ik.size, dtype=np.bool_)
741+
self._link_end = np.zeros(self._ik.size, dtype=np.bool_)
737742
self._link_start[self._i_1k] = True
738743
self._link_end[self._i_nk] = True
739744
# End roughness
@@ -1094,7 +1099,7 @@ def _configure_internals_variable(self, internal_links=4, mobile_elements=True):
10941099
# xx[:, :] = np.vstack([np.linspace(j, i + j + k, njunctions)
10951100
# for i, j, k in zip(dx_j, _dx_uk, _dx_dk)])
10961101
zz[:] = xx * _m.reshape(-1, 1) + _b0.reshape(-1, 1)
1097-
_fixed = np.ones(xx.shape, dtype=np.bool8)
1102+
_fixed = np.ones(xx.shape, dtype=np.bool_)
10981103
_xc = None
10991104
_zc = None
11001105
c = None

0 commit comments

Comments
 (0)