hydraulics.py 文件源码

def get_demand_or_head_residual(self, head, demand):

        if self.mode == 'PDD':
            minP = self.minimum_pressures
            nomP = self.nominal_pressures
            j_d = self.junction_demand
            m = self._slope_of_pdd_curve
            delta = self._pdd_smoothing_delta
            n_j = self.num_junctions
            P = head[:n_j] - self.node_elevations[:n_j]
            H = head[:n_j]
            Dact = demand[:n_j]

            self.demand_or_head_residual[:n_j] = (
                self.isolated_junction_array * H + (1.0 - self.isolated_junction_array)*(
                    (P <= minP) * (Dact - j_d*m*(P-minP)) +
                    (P > minP) * (P <= (minP + delta)) * (
                        Dact - j_d*(
                            self.pdd_poly1_coeffs_a*P**3 +
                            self.pdd_poly1_coeffs_b*P**2 +
                            self.pdd_poly1_coeffs_c*P +
                            self.pdd_poly1_coeffs_d
                        )
                    ) +
                    (P > (nomP - delta)) * (P <= nomP) * (
                        Dact - j_d*(
                            self.pdd_poly2_coeffs_a*P**3 +
                            self.pdd_poly2_coeffs_b*P**2 +
                            self.pdd_poly2_coeffs_c*P +
                            self.pdd_poly2_coeffs_d
                        )
                    ) +
                    (P > nomP) * (Dact - j_d * (m*(P-nomP) + 1.0))
                )
            )
            # for the last segment, assignment is required because 0*np.nan does not equal 0 (same with np.inf)
            last_segment = (Dact - j_d*((P-minP)/(nomP-minP))**0.5)
            last_segment[np.bitwise_not((P > (minP + delta))*(P <= (nomP - delta)))] = 0.0
            self.demand_or_head_residual[:n_j] = (self.demand_or_head_residual[:n_j] +
                                                  last_segment*(1.0-self.isolated_junction_array))
        else:
            self.demand_or_head_residual[:self.num_junctions] = (
                self.isolated_junction_array * head[:self.num_junctions] +
                (1.0 - self.isolated_junction_array) * (demand[:self.num_junctions] - self.junction_demand)
            )
        for node_id in self._tank_ids:
            self.demand_or_head_residual[node_id] = head[node_id] - self.tank_head[node_id]
        for node_id in self._reservoir_ids:
            self.demand_or_head_residual[node_id] = head[node_id] - self.reservoir_head[node_id]