Feature/circ pump branch model

src/pandapipes/component_models/abstract_models/circulation_pump.py

@@ -4,13 +4,12 @@
 
 import numpy as np
 
-from pandapipes.component_models.abstract_models.branch_wzerolength_models import \
-    BranchWZeroLengthComponent
-from pandapipes.component_models.component_toolbox import set_fixed_node_entries, \
-    get_mass_flow_at_nodes
-from pandapipes.idx_branch import D, AREA, ACTIVE
-from pandapipes.idx_node import PINIT
-from pandapipes.pf.pipeflow_setup import get_lookup
+from pandapipes.component_models.abstract_models.branch_wzerolength_models import BranchWZeroLengthComponent
+from pandapipes.component_models.component_toolbox import set_fixed_node_entries, standard_branch_wo_internals_result_lookup
+from pandapipes.idx_branch import D, AREA, PUMP_TYPE, CIRC, LOAD_VEC_BRANCHES_T, TO_NODE
+from pandapipes.idx_node import MDOTSLACKINIT, CIRC_PUMP_OCCURENCE, EXT_GRID_OCCURENCE
+from pandapipes.pf.pipeflow_setup import get_fluid
+from pandapipes.pf.result_extraction import extract_branch_results_without_internals
 
 try:
     import pandaplan.core.pplog as logging
@@ -40,7 +39,14 @@ def get_result_table(cls, net):
                 if False, returns columns as tuples also specifying the dtypes
         :rtype: (list, bool)
         """
-        return ["mdot_flow_kg_per_s", "deltap_bar"], True
+        if get_fluid(net).is_gas:
+            output = ["v_from_m_per_s", "v_to_m_per_s", "v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k",
+                      "t_to_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds", "lambda",
+                      "normfactor_from", "normfactor_to", "qext_w"]
+        else:
+            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "mdot_from_kg_per_s",
+                      "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds", "lambda", "qext_w"]
+        return output, True
 
     @classmethod
     def active_identifier(cls):
@@ -73,8 +79,8 @@ def create_pit_node_entries(cls, net, node_pit):
         # TODO: there should be a warning, if any p_bar value is not given or any of the types does
         #       not contain "p", as this should not be allowed for this component
         press = circ_pump_tbl.p_flow_bar.values
-        set_fixed_node_entries(net, node_pit, junction, circ_pump_tbl.type.values, press,
-                               circ_pump_tbl.t_flow_k.values, cls.get_connected_node_type())
+        set_fixed_node_entries(net, node_pit, junction, circ_pump_tbl.type.values, press, circ_pump_tbl.t_flow_k.values,
+                               cls.get_connected_node_type(), circ_pump=True)
         return circ_pump_tbl, press
 
     @classmethod
@@ -90,7 +96,39 @@ def create_pit_branch_entries(cls, net, branch_pit):
         circ_pump_pit = super().create_pit_branch_entries(net, branch_pit)
         circ_pump_pit[:, D] = 0.1
         circ_pump_pit[:, AREA] = circ_pump_pit[:, D] ** 2 * np.pi / 4
-        circ_pump_pit[:, ACTIVE] = False
+        circ_pump_pit[:, PUMP_TYPE] = CIRC
+        return circ_pump_pit
+
+    @classmethod
+    def adaption_after_derivatives_hydraulic(cls, net, branch_pit, node_pit, idx_lookups, options):
+        """
+        Function which creates pit branch entries with a specific table.
+        :param net: The pandapipes network
+        :type net: pandapipesNet
+        :param branch_pit:
+        :type branch_pit:
+        :return: No Output.
+        """
+        f, t = idx_lookups[cls.table_name()]
+        circ_pump_pit = branch_pit[f:t, :]
+        tn = circ_pump_pit[:, TO_NODE].astype(np.int32)
+        mask = node_pit[tn, CIRC_PUMP_OCCURENCE] == node_pit[tn, EXT_GRID_OCCURENCE]
+        node_pit[tn[mask], MDOTSLACKINIT] = 0
+        return circ_pump_pit
+
+    @classmethod
+    def adaption_after_derivatives_thermal(cls, net, branch_pit, node_pit, idx_lookups, options):
+        """
+        Function which creates pit branch entries with a specific table.
+        :param net: The pandapipes network
+        :type net: pandapipesNet
+        :param branch_pit:
+        :type branch_pit:
+        :return: No Output.
+        """
+        f, t = idx_lookups[cls.table_name()]
+        circ_pump_pit = branch_pit[f:t, :]
+        circ_pump_pit[:, LOAD_VEC_BRANCHES_T] = 0
 
     @classmethod
     def extract_results(cls, net, options, branch_results, mode):
@@ -107,34 +145,8 @@ def extract_results(cls, net, options, branch_results, mode):
         :type options:
         :return: No Output.
         """
-        circ_pump_tbl = net[cls.table_name()]
-
-        if len(circ_pump_tbl) == 0:
-            return
-
-        res_table = net["res_" + cls.table_name()]
-
-        branch_pit = net['_pit']['branch']
-        node_pit = net["_pit"]["node"]
-
-        junction_lookup = get_lookup(net, "node", "index")[
-            cls.get_connected_node_type().table_name()]
-        fn_col, tn_col = cls.from_to_node_cols()
-        # get indices in internal structure for flow_junctions in circ_pump tables which are
-        # "active"
-        flow_junctions = circ_pump_tbl[tn_col].values
-        flow_nodes = junction_lookup[flow_junctions]
-        in_service = circ_pump_tbl.in_service.values
-        p_grids = np.isin(circ_pump_tbl.type.values, ["p", "pt"]) & in_service
-        sum_mass_flows, inverse_nodes, counts = get_mass_flow_at_nodes(net, node_pit, branch_pit,
-                                                                       flow_nodes[p_grids], cls)
-
-        # positive results mean that the circ_pump feeds in, negative means that the ext grid
-        # extracts (like a load)
-        res_table["mdot_flow_kg_per_s"].values[p_grids] = - (sum_mass_flows / counts)[inverse_nodes]
 
-        return_junctions = circ_pump_tbl[fn_col].values
-        return_nodes = junction_lookup[return_junctions]
+        required_results_hyd, required_results_ht = standard_branch_wo_internals_result_lookup(net)
 
-        deltap_bar = node_pit[flow_nodes, PINIT] - node_pit[return_nodes, PINIT]
-        res_table["deltap_bar"].values[in_service] = deltap_bar[in_service]
+        extract_branch_results_without_internals(net, branch_results, required_results_hyd, required_results_ht,
+                                                 cls.table_name(), mode)

src/pandapipes/component_models/circulation_pump_mass_component.py

@@ -2,15 +2,12 @@
 # and Energy System Technology (IEE), Kassel, and University of Kassel. All rights reserved.
 # Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
 
-import numpy as np
 from numpy import dtype
 
-from pandapipes.component_models.junction_component import Junction
 from pandapipes.component_models.abstract_models.circulation_pump import CirculationPump
-from pandapipes.idx_node import LOAD
-from pandapipes.pf.internals_toolbox import _sum_by_group
-from pandapipes.pf.pipeflow_setup import get_lookup
-from pandapipes.pf.pipeflow_setup import get_net_option
+from pandapipes.component_models.junction_component import Junction
+from pandapipes.idx_branch import JAC_DERIV_DP, JAC_DERIV_DP1, JAC_DERIV_DM, MDOTINIT, \
+    LOAD_VEC_BRANCHES
 
 try:
     import pandaplan.core.pplog as logging
@@ -51,24 +48,17 @@ def active_identifier(cls):
         return "in_service"
 
     @classmethod
-    def create_pit_node_entries(cls, net, node_pit):
-        """
-        Function which creates pit node entries.
+    def create_pit_branch_entries(cls, net, branch_pit):
+        circ_pump_pit = super().create_pit_branch_entries(net, branch_pit)
+        circ_pump_pit[:, MDOTINIT] = net[cls.table_name()].mdot_flow_kg_per_s.values
 
-        :param net: The pandapipes network
-        :type net: pandapipesNet
-        :param node_pit:
-        :type node_pit:
-        :return: No Output.
-        """
-        circ_pump, _ = super().create_pit_node_entries(net, node_pit)
-
-        mf = np.nan_to_num(circ_pump.mdot_flow_kg_per_s.values)
-        mass_flow_loads = mf * circ_pump.in_service.values
-        juncts, loads_sum = _sum_by_group(get_net_option(net, "use_numba"),
-                                          circ_pump.return_junction.values, mass_flow_loads)
-        junction_idx_lookups = get_lookup(net, "node", "index")[
-            cls.get_connected_node_type().table_name()]
-        index = junction_idx_lookups[juncts]
-        node_pit[index, LOAD] += loads_sum
+    @classmethod
+    def adaption_after_derivatives_hydraulic(cls, net, branch_pit, node_pit, idx_lookups, options):
+        # set all pressure derivatives to 0 and velocity to 1; load vector must be 0, as no change
+        # of velocity is allowed during the pipeflow iteration
+        circ_pump_pit = super().adaption_after_derivatives_hydraulic(net, branch_pit, node_pit, idx_lookups, options)
+        circ_pump_pit[:, JAC_DERIV_DP] = 0
+        circ_pump_pit[:, JAC_DERIV_DP1] = 0
+        circ_pump_pit[:, JAC_DERIV_DM] = 1
+        circ_pump_pit[:, LOAD_VEC_BRANCHES] = 0
 

src/pandapipes/component_models/circulation_pump_pressure_component.py

@@ -5,8 +5,8 @@
 from numpy import dtype
 
 from pandapipes.component_models.abstract_models.circulation_pump import CirculationPump
-from pandapipes.component_models.component_toolbox import set_fixed_node_entries
 from pandapipes.component_models.junction_component import Junction
+from pandapipes.idx_branch import JAC_DERIV_DP, JAC_DERIV_DP1, PL
 
 try:
     import pandaplan.core.pplog as logging
@@ -29,14 +29,8 @@ def get_component_input(cls):
         :return:
         :rtype:
         """
-        return [("name", dtype(object)),
-                ("return_junction", "u4"),
-                ("flow_junction", "u4"),
-                ("p_flow_bar", "f8"),
-                ("t_flow_k", "f8"),
-                ("plift_bar", "f8"),
-                ("in_service", 'bool'),
-                ("type", dtype(object))]
+        return [("name", dtype(object)), ("return_junction", "u4"), ("flow_junction", "u4"), ("p_flow_bar", "f8"),
+                ("t_flow_k", "f8"), ("plift_bar", "f8"), ("in_service", 'bool'), ("type", dtype(object))]
 
     @classmethod
     def active_identifier(cls):
@@ -47,19 +41,23 @@ def get_connected_node_type(cls):
         return Junction
 
     @classmethod
-    def create_pit_node_entries(cls, net, node_pit):
+    def create_pit_branch_entries(cls, net, branch_pit):
         """
-        Function which creates pit node entries.
-
+        Function which creates pit branch entries with a specific table.
         :param net: The pandapipes network
         :type net: pandapipesNet
-        :param node_pit:
-        :type node_pit:
+        :param branch_pit:
+        :type branch_pit:
         :return: No Output.
         """
-        circ_pump, press = super().create_pit_node_entries(net, node_pit)
+        circ_pump_pit = super().create_pit_branch_entries(net, branch_pit)
+        circ_pump_pit[:, PL] = net[cls.table_name()]['plift_bar'].values
+        return circ_pump_pit
 
-        junction = circ_pump[cls.from_to_node_cols()[0]].values
-        p_in = press - circ_pump.plift_bar.values
-        set_fixed_node_entries(net, node_pit, junction, circ_pump.type.values, p_in, None,
-                               cls.get_connected_node_type(), "p")
+    @classmethod
+    def adaption_after_derivatives_hydraulic(cls, net, branch_pit, node_pit, idx_lookups, options):
+        # set all pressure derivatives to 0 and velocity to 1; load vector must be 0, as no change
+        # of velocity is allowed during the pipeflow iteration
+        circ_pump_pit = super().adaption_after_derivatives_hydraulic(net, branch_pit, node_pit, idx_lookups, options)
+        circ_pump_pit[:, JAC_DERIV_DP] = 1
+        circ_pump_pit[:, JAC_DERIV_DP1] = -1

src/pandapipes/component_models/component_toolbox.py

@@ -6,11 +6,9 @@
 import pandas as pd
 
 from pandapipes import get_fluid
-from pandapipes.constants import NORMAL_PRESSURE, TEMP_GRADIENT_KPM, AVG_TEMPERATURE_K, \
-    HEIGHT_EXPONENT
-from pandapipes.idx_branch import LOAD_VEC_NODES, FROM_NODE, TO_NODE
-from pandapipes.idx_node import EXT_GRID_OCCURENCE, EXT_GRID_OCCURENCE_T
-from pandapipes.idx_node import PINIT, NODE_TYPE, P, TINIT, NODE_TYPE_T, T, LOAD
+from pandapipes.constants import NORMAL_PRESSURE, TEMP_GRADIENT_KPM, AVG_TEMPERATURE_K, HEIGHT_EXPONENT
+from pandapipes.idx_node import EXT_GRID_OCCURENCE, EXT_GRID_OCCURENCE_T, PINIT, NODE_TYPE, P, TINIT, NODE_TYPE_T, T, \
+    MDOTSLACKINIT, JAC_DERIV_MSL, CIRC_PUMP_OCCURENCE
 from pandapipes.pf.internals_toolbox import _sum_by_group
 from pandapipes.pf.pipeflow_setup import get_net_option, get_lookup
 
@@ -23,8 +21,7 @@ def p_correction_height_air(height):
     :return:
     :rtype:
     """
-    return NORMAL_PRESSURE * np.power(1 - height * TEMP_GRADIENT_KPM / AVG_TEMPERATURE_K,
-                                      HEIGHT_EXPONENT)
+    return NORMAL_PRESSURE * np.power(1 - height * TEMP_GRADIENT_KPM / AVG_TEMPERATURE_K, HEIGHT_EXPONENT)
 
 
 def vinterp(min_vals, max_vals, lengths):
@@ -84,7 +81,7 @@ def init_results_element(net, element, output, all_float):
     """
     res_element = "res_" + element
     if all_float:
-        net[res_element] = pd.DataFrame(np.NAN, columns=output, index=net[element].index,
+        net[res_element] = pd.DataFrame(np.nan, columns=output, index=net[element].index,
                                         dtype=np.float64)
     else:
         net[res_element] = pd.DataFrame(np.zeros(0, dtype=output), index=[])
@@ -125,8 +122,7 @@ def add_new_component(net, component, overwrite=False):
         if geodata is not None:
             net.update({name + '_geodata': geodata})
             if isinstance(net[name + '_geodata'], list):
-                net[name + '_geodata'] = pd.DataFrame(np.zeros(0, dtype=net[name + '_geodata']),
-                                                      index=[])
+                net[name + '_geodata'] = pd.DataFrame(np.zeros(0, dtype=net[name + '_geodata']), index=[])
 
 
 def set_entry_check_repeat(pit, column, entry, repeat_number, repeated=True):
@@ -136,11 +132,11 @@ def set_entry_check_repeat(pit, column, entry, repeat_number, repeated=True):
         pit[:, column] = entry
 
 
-def set_fixed_node_entries(net, node_pit, junctions, eg_types, p_values, t_values, node_comp,
-                           mode="sequential"):
+def set_fixed_node_entries(net, node_pit, junctions, eg_types, p_values, t_values, node_comp, mode="sequential",
+                           circ_pump=False):
     junction_idx_lookups = get_lookup(net, "node", "index")[node_comp.table_name()]
     for eg_type in ("p", "t"):
-        if eg_type not in mode and mode != "sequential" and mode!= "bidrectional":
+        if eg_type not in mode and mode != "sequential" and mode != "bidrectional":
             continue
         if eg_type == "p":
             val_col, type_col, eg_count_col, typ, valid_types, values = \
@@ -159,41 +155,22 @@ def set_fixed_node_entries(net, node_pit, junctions, eg_types, p_values, t_value
                                     + press_sum) / (number + node_pit[index, eg_count_col])
         node_pit[index, type_col] = typ
         node_pit[index, eg_count_col] += number
-
-
-def get_mass_flow_at_nodes(net, node_pit, branch_pit, eg_nodes, comp):
-    node_uni, inverse_nodes, counts = np.unique(eg_nodes, return_counts=True, return_inverse=True)
-    eg_from_branches = np.isin(branch_pit[:, FROM_NODE], node_uni)
-    eg_to_branches = np.isin(branch_pit[:, TO_NODE], node_uni)
-    from_nodes = branch_pit[eg_from_branches, FROM_NODE]
-    to_nodes = branch_pit[eg_to_branches, TO_NODE]
-    mass_flow_from = branch_pit[eg_from_branches, LOAD_VEC_NODES]
-    mass_flow_to = branch_pit[eg_to_branches, LOAD_VEC_NODES]
-    loads = node_pit[node_uni, LOAD]
-    all_index_nodes = np.concatenate([from_nodes, to_nodes, node_uni])
-    all_mass_flows = np.concatenate([-mass_flow_from, mass_flow_to, -loads])
-    nodes, sum_mass_flows = _sum_by_group(get_net_option(net, "use_numba"), all_index_nodes,
-                                          all_mass_flows)
-    if not np.all(nodes == node_uni):
-        raise UserWarning("In component %s: Something went wrong with the mass flow balance. "
-                          "Please report this error at github." % comp.__name__)
-    return sum_mass_flows, inverse_nodes, counts
+        if eg_type == 'p':
+            node_pit[index, JAC_DERIV_MSL] = -1.
+        if circ_pump and eg_type == 'p':
+            node_pit[index, CIRC_PUMP_OCCURENCE] += number
 
 
 def standard_branch_wo_internals_result_lookup(net):
-    required_results_hyd = [
-        ("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_to_kg_per_s", "mf_to"),
-        ("mdot_from_kg_per_s", "mf_from"), ("vdot_norm_m3_per_s", "vf"), ("lambda", "lambda"),
-        ("reynolds", "reynolds")
-    ]
+    required_results_hyd = [("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_to_kg_per_s", "mf_to"),
+                            ("mdot_from_kg_per_s", "mf_from"), ("vdot_norm_m3_per_s", "vf"), ("lambda", "lambda"),
+                            ("reynolds", "reynolds")]
     required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to")]
 
     if get_fluid(net).is_gas:
-        required_results_hyd.extend([
-            ("v_from_m_per_s", "v_gas_from"), ("v_to_m_per_s", "v_gas_to"),
-            ("v_mean_m_per_s", "v_gas_mean"), ("normfactor_from", "normfactor_from"),
-            ("normfactor_to", "normfactor_to")
-        ])
+        required_results_hyd.extend(
+            [("v_from_m_per_s", "v_gas_from"), ("v_to_m_per_s", "v_gas_to"), ("v_mean_m_per_s", "v_gas_mean"),
+             ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to")])
     else:
         required_results_hyd.extend([("v_mean_m_per_s", "v_mps")])
 
@@ -218,5 +195,5 @@ def get_component_array(net, component_name, component_type="branch", mode='hydr
     if not only_active:
         return net["_pit"]["components"][component_name]
     f_all, t_all = get_lookup(net, component_type, "from_to")[component_name]
-    in_service_elm = get_lookup(net, component_type, "active_%s"%mode)[f_all:t_all]
+    in_service_elm = get_lookup(net, component_type, "active_%s" % mode)[f_all:t_all]
     return net["_pit"]["components"][component_name][in_service_elm]

src/pandapipes/component_models/ext_grid_component.py

@@ -6,9 +6,9 @@
 from numpy import dtype
 
 from pandapipes.component_models.abstract_models.node_element_models import NodeElementComponent
-from pandapipes.component_models.component_toolbox import set_fixed_node_entries, \
-    get_mass_flow_at_nodes
+from pandapipes.component_models.component_toolbox import set_fixed_node_entries
 from pandapipes.pf.pipeflow_setup import get_lookup
+from pandapipes.idx_node import MDOTSLACKINIT
 
 try:
     import pandaplan.core.pplog as logging
@@ -91,8 +91,8 @@ def extract_results(cls, net, options, branch_results, mode):
         # get indices in internal structure for junctions in ext_grid tables which are "active"
         eg_nodes = get_lookup(net, "node", "index")[cls.get_connected_node_type().table_name()][
             junction[p_grids]]
-        sum_mass_flows, inverse_nodes, counts = get_mass_flow_at_nodes(net, node_pit, branch_pit,
-                                                                       eg_nodes, cls)
+        node_uni, inverse_nodes, counts = np.unique(eg_nodes, return_counts=True, return_inverse=True)
+        sum_mass_flows = node_pit[node_uni, MDOTSLACKINIT]
 
         # positive results mean that the ext_grid feeds in, negative means that the ext grid
         # extracts (like a load)