diff --git a/.github/workflows/run_tests_develop.yml b/.github/workflows/run_tests_develop.yml
index cbe5ea65..e3a7b98d 100644
--- a/.github/workflows/run_tests_develop.yml
+++ b/.github/workflows/run_tests_develop.yml
@@ -179,9 +179,7 @@ jobs:
       - name: Check docs for Python ${{ matrix.python-version }}
         uses: e2nIEE/sphinx-action@master
         with:
-          pre-build-command: "apt-get update -y && apt-get install -y git;
-                              git --version;
-                              python -m pip install --upgrade pip;
-                              python -m pip install .[docs];"
+          pre-build-command: "python -m pip install --upgrade pip;
+                              python -m pip install .[docs]"
           build-command: "sphinx-build -b html source _build -W"
           docs-folder: "doc/"
diff --git a/.github/workflows/run_tests_master.yml b/.github/workflows/run_tests_master.yml
index d3be4c29..339bfc68 100644
--- a/.github/workflows/run_tests_master.yml
+++ b/.github/workflows/run_tests_master.yml
@@ -144,9 +144,7 @@ jobs:
       - name: Check docs for Python ${{ matrix.python-version }}
         uses: e2nIEE/sphinx-action@master
         with:
-          pre-build-command: "apt-get update -y && apt-get install -y git;
-                              git --version;
-                              python -m pip install --upgrade pip;
+          pre-build-command: "python -m pip install --upgrade pip;
                               python -m pip install .[docs]"
           build-command: "sphinx-build -b html source _build -W"
           docs-folder: "doc/"
diff --git a/src/pandapipes/component_models/__init__.py b/src/pandapipes/component_models/__init__.py
index 3b3f54c5..52d45480 100644
--- a/src/pandapipes/component_models/__init__.py
+++ b/src/pandapipes/component_models/__init__.py
@@ -17,3 +17,4 @@
 from pandapipes.component_models.flow_control_component import *
 from pandapipes.component_models.mass_storage_component import *
 from pandapipes.component_models.heat_consumer_component import *
+from pandapipes.component_models.component_toolbox import *
\ No newline at end of file
diff --git a/src/pandapipes/component_models/abstract_models/branch_wzerolength_models.py b/src/pandapipes/component_models/abstract_models/branch_wzerolength_models.py
index 31c013d4..ce6f0d71 100644
--- a/src/pandapipes/component_models/abstract_models/branch_wzerolength_models.py
+++ b/src/pandapipes/component_models/abstract_models/branch_wzerolength_models.py
@@ -2,9 +2,11 @@
 # 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 pandapipes.component_models.abstract_models.branch_wo_internals_models import \
     BranchWOInternalsComponent
-from pandapipes.idx_branch import LENGTH, K, TEXT, ALPHA
+from pandapipes.idx_branch import LENGTH, K, TEXT, ALPHA, D, AREA
 
 try:
     import pandaplan.core.pplog as logging
@@ -60,6 +62,8 @@ def create_pit_branch_entries(cls, net, branch_pit):
         branch_wzerolength_pit[:, K] = 1000
         branch_wzerolength_pit[:, TEXT] = 293.15
         branch_wzerolength_pit[:, ALPHA] = 0
+        branch_wzerolength_pit[:, D] = 0.1
+        branch_wzerolength_pit[:, AREA] = branch_wzerolength_pit[:, D] ** 2 * np.pi / 4
         return branch_wzerolength_pit
 
     @classmethod
diff --git a/src/pandapipes/component_models/abstract_models/circulation_pump.py b/src/pandapipes/component_models/abstract_models/circulation_pump.py
index c678b44a..5f7bfb4e 100644
--- a/src/pandapipes/component_models/abstract_models/circulation_pump.py
+++ b/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, BRANCH_TYPE, CIRC, LOAD_VEC_BRANCHES_T, TO_NODE
+from pandapipes.idx_node import MDOTSLACKINIT, VAR_MASS_SLACK, JAC_DERIV_MSL
+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
@@ -34,13 +33,23 @@ def get_component_input(cls):
     def get_result_table(cls, net):
         """
 
+        Gets the result table.
+
         :param net: The pandapipes network
         :type net: pandapipesNet
         :return: (columns, all_float) - the column names and whether they are all float type. Only
                 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", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s",
+                      "normfactor_from", "normfactor_to"]
+        else:
+            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s",
+                      "mdot_to_kg_per_s", "vdot_m3_per_s"]
+        output += ['deltat_k', 'qext_w']
+        return output, True
 
     @classmethod
     def active_identifier(cls):
@@ -72,10 +81,14 @@ 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())
-        return circ_pump_tbl, press
+        types = circ_pump_tbl.type.values
+        p_values = circ_pump_tbl.p_flow_bar.values
+        t_values = circ_pump_tbl.t_flow_k.values
+        index_p = set_fixed_node_entries(
+            net, node_pit, junction, types, p_values, cls.get_connected_node_type(), 'p')
+        set_fixed_node_entries(net, node_pit, junction, types, t_values, cls.get_connected_node_type(), 't')
+        node_pit[index_p, JAC_DERIV_MSL] = -1.
+        return circ_pump_tbl, p_values
 
     @classmethod
     def create_pit_branch_entries(cls, net, branch_pit):
@@ -90,7 +103,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[:, BRANCH_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, VAR_MASS_SLACK].astype(bool)
+        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 +152,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)
diff --git a/src/pandapipes/component_models/circulation_pump_mass_component.py b/src/pandapipes/component_models/circulation_pump_mass_component.py
index 39cdb217..64e3f872 100644
--- a/src/pandapipes/component_models/circulation_pump_mass_component.py
+++ b/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
 
diff --git a/src/pandapipes/component_models/circulation_pump_pressure_component.py b/src/pandapipes/component_models/circulation_pump_pressure_component.py
index fb7ff2af..9e77a549 100644
--- a/src/pandapipes/component_models/circulation_pump_pressure_component.py
+++ b/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
diff --git a/src/pandapipes/component_models/component_toolbox.py b/src/pandapipes/component_models/component_toolbox.py
index 4ba1324d..e4cbf92e 100644
--- a/src/pandapipes/component_models/component_toolbox.py
+++ b/src/pandapipes/component_models/component_toolbox.py
@@ -8,11 +8,11 @@
 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.pf.internals_toolbox import _sum_by_group
+from pandapipes.idx_branch import LOAD_VEC_NODES_FROM, LOAD_VEC_NODES_TO, FROM_NODE, TO_NODE
+from pandapipes.idx_node import (EXT_GRID_OCCURENCE, EXT_GRID_OCCURENCE_T,
+                                 PINIT, NODE_TYPE, P, TINIT, NODE_TYPE_T, T, LOAD)
 from pandapipes.pf.pipeflow_setup import get_net_option, get_lookup
+from pandapipes.pf.internals_toolbox import _sum_by_group
 
 
 def p_correction_height_air(height):
@@ -136,29 +136,36 @@ 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, types, values, node_comp, mode):
+    if not len(junctions):
+        return [], []
+
     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":
-            continue
-        if eg_type == "p":
-            val_col, type_col, eg_count_col, typ, valid_types, values = \
-                PINIT, NODE_TYPE, EXT_GRID_OCCURENCE, P, ["p", "pt"], p_values
-        else:
-            val_col, type_col, eg_count_col, typ, valid_types, values = \
-                TINIT, NODE_TYPE_T, EXT_GRID_OCCURENCE_T, T, ["t", "pt"], t_values
-        mask = np.isin(eg_types, valid_types)
-        if not np.any(mask):
-            continue
-        use_numba = get_net_option(net, "use_numba")
-        juncts, press_sum, number = _sum_by_group(use_numba, junctions[mask], values[mask],
-                                                  np.ones_like(values[mask], dtype=np.int32))
-        index = junction_idx_lookups[juncts]
-        node_pit[index, val_col] = (node_pit[index, val_col] * node_pit[index, eg_count_col]
-                                    + press_sum) / (number + node_pit[index, eg_count_col])
-        node_pit[index, type_col] = typ
-        node_pit[index, eg_count_col] += number
+    use_numba = get_net_option(net, "use_numba")
+
+    if mode == "p":
+        val_col, type_col, count_col, typ, valid_types, values = \
+            PINIT, NODE_TYPE, EXT_GRID_OCCURENCE, P, ["p", "pt"], values
+    elif mode == "t":
+        val_col, type_col, count_col, typ, valid_types, values = \
+            TINIT, NODE_TYPE_T, EXT_GRID_OCCURENCE_T, T, ["t", "pt"], values
+    else:
+        raise UserWarning(r'The mode %s is not supported. Choose either mode "p" or "t"' % mode)
+
+    mask = np.isin(types, valid_types)
+
+    juncts, val_sum, number = _sum_by_group(use_numba, junctions[mask], values[mask],
+                                            np.ones_like(values[mask], dtype=np.int32))
+
+    index = junction_idx_lookups[juncts]
+
+    node_pit[index, val_col] = (node_pit[index, val_col] * node_pit[index, count_col] + val_sum) / \
+                               (number + node_pit[index, count_col])
+
+    node_pit[index, count_col] += number
+    node_pit[index, type_col] = typ
+
+    return index
 
 
 def get_mass_flow_at_nodes(net, node_pit, branch_pit, eg_nodes, comp):
@@ -167,8 +174,8 @@ def get_mass_flow_at_nodes(net, node_pit, branch_pit, eg_nodes, comp):
     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]
+    mass_flow_from = branch_pit[eg_from_branches, LOAD_VEC_NODES_FROM]
+    mass_flow_to = branch_pit[eg_to_branches, LOAD_VEC_NODES_TO]
     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])
@@ -183,18 +190,17 @@ def get_mass_flow_at_nodes(net, node_pit, branch_pit, eg_nodes, comp):
 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"), ("lambda", "lambda"), ("reynolds", "reynolds")
+        ("mdot_from_kg_per_s", "mf_from")
     ]
     required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
 
     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_from", "normfactor_from"),
             ("normfactor_to", "normfactor_to"), ("vdot_norm_m3_per_s", "vf")
         ])
     else:
-        required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
+        required_results_hyd.extend([("vdot_m3_per_s", "vf")])
 
     return required_results_hyd, required_results_ht
 
diff --git a/src/pandapipes/component_models/compressor_component.py b/src/pandapipes/component_models/compressor_component.py
index 3d16b7b5..84c0e8a1 100644
--- a/src/pandapipes/component_models/compressor_component.py
+++ b/src/pandapipes/component_models/compressor_component.py
@@ -40,9 +40,6 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         compressor_pit = super(Pump, cls).create_pit_branch_entries(net, branch_pit)
-
-        compressor_pit[:, D] = 0.1
-        compressor_pit[:, AREA] = compressor_pit[:, D] ** 2 * np.pi / 4
         compressor_pit[:, LC] = 0
 
     @classmethod
diff --git a/src/pandapipes/component_models/ext_grid_component.py b/src/pandapipes/component_models/ext_grid_component.py
index 2e1d9e79..e88590a0 100644
--- a/src/pandapipes/component_models/ext_grid_component.py
+++ b/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, VAR_MASS_SLACK, JAC_DERIV_MSL
 
 try:
     import pandaplan.core.pplog as logging
@@ -55,11 +55,15 @@ def create_pit_node_entries(cls, net, node_pit):
         ext_grids = ext_grids[ext_grids[cls.active_identifier()].values]
 
         junction = ext_grids[cls.get_node_col()].values
-        press = ext_grids.p_bar.values
-        set_fixed_node_entries(net, node_pit, junction, ext_grids.type.values, press,
-                               ext_grids.t_k.values, cls.get_connected_node_type())
-
-        return ext_grids, press
+        types = ext_grids.type.values
+        p_values = ext_grids.p_bar.values
+        t_values = ext_grids.t_k.values
+        index_p = set_fixed_node_entries(
+            net, node_pit, junction, types, p_values, cls.get_connected_node_type(), 'p')
+        set_fixed_node_entries(net, node_pit, junction, types, t_values, cls.get_connected_node_type(), 't')
+        node_pit[index_p, JAC_DERIV_MSL] = -1.
+        node_pit[index_p, VAR_MASS_SLACK] = True
+        return ext_grids, p_values
 
     @classmethod
     def extract_results(cls, net, options, branch_results, mode):
@@ -91,8 +95,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)
diff --git a/src/pandapipes/component_models/flow_control_component.py b/src/pandapipes/component_models/flow_control_component.py
index e9cd4059..57e2d08f 100644
--- a/src/pandapipes/component_models/flow_control_component.py
+++ b/src/pandapipes/component_models/flow_control_component.py
@@ -49,8 +49,6 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         fc_branch_pit = super().create_pit_branch_entries(net, branch_pit)
-        fc_branch_pit[:, D] = net[cls.table_name()].diameter_m.values
-        fc_branch_pit[:, AREA] = fc_branch_pit[:, D] ** 2 * np.pi / 4
         fc_branch_pit[:, MDOTINIT] = net[cls.table_name()].controlled_mdot_kg_per_s.values
 
     @classmethod
@@ -105,7 +103,6 @@ def get_component_input(cls):
                 ("from_junction", "u4"),
                 ("to_junction", "u4"),
                 ("controlled_mdot_kg_per_s", "f8"),
-                ("diameter_m", "f8"),
                 ("control_active", "bool"),
                 ("in_service", 'bool'),
                 ("type", dtype(object))]
@@ -123,12 +120,12 @@ def get_result_table(cls, net):
         :rtype: (list, bool)
         """
         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",
+            output = ["p_from_bar", "p_to_bar",
                       "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s",
                       "vdot_norm_m3_per_s", "reynolds", "lambda", "normfactor_from",
                       "normfactor_to"]
         else:
-            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
+            output = ["p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
                       "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True
diff --git a/src/pandapipes/component_models/heat_consumer_component.py b/src/pandapipes/component_models/heat_consumer_component.py
index d3f574ef..425d9e6b 100644
--- a/src/pandapipes/component_models/heat_consumer_component.py
+++ b/src/pandapipes/component_models/heat_consumer_component.py
@@ -65,8 +65,6 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         hc_pit = super().create_pit_branch_entries(net, branch_pit)
-        hc_pit[:, D] = net[cls.table_name()].diameter_m.values
-        hc_pit[:, AREA] = hc_pit[:, D] ** 2 * np.pi / 4
         qext = net[cls.table_name()].qext_w.values
         hc_pit[~np.isnan(qext), QEXT] = qext[~np.isnan(qext)]
         mdot = net[cls.table_name()].controlled_mdot_kg_per_s.values
@@ -207,7 +205,7 @@ def get_component_input(cls):
         :rtype:
         """
         return [("name", dtype(object)), ("from_junction", "u4"), ("to_junction", "u4"), ("qext_w", "f8"),
-                ("controlled_mdot_kg_per_s", "f8"), ("deltat_k", "f8"), ("treturn_k", "f8"), ("diameter_m", "f8"),
+                ("controlled_mdot_kg_per_s", "f8"), ("deltat_k", "f8"), ("treturn_k", "f8"),
                 ("in_service", "bool"), ("type", dtype(object))]
 
     @classmethod
@@ -223,11 +221,11 @@ def get_result_table(cls, net):
         :rtype: (list, bool)
         """
         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",
+            output = ["p_from_bar", "p_to_bar", "t_from_k",
                       "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds", "lambda",
                       "normfactor_from", "normfactor_to"]
         else:
-            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s",
+            output = ["p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s",
                       "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds", "lambda"]
         output += ['deltat_k', 'qext_w']
         return output, True
diff --git a/src/pandapipes/component_models/heat_exchanger_component.py b/src/pandapipes/component_models/heat_exchanger_component.py
index 6cb75ca0..4bc48aec 100644
--- a/src/pandapipes/component_models/heat_exchanger_component.py
+++ b/src/pandapipes/component_models/heat_exchanger_component.py
@@ -52,8 +52,6 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         heat_exchanger_pit = super().create_pit_branch_entries(net, branch_pit)
-        heat_exchanger_pit[:, D] = net[cls.table_name()].diameter_m.values
-        heat_exchanger_pit[:, AREA] = heat_exchanger_pit[:, D] ** 2 * np.pi / 4
         heat_exchanger_pit[:, LC] = net[cls.table_name()].loss_coefficient.values
         heat_exchanger_pit[:, QEXT] = net[cls.table_name()].qext_w.values
 
@@ -105,12 +103,12 @@ def get_result_table(cls, net):
         :rtype: (list, bool)
         """
         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",
+            output = ["p_from_bar", "p_to_bar",
                       "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s",
                       "vdot_norm_m3_per_s", "reynolds", "lambda", "normfactor_from",
                       "normfactor_to"]
         else:
-            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
+            output = ["p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
                       "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True
diff --git a/src/pandapipes/component_models/pressure_control_component.py b/src/pandapipes/component_models/pressure_control_component.py
index 60fc69c1..20c5d2b3 100644
--- a/src/pandapipes/component_models/pressure_control_component.py
+++ b/src/pandapipes/component_models/pressure_control_component.py
@@ -7,10 +7,11 @@
 
 from pandapipes.component_models.abstract_models.branch_wzerolength_models import \
     BranchWZeroLengthComponent
+from pandapipes.component_models import standard_branch_wo_internals_result_lookup
 from pandapipes.component_models.junction_component import Junction
 from pandapipes.idx_branch import D, AREA, \
-    JAC_DERIV_DP, JAC_DERIV_DP1, JAC_DERIV_DM, BRANCH_TYPE, LOSS_COEFFICIENT as LC
-from pandapipes.idx_node import PINIT, NODE_TYPE, PC
+    JAC_DERIV_DP, JAC_DERIV_DP1, JAC_DERIV_DM, BRANCH_TYPE, LOSS_COEFFICIENT as LC, PC as PC_BRANCH
+from pandapipes.idx_node import PINIT, NODE_TYPE, PC as PC_NODE
 from pandapipes.pf.pipeflow_setup import get_lookup
 from pandapipes.pf.result_extraction import extract_branch_results_without_internals
 from pandapipes.properties.fluids import get_fluid
@@ -46,7 +47,7 @@ def create_pit_node_entries(cls, net, node_pit):
         junction_idx_lookups = get_lookup(net, "node", "index")[
             cls.get_connected_node_type().table_name()]
         index_pc = junction_idx_lookups[juncts]
-        node_pit[index_pc, NODE_TYPE] = PC
+        node_pit[index_pc, NODE_TYPE] = PC_NODE
         node_pit[index_pc, PINIT] = press
 
     @classmethod
@@ -60,9 +61,7 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         pc_pit = super().create_pit_branch_entries(net, branch_pit)
-        pc_pit[:, D] = 0.1
-        pc_pit[:, AREA] = pc_pit[:, D] ** 2 * np.pi / 4
-        pc_pit[net[cls.table_name()].control_active.values, BRANCH_TYPE] = PC
+        pc_pit[net[cls.table_name()].control_active.values, BRANCH_TYPE] = PC_BRANCH
         pc_pit[:, LC] = net[cls.table_name()].loss_coefficient.values
 
     @classmethod
@@ -74,7 +73,7 @@ def adaption_after_derivatives_hydraulic(cls, net, branch_pit, node_pit, idx_loo
         # set all PC branches to derivatives to 0
         f, t = idx_lookups[cls.table_name()]
         press_pit = branch_pit[f:t, :]
-        pc_branch = press_pit[:, BRANCH_TYPE] == PC
+        pc_branch = press_pit[:, BRANCH_TYPE] == PC_BRANCH
         press_pit[pc_branch, JAC_DERIV_DP] = 0
         press_pit[pc_branch, JAC_DERIV_DP1] = 0
         press_pit[pc_branch, JAC_DERIV_DM] = 0
@@ -94,20 +93,8 @@ def extract_results(cls, net, options, branch_results, mode):
         :type options:
         :return: No Output.
         """
-        required_results_hyd = [
-            ("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_from_kg_per_s", "mf_from"),
-            ("mdot_to_kg_per_s", "mf_to")
-        ]
-        required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
 
-        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"),
-                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to"),
-                ("vdot_norm_m3_per_s", "vf")
-            ])
-        else:
-            required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
+        required_results_hyd, required_results_ht = standard_branch_wo_internals_result_lookup(net)
 
         extract_branch_results_without_internals(net, branch_results, required_results_hyd,
                                                  required_results_ht, cls.table_name(), mode)
@@ -150,11 +137,11 @@ def get_result_table(cls, net):
         :rtype: (list, bool)
         """
         if get_fluid(net).is_gas:
-            output = ["v_from_m_per_s", "v_to_m_per_s", "p_from_bar", "p_to_bar",
+            output = ["p_from_bar", "p_to_bar",
                       "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s",
                       "vdot_norm_m3_per_s", "normfactor_from", "normfactor_to"]
         else:
-            output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
+            output = ["p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
                       "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s"]
         output += ["deltap_bar"]
         return output, True
diff --git a/src/pandapipes/component_models/pump_component.py b/src/pandapipes/component_models/pump_component.py
index 3f18cbbb..59cc44ae 100644
--- a/src/pandapipes/component_models/pump_component.py
+++ b/src/pandapipes/component_models/pump_component.py
@@ -9,6 +9,7 @@
 
 from pandapipes.component_models.abstract_models.branch_wzerolength_models import \
     BranchWZeroLengthComponent
+from pandapipes.component_models import standard_branch_wo_internals_result_lookup
 from pandapipes.component_models.component_toolbox import get_component_array
 from pandapipes.component_models.junction_component import Junction
 from pandapipes.constants import NORMAL_TEMPERATURE, NORMAL_PRESSURE, R_UNIVERSAL, P_CONVERSION
@@ -61,8 +62,6 @@ def create_pit_branch_entries(cls, net, branch_pit):
         :return: No Output.
         """
         pump_pit = super().create_pit_branch_entries(net, branch_pit)
-        pump_pit[:, D] = 0.1
-        pump_pit[:, AREA] = pump_pit[:, D] ** 2 * np.pi / 4
         pump_pit[:, LC] = 0
 
     @classmethod
@@ -134,26 +133,14 @@ def extract_results(cls, net, options, branch_results, mode):
         :type mode:
         :return: No Output.
         """
-        calc_compr_pow = options['calc_compression_power']
 
-        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"), ("deltap_bar", "pl"),
-        ]
-        required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
-
-        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"),
-                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to"),
-                ("vdot_norm_m3_per_s", "vf")
-            ])
-        else:
-            required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
+        required_results_hyd, required_results_ht = standard_branch_wo_internals_result_lookup(net)
+        required_results_hyd.extend([("deltap_bar", "pl")])
 
         extract_branch_results_without_internals(net, branch_results, required_results_hyd,
                                                  required_results_ht, cls.table_name(), mode)
 
+        calc_compr_pow = options['calc_compression_power']
         if calc_compr_pow:
             f, t = get_lookup(net, "branch", "from_to")[cls.table_name()]
             from_nodes = branch_results["from_nodes"][f:t]
@@ -218,13 +205,11 @@ def get_result_table(cls, net):
 
         if get_fluid(net).is_gas:
             output = ["deltap_bar",
-                      "v_from_m_per_s", "v_to_m_per_s",
                       "p_from_bar", "p_to_bar",
                       "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s",
                       "vdot_norm_m3_per_s", "normfactor_from", "normfactor_to"]
-            # TODO: inwieweit sind diese Angaben bei imaginärem Durchmesser sinnvoll?
         else:
-            output = ["deltap_bar", "v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k",
+            output = ["deltap_bar", "p_from_bar", "p_to_bar", "t_from_k",
                       "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s"]
         if calc_compr_pow:
             output += ["compr_power_mw"]
diff --git a/src/pandapipes/component_models/valve_component.py b/src/pandapipes/component_models/valve_component.py
index 1fffac7f..5932d7ab 100644
--- a/src/pandapipes/component_models/valve_component.py
+++ b/src/pandapipes/component_models/valve_component.py
@@ -71,6 +71,11 @@ def get_component_input(cls):
     def extract_results(cls, net, options, branch_results, mode):
         required_results_hyd, required_results_ht = standard_branch_wo_internals_result_lookup(net)
 
+        required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("lambda", "lambda"), ("reynolds", "reynolds")])
+
+        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")])
+
         extract_branch_results_without_internals(net, branch_results, required_results_hyd,
                                                  required_results_ht, cls.table_name(), mode)
 
diff --git a/src/pandapipes/create.py b/src/pandapipes/create.py
index c3347359..dbf233f6 100644
--- a/src/pandapipes/create.py
+++ b/src/pandapipes/create.py
@@ -325,7 +325,7 @@ def create_ext_grid(net, junction, p_bar=None, t_k=None, type="auto", name=None,
     return index
 
 
-def create_heat_exchanger(net, from_junction, to_junction, diameter_m, qext_w, loss_coefficient=0, name=None,
+def create_heat_exchanger(net, from_junction, to_junction, qext_w, loss_coefficient=0, name=None,
                           index=None, in_service=True, type="heat_exchanger", **kwargs):
     """
     Creates a heat exchanger element in net["heat_exchanger"] from heat exchanger parameters.
@@ -338,8 +338,6 @@ def create_heat_exchanger(net, from_junction, to_junction, diameter_m, qext_w, l
     :param to_junction: ID of the junction on the other side which the heat exchanger will be\
             connected with
     :type to_junction: int
-    :param diameter_m: The heat exchanger inner diameter in [m]
-    :type diameter_m: float
     :param qext_w: External heat flux in [W]. If positive, heat is derived from the network. If
             negative, heat is being fed into the network from a heat source.
     :type qext_w: float
@@ -360,15 +358,18 @@ def create_heat_exchanger(net, from_junction, to_junction, diameter_m, qext_w, l
     :rtype: int
 
     :Example:
-        >>> create_heat_exchanger(net, from_junction=0, to_junction=1,
-        >>>                       diameter_m=40e-3, qext_w=2000)
+        >>> create_heat_exchanger(net, from_junction=0, to_junction=1, qext_w=2000)
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the compoent table for postprocessing purposes by you if required.')
+
     add_new_component(net, HeatExchanger)
 
     index = _get_index_with_check(net, "heat_exchanger", index, "heat exchanger")
     _check_branch(net, "Heat exchanger", index, from_junction, to_junction)
 
-    v = {"name": name, "from_junction": from_junction, "to_junction": to_junction, "diameter_m": diameter_m,
+    v = {"name": name, "from_junction": from_junction, "to_junction": to_junction,
          "qext_w": qext_w, "loss_coefficient": loss_coefficient, "in_service": bool(in_service), "type": type}
     _set_entries(net, "heat_exchanger", index, **v, **kwargs)
 
@@ -959,7 +960,7 @@ def create_pressure_control(net, from_junction, to_junction, controlled_junction
     return index
 
 
-def create_flow_control(net, from_junction, to_junction, controlled_mdot_kg_per_s, diameter_m, control_active=True,
+def create_flow_control(net, from_junction, to_junction, controlled_mdot_kg_per_s, control_active=True,
                         name=None, index=None, in_service=True, type="fc", **kwargs):
     """
     Adds one flow control with a constant mass flow in table net["flow_control"].
@@ -974,9 +975,6 @@ def create_flow_control(net, from_junction, to_junction, controlled_mdot_kg_per_
     :type to_junction: int
     :param controlled_mdot_kg_per_s: Mass flow set point
     :type controlled_mdot_kg_per_s: float
-    :param diameter_m: Measure of the diameter to derive the cross-sectional area (important for \
-            the velocity calculation)
-    :type diameter_m: float
     :param control_active: Variable to state whether the flow control is active (otherwise \
             similar to open valve)
     :type control_active: bool, default True
@@ -998,6 +996,9 @@ def create_flow_control(net, from_junction, to_junction, controlled_mdot_kg_per_
         >>> create_flow_control(net, 0, 1, 0.5, 0.8)
 
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the compoent table for postprocessing purposes by you if required.')
 
     add_new_component(net, FlowControlComponent)
 
@@ -1007,13 +1008,13 @@ def create_flow_control(net, from_junction, to_junction, controlled_mdot_kg_per_
     _check_branch(net, "FlowControl", index, from_junction, to_junction)
 
     _set_entries(net, "flow_control", index, name=name, from_junction=from_junction, to_junction=to_junction,
-                 controlled_mdot_kg_per_s=controlled_mdot_kg_per_s, diameter_m=diameter_m,
+                 controlled_mdot_kg_per_s=controlled_mdot_kg_per_s,
                  control_active=bool(control_active), in_service=bool(in_service), type=type, **kwargs)
 
     return index
 
 
-def create_heat_consumer(net, from_junction, to_junction, diameter_m, qext_w=None, controlled_mdot_kg_per_s=None,
+def create_heat_consumer(net, from_junction, to_junction, qext_w=None, controlled_mdot_kg_per_s=None,
                          deltat_k=None, treturn_k=None, name=None, index=None, in_service=True, type="heat_consumer",
                          **kwargs):
     """
@@ -1027,8 +1028,6 @@ def create_heat_consumer(net, from_junction, to_junction, diameter_m, qext_w=Non
     :param to_junction: ID of the junction on the other side which the heat consumer will be \
         connected with
     :type to_junction: int
-    :param diameter_m: The heat consumer inner diameter in [m] - only for result calculation
-    :type diameter_m: float
     :param qext_w: External heat flux in [W]. If positive, heat is extracted from the network. If \
         negative, heat is being fed into the network from a heat source.
     :type qext_w: float, default None
@@ -1057,6 +1056,10 @@ def create_heat_consumer(net, from_junction, to_junction, diameter_m, qext_w=Non
         >>> create_heat_consumer(net,from_junction=0, to_junction=1, diameter_m=40e-3, qext_w=20000,
         >>>                     controlled_mdot_kg_per_s=0.4, name="heat_consumer1")
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the compoent table for postprocessing purposes by you if required.')
+
     if ((controlled_mdot_kg_per_s is None) + (qext_w is None) + (deltat_k is None) + (treturn_k is None) != 2):
         raise AttributeError(r"Define exactly two varibales from 'controlled_mdot_kg_per_s', "
                              r"'qext_w' and 'deltat_k' or 'treturn_k' different from None.")
@@ -1069,7 +1072,7 @@ def create_heat_consumer(net, from_junction, to_junction, diameter_m, qext_w=Non
     index = _get_index_with_check(net, "heat_consumer", index, "heat consumer")
     _check_branch(net, "Heat consumer", index, from_junction, to_junction)
 
-    v = {"name": name, "from_junction": from_junction, "to_junction": to_junction, "diameter_m": diameter_m,
+    v = {"name": name, "from_junction": from_junction, "to_junction": to_junction,
          "qext_w": qext_w, "controlled_mdot_kg_per_s": controlled_mdot_kg_per_s, "deltat_k": deltat_k,
          "treturn_k": treturn_k, "in_service": bool(in_service), "type": type}
     _set_entries(net, "heat_consumer", index, **v, **kwargs)
@@ -1582,7 +1585,7 @@ def create_pressure_controls(net, from_junctions, to_junctions, controlled_junct
     return index
 
 
-def create_flow_controls(net, from_junctions, to_junctions, controlled_mdot_kg_per_s, diameter_m, control_active=True,
+def create_flow_controls(net, from_junctions, to_junctions, controlled_mdot_kg_per_s, control_active=True,
                          name=None, index=None, in_service=True, type="fc", **kwargs):
     """
     Convenience function for creating many flow controls at once. Parameters 'from_junctions'\
@@ -1599,9 +1602,6 @@ def create_flow_controls(net, from_junctions, to_junctions, controlled_mdot_kg_p
     :type to_junctions: Iterable(int)
     :param controlled_mdot_kg_per_s: Mass flow set points
     :type controlled_mdot_kg_per_s: Iterable or float
-    :param diameter_m: Measure of the diameter to derive the cross-sectional area (important for \
-            the velocity calculation)
-    :type diameter_m: Iterable or float
     :param control_active: Variable to state whether the flow control is active (otherwise \
             similar to open valve)
     :type control_active: bool, default True
@@ -1620,12 +1620,15 @@ def create_flow_controls(net, from_junctions, to_junctions, controlled_mdot_kg_p
     :rtype: array(int)
 
     :Example:
-        Create two flow controllers with 0.8 m diameter between junction 0 and 1 with 0.5 kg/s and
+        Create two flow controllers between junction 0 and 1 with 0.5 kg/s and
         junction 2 and 4 with 0.9 kg/s, respectively.
 
-        >>> create_flow_controls(net, [0, 2], [1, 4], [0.5, 0.9], [0.8, 0.8])
+        >>> create_flow_controls(net, [0, 2], [1, 4], [0.5, 0.9])
 
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the compoent table for postprocessing purposes by you if required.')
 
     add_new_component(net, FlowControlComponent)
 
@@ -1633,14 +1636,14 @@ def create_flow_controls(net, from_junctions, to_junctions, controlled_mdot_kg_p
     _check_branches(net, from_junctions, to_junctions, "flow_control")
 
     entries = {"name": name, "from_junction": from_junctions, "to_junction": to_junctions,
-               "controlled_mdot_kg_per_s": controlled_mdot_kg_per_s, "diameter_m": diameter_m,
+               "controlled_mdot_kg_per_s": controlled_mdot_kg_per_s,
                "control_active": control_active, "in_service": in_service, "type": type}
     _set_multiple_entries(net, "flow_control", index, **entries, **kwargs)
 
     return index
 
 
-def create_heat_exchangers(net, from_junctions, to_junctions, diameter_m, qext_w, loss_coefficient=0, name=None,
+def create_heat_exchangers(net, from_junctions, to_junctions, qext_w, loss_coefficient=0, name=None,
                            index=None, in_service=True, type="heat_exchanger", **kwargs):
     """
     Convenience function for creating many heat exchangers at once. Parameters 'from_junctions'\
@@ -1655,8 +1658,6 @@ def create_heat_exchangers(net, from_junctions, to_junctions, diameter_m, qext_w
     :param to_junctions: IDs of the junctions on the other side the heat exchangers will be\
             connected with
     :type to_junctions: Iterable(int)
-    :param diameter_m: The heat exchangers inner diameter in [m]
-    :type diameter_m: Iterable(float) or float
     :param qext_w: External heat flux in [W]. If positive, heat is derived from the network. If
             negative, heat is being fed into the network from a heat source.
     :type qext_w: Iterable(float) or float
@@ -1679,22 +1680,25 @@ def create_heat_exchangers(net, from_junctions, to_junctions, diameter_m, qext_w
     :rtype: Iterable(int), default None
 
     :Example:
-        >>> create_heat_exchangers(net, from_junctions=[0,1], to_junctions=[2,3],
-        >>>                       diameter_m=40e-3, qext_w=2000)
+        >>> create_heat_exchangers(net, from_junctions=[0,1], to_junctions=[2,3], qext_w=2000)
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the component table for postprocessing purposes by you if required.')
+
     add_new_component(net, HeatExchanger)
 
     index = _get_multiple_index_with_check(net, "heat_exchanger", index, len(from_junctions))
     _check_branches(net, from_junctions, to_junctions, "heat_exchanger")
 
-    entries = {"name": name, "from_junction": from_junctions, "to_junction": to_junctions, "diameter_m": diameter_m,
+    entries = {"name": name, "from_junction": from_junctions, "to_junction": to_junctions,
                "qext_w": qext_w, "loss_coefficient": loss_coefficient, "in_service": bool(in_service), "type": type}
     _set_multiple_entries(net, "heat_exchanger", index, **entries, **kwargs)
 
     return index
 
 
-def create_heat_consumers(net, from_junctions, to_junctions, diameter_m, qext_w=None, controlled_mdot_kg_per_s=None,
+def create_heat_consumers(net, from_junctions, to_junctions, qext_w=None, controlled_mdot_kg_per_s=None,
                           deltat_k=None, treturn_k=None, name=None, index=None, in_service=True, type="heat_consumer",
                           **kwargs):
     """
@@ -1708,8 +1712,6 @@ def create_heat_consumers(net, from_junctions, to_junctions, diameter_m, qext_w=
     :param to_junctions: IDs of the junctions on the other side which the heat consumers will be \
         connected with
     :type to_junctions: Iterable(int)
-    :param diameter_m: The heat consumers' inner diameter in [m] - only for result calculation
-    :type diameter_m: Iterable(float) or float
     :param qext_w: External heat flux in [W]. If positive, heat is extracted from the network. If \
         negative, heat is being fed into the network from a heat source.
     :type qext_w: Iterable(float) or float, default None
@@ -1736,9 +1738,12 @@ def create_heat_consumers(net, from_junctions, to_junctions, diameter_m, qext_w=
     :rtype: int
 
     :Example:
-        >>> create_heat_consumers(net,from_junctions=[0, 3], to_junctions=[1, 5], diameter_m=40e-3,
+        >>> create_heat_consumers(net,from_junctions=[0, 3], to_junctions=[1, 5],
         >>>                       qext_w=20000, controlled_mdot_kg_per_s=[0.5, 0.9])
     """
+    if 'diameter_m' in kwargs:
+        logger.warning(r'diameter_m is deprecated as it has no effect on the calculation and results. Nonetheless, '
+                       r'it will be stored in the compoent table for postprocessing purposes by you if required.')
     check_vars = [controlled_mdot_kg_per_s, qext_w, deltat_k, treturn_k]
     var_sums = np.zeros([4, len(from_junctions)])
     for i, cv in enumerate(check_vars):
@@ -1756,7 +1761,7 @@ def create_heat_consumers(net, from_junctions, to_junctions, diameter_m, qext_w=
     index = _get_multiple_index_with_check(net, "heat_consumer", index, len(from_junctions))
     _check_branches(net, from_junctions, to_junctions, "heat_consumer")
 
-    entries = {"name": name, "from_junction": from_junctions, "to_junction": to_junctions, "diameter_m": diameter_m,
+    entries = {"name": name, "from_junction": from_junctions, "to_junction": to_junctions,
                "qext_w": qext_w, "controlled_mdot_kg_per_s": controlled_mdot_kg_per_s, "deltat_k": deltat_k,
                "treturn_k": treturn_k, "in_service": bool(in_service), "type": type}
     _set_multiple_entries(net, "heat_consumer", index, **entries, **kwargs)
diff --git a/src/pandapipes/idx_branch.py b/src/pandapipes/idx_branch.py
index 80eb33d8..58cc0f9f 100644
--- a/src/pandapipes/idx_branch.py
+++ b/src/pandapipes/idx_branch.py
@@ -3,40 +3,50 @@
 # Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
 
 # branch types
-# no types defined
+CIRC = 1  # Circ pump branch
+PC = 2  # Pressure controller branch
 
 # branch indices
 TABLE_IDX = 0  # number of the table that this branch belongs to
 ELEMENT_IDX = 1  # index of the element that this branch belongs to (within the given table)
-FROM_NODE = 2  # f, from bus number
-TO_NODE = 3  # t, to bus number
-ACTIVE = 4
-LENGTH = 5  # Pipe length in [m]
-D = 6  # Diameter in [m]
-AREA = 7  # Area in [m²]
-K = 8  # Pipe roughness in [m]
-MDOTINIT = 9  # mass in  [m/s]
-RE = 10 # Reynolds number
-LAMBDA = 11  # Lambda
-JAC_DERIV_DM = 12  # Slot for the derivative by mass
-JAC_DERIV_DP = 13  # Slot for the derivative by pressure from_node
-JAC_DERIV_DP1 = 14  # Slot for the derivative by pressure to_node
-LOAD_VEC_BRANCHES = 15  # Slot for the load vector for the branches
-JAC_DERIV_DM_NODE = 16  # Slot for the derivative by mass for the nodes connected to branch
-LOAD_VEC_NODES = 17  # Slot for the load vector of the nodes connected to branch
-LOSS_COEFFICIENT = 18
-ALPHA = 19  # Slot for heat transfer coefficient
-JAC_DERIV_DT = 20
-JAC_DERIV_DTOUT = 21
-LOAD_VEC_BRANCHES_T = 22
-TOUTINIT = 23  # Internal slot for outlet pipe temperature
-JAC_DERIV_DT_NODE = 24  # Slot for the derivative fpr T for the nodes connected to branch
-LOAD_VEC_NODES_T = 25 # Slot for the load vector of the nodes connected to branch
-FROM_NODE_T_SWITCHED = 26 # flag to indicate if the from and to node are switched in the thermal calculation
-QEXT = 27  # heat input into the branch [W]
-TEXT = 28 # temperature of surrounding [K]
-PL = 29 # Pressure lift [bar]
-TL = 30 # Temperature lift [K]
-BRANCH_TYPE = 31  # branch type relevant for the pressure controller
+BRANCH_TYPE = 2  # branch type relevant for the pressure controller
+PUMP_TYPE = 3
+FROM_NODE = 4  # f, from bus number
+TO_NODE = 5  # t, to bus number
+ACTIVE = 6
+LENGTH = 7  # Pipe length in [m]
+D = 8  # Diameter in [m]
+AREA = 9  # Area in [m²]
+K = 10  # Pipe roughness in [m]
+RE = 11 # Reynolds number
+LAMBDA = 12  # Lambda
+LOSS_COEFFICIENT = 13
+ALPHA = 14  # Slot for heat transfer coefficient
+QEXT = 15  # heat input into the branch [W]
+TEXT = 16 # temperature of surrounding [K]
+PL = 17 # Pressure lift [bar]
+TL = 18 # Temperature lift [K]
 
-branch_cols = 32
+MDOTINIT = 19  # mass in  [m/s]
+MDOTINIT_T = 20
+FROM_NODE_T_SWITCHED = 21 # flag to indicate if the from and to node are switched in the thermal calculation
+TOUTINIT = 22  # Internal slot for outlet pipe temperature
+
+JAC_DERIV_DM = 23  # Slot for the derivative by mass
+JAC_DERIV_DP = 24  # Slot for the derivative by pressure from_node
+JAC_DERIV_DP1 = 25  # Slot for the derivative by pressure to_node
+JAC_DERIV_DM_NODE = 26  # Slot for the derivative by mass for the nodes connected to branch
+LOAD_VEC_BRANCHES = 27  # Slot for the load vector for the branches
+LOAD_VEC_NODES_FROM = 28  # Slot for the load vector of the from nodes connected to branch
+LOAD_VEC_NODES_TO = 29  # Slot for the load vector of the to nodes connected to branch
+
+
+JAC_DERIV_DT = 30
+JAC_DERIV_DTOUT = 31
+JAC_DERIV_DT_NODE = 32  # Slot for the node equation derivative of T for the nodes branch is connected to
+JAC_DERIV_DTOUT_NODE = 33  # Slot for the node equation derivative of T for the corresponding branch
+LOAD_VEC_BRANCHES_T = 34
+LOAD_VEC_NODES_FROM_T = 35 # Slot for the load vector of the from nodes connected to branch
+LOAD_VEC_NODES_TO_T = 36 # Slot for the load vector of the to nodes connected to branch
+
+branch_cols = 37
diff --git a/src/pandapipes/idx_node.py b/src/pandapipes/idx_node.py
index c4394f80..194fbba9 100644
--- a/src/pandapipes/idx_node.py
+++ b/src/pandapipes/idx_node.py
@@ -5,23 +5,28 @@
 # node types
 P = 1  # Reference node, pressure is fixed
 L = 2  # All other nodes
-T = 10  # Reference node with fixed temperature, otherwise 0
-PC = 20  # Controlled node with fixed pressure p
-NONE = 3  # None
+T = 3  # Reference node with fixed temperature, otherwise 0
+PC = 4  # Controlled node with fixed pressure p
 
 # node indices
 TABLE_IDX = 0  # number of the table that this node belongs to
 ELEMENT_IDX = 1  # index of the element that this node belongs to (within the given table)
 NODE_TYPE = 2  # junction type
-ACTIVE = 3
-PINIT = 4
-LOAD = 5
-HEIGHT = 6
-TINIT = 7
-PAMB = 8 # Ambient pressure in [bar]
-LOAD_T = 9  # Heat power drawn in [W]
-NODE_TYPE_T = 10
-EXT_GRID_OCCURENCE = 11
-EXT_GRID_OCCURENCE_T = 12
+NODE_TYPE_T = 3
+ACTIVE = 4
+HEIGHT = 5
+PAMB = 6 # Ambient pressure in [bar]
+LOAD = 7
+LOAD_T = 8  # Heat power drawn in [W]
+EXT_GRID_OCCURENCE = 9
+EXT_GRID_OCCURENCE_T = 10
+INFEED = 11
+VAR_MASS_SLACK = 12 #required as slack do not necesseraly allow mass different from zero
 
-node_cols = 13
+PINIT = 13
+MDOTSLACKINIT = 14
+TINIT = 15
+
+JAC_DERIV_MSL = 16
+
+node_cols = 17
\ No newline at end of file
diff --git a/src/pandapipes/io/file_io.py b/src/pandapipes/io/file_io.py
index 15766b03..40ab3b5b 100644
--- a/src/pandapipes/io/file_io.py
+++ b/src/pandapipes/io/file_io.py
@@ -5,11 +5,10 @@
 import json
 import os
 import pickle
-from functools import partial
 
 from pandapower.io_utils import PPJSONEncoder, to_dict_with_coord_transform, \
     get_raw_data_from_pickle, transform_net_with_df_and_geo, PPJSONDecoder
-from pandapower.io_utils import pp_hook, encrypt_string, decrypt_string
+from pandapower.io_utils import encrypt_string, decrypt_string
 
 from pandapipes.io.convert_format import convert_format
 from pandapipes.io.io_utils import isinstance_partial, FromSerializableRegistryPpipe
diff --git a/src/pandapipes/pf/build_system_matrix.py b/src/pandapipes/pf/build_system_matrix.py
index a6a58cf2..5f0217e1 100644
--- a/src/pandapipes/pf/build_system_matrix.py
+++ b/src/pandapipes/pf/build_system_matrix.py
@@ -5,11 +5,13 @@
 import numpy as np
 from scipy.sparse import csr_matrix
 
-from pandapipes.idx_branch import FROM_NODE, TO_NODE, JAC_DERIV_DM, JAC_DERIV_DP, JAC_DERIV_DP1, \
-    JAC_DERIV_DM_NODE, LOAD_VEC_NODES, LOAD_VEC_BRANCHES, JAC_DERIV_DT, JAC_DERIV_DTOUT, \
-    JAC_DERIV_DT_NODE, LOAD_VEC_NODES_T, LOAD_VEC_BRANCHES_T, BRANCH_TYPE
-from pandapipes.idx_node import LOAD, TINIT
-from pandapipes.idx_node import P, PC, NODE_TYPE, T, NODE_TYPE_T
+from pandapipes.idx_branch import (FROM_NODE, TO_NODE, JAC_DERIV_DM, JAC_DERIV_DP, JAC_DERIV_DP1, \
+    JAC_DERIV_DM_NODE, LOAD_VEC_NODES_FROM, LOAD_VEC_NODES_TO, LOAD_VEC_BRANCHES, JAC_DERIV_DT, JAC_DERIV_DTOUT,
+    CIRC, PC as PC_BRANCH, JAC_DERIV_DT_NODE, JAC_DERIV_DTOUT_NODE, LOAD_VEC_NODES_FROM_T, LOAD_VEC_NODES_TO_T, \
+    LOAD_VEC_BRANCHES_T, BRANCH_TYPE)
+
+from pandapipes.idx_node import (P, PC as PC_NODE, NODE_TYPE, T, NODE_TYPE_T, LOAD, LOAD_T, INFEED,
+                                 MDOTSLACKINIT, JAC_DERIV_MSL)
 from pandapipes.pf.internals_toolbox import _sum_by_group_sorted, _sum_by_group, \
     get_from_nodes_corrected, get_to_nodes_corrected
 from pandapipes.pf.pipeflow_setup import get_net_option
@@ -32,15 +34,16 @@ def build_system_matrix(net, branch_pit, node_pit, heat_mode):
     """
     update_option = get_net_option(net, "only_update_hydraulic_matrix")
     update_only = update_option and "hydraulic_data_sorting" in net["_internal_data"] \
-        and "hydraulic_matrix" in net["_internal_data"]
+                  and "hydraulic_matrix" in net["_internal_data"]
     use_numba = get_net_option(net, "use_numba")
 
     len_b = len(branch_pit)
     len_n = len(node_pit)
     branch_matrix_indices = np.arange(len_b) + len_n
-    ntyp_col, slack_type, pc_type, num_der = \
-        (NODE_TYPE, P, PC, 3) if not heat_mode else (NODE_TYPE_T, T, PC, 2)
-    pc_nodes = np.where(node_pit[:, ntyp_col] == pc_type)[0]
+    ntyp_col, slack_type, pcn_type, pcb_type, branch_type, num_der = \
+        (NODE_TYPE, P, PC_NODE, PC_BRANCH, BRANCH_TYPE, 3) \
+            if not heat_mode else (NODE_TYPE_T, T, None, None, BRANCH_TYPE, 2)
+    pc_nodes = np.where(node_pit[:, ntyp_col] == pcn_type)[0]
 
     if not heat_mode:
         fn = branch_pit[:, FROM_NODE].astype(np.int32)
@@ -48,116 +51,168 @@ def build_system_matrix(net, branch_pit, node_pit, heat_mode):
     else:
         fn = get_from_nodes_corrected(branch_pit)
         tn = get_to_nodes_corrected(branch_pit)
-    not_slack_fn_branch_mask = node_pit[fn, ntyp_col] != slack_type
-    not_slack_tn_branch_mask = node_pit[tn, ntyp_col] != slack_type
-    pc_branch_mask = branch_pit[:, BRANCH_TYPE] == pc_type
+    pc_branch_mask = branch_pit[:, branch_type] == pcb_type
     slack_nodes = np.where(node_pit[:, ntyp_col] == slack_type)[0]
     pc_matrix_indices = branch_matrix_indices[pc_branch_mask]
 
+    # size of the matrix
     if not heat_mode:
+        len_sl = len(slack_nodes)
+        slack_mass_matrix_indices = np.arange(len_sl) + len_b + len_n
+        slack_masses_from, slack_branches_from = np.where(branch_pit[:, FROM_NODE] == slack_nodes[:, None])
+        slack_masses_to, slack_branches_to = np.where(branch_pit[:, TO_NODE] == slack_nodes[:, None])
+        not_slack_fn_branch_mask = node_pit[fn, ntyp_col] != slack_type
+        not_slack_tn_branch_mask = node_pit[tn, ntyp_col] != slack_type
         len_fn_not_slack = np.sum(not_slack_fn_branch_mask)
         len_tn_not_slack = np.sum(not_slack_tn_branch_mask)
         len_fn1 = num_der * len_b + len_fn_not_slack
         len_tn1 = len_fn1 + len_tn_not_slack
         len_pc = len_tn1 + pc_nodes.shape[0]
-        full_len = len_pc + slack_nodes.shape[0]
+        len_slack = len_pc + slack_nodes.shape[0]
+        len_fsb = len_slack + len(slack_branches_from)
+        len_tsb = len_fsb + len(slack_branches_to)
+        full_len = len_tsb + slack_nodes.shape[0]
     else:
-        inc_flow_sum = np.zeros(len(node_pit[:, LOAD]))
-        tn_unique_der, tn_sums_der = _sum_by_group(use_numba, tn, branch_pit[:, JAC_DERIV_DT_NODE])
-        inc_flow_sum[tn_unique_der] += tn_sums_der
-        len_fn1 = num_der * len_b + len(tn_unique_der)
-        len_tn1 = len_fn1 + len_b
-        full_len = len_tn1 + slack_nodes.shape[0]
+        len_sl = 0
+        not_slack_branch_mask = branch_pit[:, branch_type] != CIRC
+        len_tn_not_slack = np.sum(not_slack_branch_mask)
+        infeed_node = np.arange(len_n)[node_pit[:, INFEED].astype(np.bool_)]
+        len_tn1 = num_der * len_b + len_tn_not_slack
+        len_tn2 = len_tn1 + len_tn_not_slack
+        full_len = len_tn2 + slack_nodes.shape[0]
 
     system_data = np.zeros(full_len, dtype=np.float64)
 
+    # entries in the matrix
     if not heat_mode:
-        # pdF_dv
+
+        # branch equations
+        # ----------------
+        # branch_dF_dm
         system_data[:len_b] = branch_pit[:, JAC_DERIV_DM]
-        # pdF_dpi
+        # branch_dF_dp_from
         system_data[len_b:2 * len_b] = branch_pit[:, JAC_DERIV_DP]
-        # pdF_dpi1
+        # branch_dF_dp_to
         system_data[2 * len_b:3 * len_b] = branch_pit[:, JAC_DERIV_DP1]
-        # jdF_dv_from_nodes
+
+        # node equations
+        # --------------
+        # from_node_dF_dm
         system_data[3 * len_b:len_fn1] = branch_pit[not_slack_fn_branch_mask, JAC_DERIV_DM_NODE] * (-1)
-        # jdF_dv_to_nodes
+        # to_node_dF_dm
         system_data[len_fn1:len_tn1] = branch_pit[not_slack_tn_branch_mask, JAC_DERIV_DM_NODE]
-        # pc_nodes and p_nodes
-        system_data[len_tn1:] = 1
+
+        # fixed pressure equations
+        # ------------------------
+        # pc_nodes and slack_nodes
+        system_data[len_tn1:len_slack] = 1
+
+        # mass flow slack equation
+        # --------------
+        # from_slack_dF_dm
+        system_data[len_slack:len_fsb] = branch_pit[slack_branches_from, JAC_DERIV_DM_NODE] * (-1)
+        # to_slack_dF_dm
+        system_data[len_fsb:len_tsb] = branch_pit[slack_branches_to, JAC_DERIV_DM_NODE]
+        # slackmass_dF_dmslack
+        system_data[len_tsb:] = node_pit[slack_nodes, JAC_DERIV_MSL]
     else:
+
+        # branch equations
+        # ----------------
+        # branch_dF_dT_from
         system_data[:len_b] = branch_pit[:, JAC_DERIV_DT]
-        # pdF_dpi1
+        # branch_dF_dT_out
         system_data[len_b:2 * len_b] = branch_pit[:, JAC_DERIV_DTOUT]
-        # jdF_dv_from_nodes
-        system_data[2 * len_b:len_fn1] = inc_flow_sum[tn_unique_der]
-        # jdF_dv_to_nodes
-        data = branch_pit[:, JAC_DERIV_DT_NODE] * (-1)
-        rows = tn
-        index = np.argsort(rows)
-        data = data[index]
-        system_data[len_fn1:len_fn1 + len_b] = data
-        system_data[len_fn1 + len_b:] = 1
 
+        # node equations
+        # --------------
+        # node_dF_dT_to
+        system_data[2 * len_b:len_tn1] = branch_pit[not_slack_branch_mask, JAC_DERIV_DT_NODE]
+        # node_dF_dT_out
+        system_data[len_tn1:len_tn2] = branch_pit[not_slack_branch_mask, JAC_DERIV_DTOUT_NODE]
+
+        # fixed temperature equations
+        # ---------------------------
+        # t_nodes
+        system_data[len_tn2:] = 1
+
+    # position in the matrix
     if not update_only:
         system_cols = np.zeros(full_len, dtype=np.int32)
         system_rows = np.zeros(full_len, dtype=np.int32)
 
         if not heat_mode:
-            # pdF_dv
+
+            # branch equations
+            # ----------------
+            # branch_dF_dm
             system_cols[:len_b] = branch_matrix_indices
             system_rows[:len_b] = branch_matrix_indices
-
-            # pdF_dpi
+            # branch_dF_dp_from
             system_cols[len_b:2 * len_b] = fn
             system_rows[len_b:2 * len_b] = branch_matrix_indices
-
-            # pdF_dpi1
+            # branch_dF_dp_to
             system_cols[2 * len_b:3 * len_b] = tn
             system_rows[2 * len_b:3 * len_b] = branch_matrix_indices
 
-            # jdF_dv_from_nodes
+            # node equations
+            # --------------
+            # from_node_dF_dm
             system_cols[3 * len_b:len_fn1] = branch_matrix_indices[not_slack_fn_branch_mask]
             system_rows[3 * len_b:len_fn1] = fn[not_slack_fn_branch_mask]
-
-            # jdF_dv_to_nodes
+            # to_node_dF_dm
             system_cols[len_fn1:len_tn1] = branch_matrix_indices[not_slack_tn_branch_mask]
             system_rows[len_fn1:len_tn1] = tn[not_slack_tn_branch_mask]
 
+            # fixed pressure equations
+            # -----------------------
             # pc_nodes
             system_cols[len_tn1:len_pc] = pc_nodes
             system_rows[len_tn1:len_pc] = pc_matrix_indices
+            # slack_nodes
+            system_cols[len_pc:len_slack] = slack_nodes
+            system_rows[len_pc:len_slack] = slack_nodes
+
+            # mass flow slack equation
+            # --------------
+            # from_slack_dF_dm
+            system_cols[len_slack:len_fsb] = branch_matrix_indices[slack_branches_from]
+            system_rows[len_slack:len_fsb] = slack_mass_matrix_indices[slack_masses_from]
+            # to_slack_dF_dm
+            system_cols[len_fsb:len_tsb] = branch_matrix_indices[slack_branches_to]
+            system_rows[len_fsb:len_tsb] = slack_mass_matrix_indices[slack_masses_to]
+            # slackmass_dF_dmslack
+            system_cols[len_tsb:] = slack_mass_matrix_indices
+            system_rows[len_tsb:] = slack_mass_matrix_indices
 
-            # p_nodes
-            system_cols[len_pc:] = slack_nodes
-            system_rows[len_pc:] = slack_nodes
         else:
-            # pdF_dTfromnode
+            # branch equations
+            # ----------------
+            # branch_dF_dT_from
             system_cols[:len_b] = fn
             system_rows[:len_b] = branch_matrix_indices
-
-            # pdF_dTout
+            # branch_dF_dT_out
             system_cols[len_b:2 * len_b] = branch_matrix_indices
             system_rows[len_b:2 * len_b] = branch_matrix_indices
 
-            # t_nodes
-            system_cols[len_fn1 + len_b:] = slack_nodes
-            system_rows[len_fn1 + len_b:] = np.arange(0, len(slack_nodes))
-
-            # jdF_dTnode_
-            tn_unique_idx = np.unique(tn, return_index=True)
-            system_cols[2 * len_b:len_fn1] = tn_unique_idx[0]
-            system_rows[2 * len_b:len_fn1] = len(slack_nodes) + np.arange(0, len(tn_unique_der))
+            # node equations
+            # --------------
+            # node_dF_dT_to
+            system_cols[2 * len_b:len_tn1] = tn[not_slack_branch_mask]
+            system_rows[2 * len_b:len_tn1] = tn[not_slack_branch_mask]
+            # node_dF_dT_out
+            system_cols[len_tn1:len_tn2] = branch_matrix_indices[not_slack_branch_mask]
+            system_rows[len_tn1:len_tn2] = tn[not_slack_branch_mask]
 
-            # jdF_dTout
-            branch_order = np.argsort(tn)
-            tn_uni, tn_uni_counts = np.unique(tn[branch_order], return_counts=True)
-            row_index = np.repeat(np.arange(len(tn_uni)), tn_uni_counts)
-            system_cols[len_fn1:len_fn1 + len_b] = branch_matrix_indices[branch_order]
-            system_rows[len_fn1:len_fn1 + len_b] = len(slack_nodes) + row_index
+            # fixed temperature equations
+            # ---------------------------
+            # t_nodes (overwrites only infeeding nodes' equation)
+            system_cols[len_tn2:] = slack_nodes
+            system_rows[len_tn2:] = infeed_node
 
         if not update_option:
             system_matrix = csr_matrix((system_data, (system_rows, system_cols)),
-                                       shape=(len_n + len_b, len_n + len_b))
+                                       shape=(len_n + len_b + len_sl, len_n + len_b + len_sl))
 
         else:
             data_order = np.lexsort([system_cols, system_rows])
@@ -165,12 +220,12 @@ def build_system_matrix(net, branch_pit, node_pit, heat_mode):
             system_cols = system_cols[data_order]
             system_rows = system_rows[data_order]
 
-            row_counter = np.zeros(len_b + len_n + 1, dtype=np.int32)
+            row_counter = np.zeros(len_b + len_n + len_sl + 1, dtype=np.int32)
             unique_rows, row_counts = _sum_by_group_sorted(system_rows, np.ones_like(system_rows))
             row_counter[unique_rows + 1] += row_counts
             ptr = row_counter.cumsum()
             system_matrix = csr_matrix((system_data, system_cols, ptr),
-                                       shape=(len_n + len_b, len_n + len_b))
+                                       shape=(len_n + len_b + len_sl, len_n + len_b + len_sl))
             net["_internal_data"]["hydraulic_data_sorting"] = data_order
             net["_internal_data"]["hydraulic_matrix"] = system_matrix
     else:
@@ -179,24 +234,36 @@ def build_system_matrix(net, branch_pit, node_pit, heat_mode):
         system_matrix = net["_internal_data"]["hydraulic_matrix"]
         system_matrix.data = system_data
 
+    # load vector on the right side
     if not heat_mode:
-        load_vector = np.empty(len_n + len_b)
-        load_vector[len_n:] = branch_pit[:, LOAD_VEC_BRANCHES]
+        load_vector = np.empty(len_n + len_b + len_sl)
+        load_vector[len_n:len_b + len_n] = branch_pit[:, LOAD_VEC_BRANCHES]
         load_vector[:len_n] = node_pit[:, LOAD] * (-1)
-        fn_unique, fn_sums = _sum_by_group(use_numba, fn, branch_pit[:, LOAD_VEC_NODES])
-        tn_unique, tn_sums = _sum_by_group(use_numba, tn, branch_pit[:, LOAD_VEC_NODES])
+        fn_unique, fn_sums = _sum_by_group(use_numba, fn, branch_pit[:, LOAD_VEC_NODES_FROM])
+        tn_unique, tn_sums = _sum_by_group(use_numba, tn, branch_pit[:, LOAD_VEC_NODES_TO])
         load_vector[fn_unique] -= fn_sums
         load_vector[tn_unique] += tn_sums
         load_vector[slack_nodes] = 0
         load_vector[pc_matrix_indices] = 0
+
+        load_vector[slack_mass_matrix_indices] = node_pit[slack_nodes, LOAD] * (-1)
+        fsb_unique, fsb_sums = _sum_by_group(use_numba, slack_masses_from,
+                                             branch_pit[slack_branches_from, LOAD_VEC_NODES_FROM])
+        tsb_unique, tsb_sums = _sum_by_group(use_numba, slack_masses_to,
+                                             branch_pit[slack_branches_to, LOAD_VEC_NODES_TO])
+        load_vector[slack_mass_matrix_indices[fsb_unique]] -= fsb_sums
+        load_vector[slack_mass_matrix_indices[tsb_unique]] += tsb_sums
+        load_vector[slack_mass_matrix_indices] -= node_pit[slack_nodes, MDOTSLACKINIT]
     else:
-        tn_unique, tn_sums = _sum_by_group(use_numba, tn, branch_pit[:, LOAD_VEC_NODES_T])
         load_vector = np.zeros(len_n + len_b)
-        load_vector[len(slack_nodes) + np.arange(0, len(tn_unique_der))] += tn_sums
-        load_vector[len(slack_nodes) + np.arange(0, len(tn_unique_der))] -= tn_sums_der * node_pit[
-            tn_unique_der, TINIT]
-        load_vector[0:len(slack_nodes)] = 0.
-
         load_vector[len_n:] = branch_pit[:, LOAD_VEC_BRANCHES_T]
+        load_vector[:len_n] = node_pit[:, LOAD_T] * (-1)
+        # This approach can be used if you consider the effect of sources with given temperature
+        # fn_unique, fn_sums = _sum_by_group(use_numba, fn, branch_pit[:, LOAD_VEC_NODES_FROM_T])
+        fn_unique, fn_sums = _sum_by_group(use_numba, tn, branch_pit[:, LOAD_VEC_NODES_FROM_T])
+        tn_unique, tn_sums = _sum_by_group(use_numba, tn, branch_pit[:, LOAD_VEC_NODES_TO_T])
+        load_vector[fn_unique] -= fn_sums
+        load_vector[tn_unique] += tn_sums
+        load_vector[infeed_node] = 0
 
     return system_matrix, load_vector
diff --git a/src/pandapipes/pf/derivative_calculation.py b/src/pandapipes/pf/derivative_calculation.py
index 8a47fa0e..4b1fe058 100644
--- a/src/pandapipes/pf/derivative_calculation.py
+++ b/src/pandapipes/pf/derivative_calculation.py
@@ -2,11 +2,12 @@
 
 from pandapipes.constants import NORMAL_TEMPERATURE
 from pandapipes.idx_branch import LENGTH, D, K, RE, LAMBDA, LOAD_VEC_BRANCHES, \
-    JAC_DERIV_DM, JAC_DERIV_DP, JAC_DERIV_DP1, LOAD_VEC_NODES, JAC_DERIV_DM_NODE, \
-    FROM_NODE, TO_NODE, TOUTINIT, TEXT, AREA, ALPHA, TL, QEXT, LOAD_VEC_NODES_T, \
-    LOAD_VEC_BRANCHES_T, JAC_DERIV_DT, JAC_DERIV_DTOUT, JAC_DERIV_DT_NODE, MDOTINIT
-from pandapipes.idx_node import TINIT as TINIT_NODE
-from pandapipes.pf.internals_toolbox import get_from_nodes_corrected
+    JAC_DERIV_DM, JAC_DERIV_DP, JAC_DERIV_DP1, LOAD_VEC_NODES_FROM, LOAD_VEC_NODES_TO, JAC_DERIV_DM_NODE, \
+    FROM_NODE, TO_NODE, TOUTINIT, TEXT, AREA, ALPHA, TL, QEXT, LOAD_VEC_NODES_FROM_T, LOAD_VEC_NODES_TO_T,\
+    LOAD_VEC_BRANCHES_T, JAC_DERIV_DT, JAC_DERIV_DTOUT, JAC_DERIV_DTOUT_NODE, \
+    JAC_DERIV_DT_NODE, MDOTINIT, BRANCH_TYPE, CIRC
+from pandapipes.idx_node import TINIT as TINIT_NODE, INFEED
+from pandapipes.pf.internals_toolbox import get_from_nodes_corrected, get_to_nodes_corrected
 from pandapipes.properties.fluids import get_fluid
 from pandapipes.properties.properties_toolbox import get_branch_real_density, get_branch_real_eta, \
     get_branch_cp
@@ -52,8 +53,9 @@ def calculate_derivatives_hydraulic(net, branch_pit, node_pit, options):
             from pandapipes.pf.derivative_toolbox import derivatives_hydraulic_incomp_np \
                 as derivatives_hydraulic_incomp
 
-        load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1 = derivatives_hydraulic_incomp(
-            branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs, height_difference, rho)
+        load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1 = (
+            derivatives_hydraulic_incomp(
+            branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs, height_difference, rho))
     else:
         if options["use_numba"]:
             from pandapipes.pf.derivative_toolbox_numba import derivatives_hydraulic_comp_numba \
@@ -69,40 +71,64 @@ def calculate_derivatives_hydraulic(net, branch_pit, node_pit, options):
         # TODO: this might not be required
         der_comp = fluid.get_der_compressibility() * der_p_m
         der_comp1 = fluid.get_der_compressibility() * der_p_m1
-        load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1 = derivatives_hydraulic_comp(
+        load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1 = (
+            derivatives_hydraulic_comp(
             node_pit, branch_pit, lambda_, der_lambda, p_init_i_abs, p_init_i1_abs, height_difference,
-            comp_fact, der_comp, der_comp1, rho, rho_n)
+            comp_fact, der_comp, der_comp1, rho, rho_n))
 
     branch_pit[:, LOAD_VEC_BRANCHES] = load_vec
     branch_pit[:, JAC_DERIV_DM] = df_dm
     branch_pit[:, JAC_DERIV_DP] = df_dp
     branch_pit[:, JAC_DERIV_DP1] = df_dp1
-    branch_pit[:, LOAD_VEC_NODES] = load_vec_nodes
+    branch_pit[:, LOAD_VEC_NODES_FROM] = load_vec_nodes_from
+    branch_pit[:, LOAD_VEC_NODES_TO] = load_vec_nodes_to
     branch_pit[:, JAC_DERIV_DM_NODE] = df_dm_nodes
 
 
 def calculate_derivatives_thermal(net, branch_pit, node_pit, _):
     fluid = get_fluid(net)
     cp = get_branch_cp(fluid, node_pit, branch_pit)
-    m_init = np.abs(branch_pit[:, MDOTINIT])
+    m_init_i = np.abs(branch_pit[:, MDOTINIT])
+    m_init_i1 = np.abs(branch_pit[:, MDOTINIT])
     from_nodes = get_from_nodes_corrected(branch_pit)
+    to_nodes = get_to_nodes_corrected(branch_pit)
     t_init_i = node_pit[from_nodes, TINIT_NODE]
     t_init_i1 = branch_pit[:, TOUTINIT]
+    t_init_n = node_pit[to_nodes, TINIT_NODE]
+    cp_n = fluid.get_heat_capacity(t_init_n)
+    cp_i1 = fluid.get_heat_capacity(t_init_i1)
     t_amb = branch_pit[:, TEXT]
     length = branch_pit[:, LENGTH]
     alpha = branch_pit[:, ALPHA] * np.pi * branch_pit[:, D]
     tl = branch_pit[:, TL]
     qext = branch_pit[:, QEXT]
     t_m = (t_init_i1 + t_init_i) / 2
+    m_m = (m_init_i + m_init_i1) / 2
+    no_cp = branch_pit[:, BRANCH_TYPE] != CIRC
+    infeed_node = np.setdiff1d(from_nodes[no_cp], to_nodes[no_cp])
 
-    branch_pit[:, LOAD_VEC_BRANCHES_T] = \
-        -(cp * m_init * (-t_init_i + t_init_i1 - tl) - alpha * (t_amb - t_m) * length + qext)
+    branch_pit[:, JAC_DERIV_DT] = - cp * m_m + alpha / 2 * length
+    branch_pit[:, JAC_DERIV_DTOUT] = cp * m_m + alpha / 2 * length
+    branch_pit[:, LOAD_VEC_BRANCHES_T] = cp * m_m * (-t_init_i + t_init_i1 - tl) - alpha * (
+                t_amb - t_m) * length + qext
 
-    branch_pit[:, JAC_DERIV_DT] = - cp * m_init + alpha / 2 * length
-    branch_pit[:, JAC_DERIV_DTOUT] = cp * m_init + alpha / 2 * length
+    branch_pit[:, JAC_DERIV_DT_NODE] = - m_init_i * cp_n
+    branch_pit[:, JAC_DERIV_DTOUT_NODE] = m_init_i1 * cp_i1
+    branch_pit[:, LOAD_VEC_NODES_FROM_T] = m_init_i1 * t_init_n * cp_n
+    branch_pit[:, LOAD_VEC_NODES_TO_T] = m_init_i1 * t_init_i1 * cp_i1
 
-    branch_pit[:, JAC_DERIV_DT_NODE] = m_init
-    branch_pit[:, LOAD_VEC_NODES_T] = m_init * t_init_i1
+    # This approach can be used if you consider the effect of sources with given temperature (checkout issue #656)
+
+    # branch_pit[:, LOAD_VEC_NODES_FROM_T] = m_init_i * t_init_i * cp_i
+    # --> cp_i is calculated by fluid.get_heat_capacity(t_init_i)
+    # branch_pit[:, LOAD_VEC_NODES_TO_T] = m_init_i1 * t_init_i1 * cp_i1
+    # --> still missing is the derivative of loads
+    # t_init = node_pit[:, TINIT_NODE]
+    # cp_n = fluid.get_heat_capacity(t_init)
+    # node_pit[:, LOAD_T] = cp_n * node_pit[:, LOAD] * t_init
+
+    node_pit[:, INFEED] = False
+    node_pit[infeed_node, INFEED] = True
 
 
 def get_derived_values(node_pit, from_nodes, to_nodes, use_numba):
diff --git a/src/pandapipes/pf/derivative_toolbox.py b/src/pandapipes/pf/derivative_toolbox.py
index 2b765663..d627560f 100644
--- a/src/pandapipes/pf/derivative_toolbox.py
+++ b/src/pandapipes/pf/derivative_toolbox.py
@@ -34,9 +34,10 @@ def derivatives_hydraulic_incomp_np(branch_pit, der_lambda, p_init_i_abs, p_init
 
     df_dm_nodes = np.ones_like(der_lambda)
     
-    load_vec_nodes = branch_pit[:, MDOTINIT]
+    load_vec_nodes_from = branch_pit[:, MDOTINIT]
+    load_vec_nodes_to = branch_pit[:, MDOTINIT]
 
-    return load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1
+    return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1
 
 
 def derivatives_hydraulic_comp_np(node_pit, branch_pit, lambda_, der_lambda, p_init_i_abs, p_init_i1_abs,
@@ -65,8 +66,11 @@ def derivatives_hydraulic_comp_np(node_pit, branch_pit, lambda_, der_lambda, p_i
                - normal_term * comp_fact * m_init2 * friction_term * p_sum_div * tm
 
     df_dm_nodes = np.ones_like(lambda_)
-    load_vec_nodes = branch_pit[:, MDOTINIT]
-    return load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1
+
+    load_vec_nodes_from = branch_pit[:, MDOTINIT]
+    load_vec_nodes_to = branch_pit[:, MDOTINIT]
+
+    return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1
 
 
 def calc_lambda_nikuradse_incomp_np(m, d, k, eta, area):
diff --git a/src/pandapipes/pf/derivative_toolbox_numba.py b/src/pandapipes/pf/derivative_toolbox_numba.py
index 2893583f..69be5c7f 100644
--- a/src/pandapipes/pf/derivative_toolbox_numba.py
+++ b/src/pandapipes/pf/derivative_toolbox_numba.py
@@ -23,7 +23,8 @@ def derivatives_hydraulic_incomp_numba(branch_pit, der_lambda, p_init_i_abs, p_i
     df_dm = np.zeros_like(der_lambda)
     df_dp = np.ones_like(der_lambda)
     df_dp1 = np.ones_like(der_lambda) * (-1)
-    load_vec_nodes = np.zeros_like(der_lambda)
+    load_vec_nodes_from = np.zeros_like(der_lambda)
+    load_vec_nodes_to = np.zeros_like(der_lambda)
     df_dm_nodes = np.ones_like(der_lambda)
 
     for i in range(le):
@@ -40,8 +41,9 @@ def derivatives_hydraulic_incomp_numba(branch_pit, der_lambda, p_init_i_abs, p_i
 
         load_vec[i] = p_diff + branch_pit[i][PL] + const_height - const_term * m_init2 * friction_term
 
-        load_vec_nodes[i] = branch_pit[i][MDOTINIT]
-    return load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1
+        load_vec_nodes_from[i] = branch_pit[i][MDOTINIT]
+        load_vec_nodes_to[i] = branch_pit[i][MDOTINIT]
+    return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1
 
 
 @jit((float64[:, :], float64[:, :], float64[:], float64[:], float64[:], float64[:], float64[:], float64[:],
@@ -53,10 +55,10 @@ def derivatives_hydraulic_comp_numba(node_pit, branch_pit, lambda_, der_lambda,
     df_dm = np.zeros_like(lambda_)
     df_dp = np.zeros_like(lambda_)
     df_dp1 = np.zeros_like(lambda_)
-    load_vec_nodes = np.zeros_like(der_lambda)
+    load_vec_nodes_from = np.zeros_like(der_lambda)
+    load_vec_nodes_to = np.zeros_like(der_lambda)
     df_dm_nodes = np.ones_like(der_lambda)
     from_nodes = branch_pit[:, FROM_NODE].astype(np.int32)
-    to_nodes = branch_pit[:, TO_NODE].astype(np.int32)
 
     # Formulas for gas pressure loss according to laminar version
     for i in range(le):
@@ -85,8 +87,9 @@ def derivatives_hydraulic_comp_numba(node_pit, branch_pit, lambda_, der_lambda,
         df_dm[i] = -1. * normal_term * comp_fact[i] * p_sum_div * tm * (2 * m_init_abs * friction_term \
             + np.divide(der_lambda[i] * branch_pit[i][LENGTH] * m_init2, branch_pit[i][D]))
 
-        load_vec_nodes[i] = branch_pit[i][MDOTINIT]
-    return load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1
+        load_vec_nodes_from[i] = branch_pit[i][MDOTINIT]
+        load_vec_nodes_to[i] = branch_pit[i][MDOTINIT]
+    return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1
 
 
 @jit((float64[:], float64[:], float64[:], float64[:], float64[:]), nopython=True)
diff --git a/src/pandapipes/pf/pipeflow_setup.py b/src/pandapipes/pf/pipeflow_setup.py
index 2e3466f7..cb4e3fe6 100644
--- a/src/pandapipes/pf/pipeflow_setup.py
+++ b/src/pandapipes/pf/pipeflow_setup.py
@@ -12,7 +12,7 @@
 from pandapipes.idx_branch import FROM_NODE, TO_NODE, branch_cols, MDOTINIT, \
     ACTIVE as ACTIVE_BR
 from pandapipes.idx_node import NODE_TYPE, P, NODE_TYPE_T, node_cols, T, ACTIVE as ACTIVE_ND, \
-    TABLE_IDX as TABLE_IDX_ND, ELEMENT_IDX as ELEMENT_IDX_ND
+    TABLE_IDX as TABLE_IDX_ND, ELEMENT_IDX as ELEMENT_IDX_ND, INFEED
 from pandapipes.pf.internals_toolbox import _sum_by_group
 from pandapipes.properties.fluids import get_fluid
 
@@ -773,6 +773,13 @@ def reduce_pit(net, mode="hydraulics"):
         net["_lookups"]["%s_from_to_active_%s" % (el, mode)] = from_to_active_lookup
 
 
+def check_infeed_number(node_pit):
+    slack_nodes = node_pit[:, NODE_TYPE_T] == T
+    infeed_nodes = node_pit[:, INFEED]
+    if sum(infeed_nodes) != sum(slack_nodes):
+        raise PipeflowNotConverged(r'The number of infeeding nodes and slacks do not match')
+
+
 class PipeflowNotConverged(ppException):
     """
     Exception being raised in case pipeflow did not converge.
diff --git a/src/pandapipes/pipeflow.py b/src/pandapipes/pipeflow.py
index da0cd4d4..ebf0a033 100644
--- a/src/pandapipes/pipeflow.py
+++ b/src/pandapipes/pipeflow.py
@@ -7,14 +7,14 @@
 from scipy.sparse.linalg import spsolve
 
 from pandapipes.idx_branch import MDOTINIT, TOUTINIT, FROM_NODE_T_SWITCHED
-from pandapipes.idx_node import PINIT, TINIT
+from pandapipes.idx_node import PINIT, TINIT, MDOTSLACKINIT, NODE_TYPE, P
 from pandapipes.pf.build_system_matrix import build_system_matrix
 from pandapipes.pf.derivative_calculation import (calculate_derivatives_hydraulic,
                                                   calculate_derivatives_thermal)
 from pandapipes.pf.pipeflow_setup import (
     get_net_option, get_net_options, set_net_option, init_options, create_internal_results,
     write_internal_results, get_lookup, create_lookups, initialize_pit, reduce_pit,
-    set_user_pf_options, init_all_result_tables, identify_active_nodes_branches,
+    set_user_pf_options, init_all_result_tables, identify_active_nodes_branches, check_infeed_number,
     PipeflowNotConverged
 )
 from pandapipes.pf.result_extraction import extract_all_results, extract_results_active_pit
@@ -172,10 +172,10 @@ def hydraulics(net):
     reduce_pit(net, mode="hydraulics")
     if not get_net_option(net, "reuse_internal_data") or "_internal_data" not in net:
         net["_internal_data"] = dict()
-    solver_vars = ['mdot', 'p']
-    tol_p, tol_m = get_net_options(net, 'tol_m', 'tol_p')
-    newton_raphson(net, solve_hydraulics, 'hydraulics', solver_vars, [tol_m, tol_p],
-                   ['branch', 'node'], 'max_iter_hyd')
+    solver_vars = ['mdot', 'p', 'mdotslack']
+    tol_p, tol_m, tol_msl = get_net_options(net, 'tol_m', 'tol_p', 'tol_m')
+    newton_raphson(net, solve_hydraulics, 'hydraulics', solver_vars, [tol_m, tol_p, tol_msl],
+                   ['branch', 'node', 'node'], 'max_iter_hyd')
     if net.converged:
         set_user_pf_options(net, hyd_flag=True)
 
@@ -246,13 +246,17 @@ def solve_hydraulics(net):
 
     m_init_old = branch_pit[:, MDOTINIT].copy()
     p_init_old = node_pit[:, PINIT].copy()
+    slack_nodes = np.where(node_pit[:, NODE_TYPE] == P)[0]
+    msl_init_old = node_pit[slack_nodes, MDOTSLACKINIT].copy()
 
     x = spsolve(jacobian, epsilon)
 
-    branch_pit[:, MDOTINIT] -= x[len(node_pit):]
+    branch_pit[:, MDOTINIT] -= x[len(node_pit):len(node_pit) + len(branch_pit)]
     node_pit[:, PINIT] -= x[:len(node_pit)] * options["alpha"]
+    node_pit[slack_nodes, MDOTSLACKINIT] -= x[len(node_pit) + len(branch_pit):]
 
-    return [branch_pit[:, MDOTINIT], m_init_old, node_pit[:, PINIT], p_init_old], epsilon
+    return [branch_pit[:, MDOTINIT], m_init_old, node_pit[:, PINIT], p_init_old, msl_init_old,
+            node_pit[slack_nodes, MDOTSLACKINIT]], epsilon
 
 
 def solve_temperature(net):
@@ -282,6 +286,8 @@ def solve_temperature(net):
     calculate_derivatives_thermal(net, branch_pit, node_pit, options)
     for comp in net['component_list']:
         comp.adaption_after_derivatives_thermal(net, branch_pit, node_pit, branch_lookups, options)
+    check_infeed_number(node_pit)
+
     jacobian, epsilon = build_system_matrix(net, branch_pit, node_pit, True)
 
     t_init_old = node_pit[:, TINIT].copy()
@@ -289,8 +295,8 @@ def solve_temperature(net):
 
     x = spsolve(jacobian, epsilon)
 
-    node_pit[:, TINIT] += x[:len(node_pit)] * options["alpha"]
-    branch_pit[:, TOUTINIT] += x[len(node_pit):]
+    node_pit[:, TINIT] -= x[:len(node_pit)] * options["alpha"]
+    branch_pit[:, TOUTINIT] -= x[len(node_pit):]
 
     return [branch_pit[:, TOUTINIT], t_out_old, node_pit[:, TINIT], t_init_old], epsilon
 
@@ -337,9 +343,9 @@ def finalize_iteration(net, niter, residual_norm, nonlinear_method, errors, tols
             return
     elif nonlinear_method != "constant":
         logger.warning("No proper nonlinear method chosen. Using constant settings.")
-    converged = True
     for error, var, tol in zip(errors.values(), solver_vars, tols):
-        converged = converged and error[niter] <= tol
+        converged = error[niter] <= tol
+        if not converged: break
         logger.debug("error_%s: %s" % (var, error[niter]))
     net.converged = converged and residual_norm <= tol_res
 
diff --git a/src/pandapipes/plotting/generic_geodata.py b/src/pandapipes/plotting/generic_geodata.py
index 51bc2cd2..9547d09d 100644
--- a/src/pandapipes/plotting/generic_geodata.py
+++ b/src/pandapipes/plotting/generic_geodata.py
@@ -104,7 +104,7 @@ def create_generic_coordinates(net, mg=None, library="igraph", geodata_table="ju
         raise ValueError("Unknown library %s - chose 'igraph' or 'networkx'" % library)
 
     if len(coords):
-        net[geodata_table].x = coords[1]
-        net[geodata_table].y = coords[0]
+        net[geodata_table].loc[:, 'x'] = coords[1]
+        net[geodata_table].loc[:, 'y'] = coords[0]
         net[geodata_table].index = net.junction.index if junctions is None else junctions
     return net
diff --git a/src/pandapipes/test/api/release_cycle/release_control_test_network.py b/src/pandapipes/test/api/release_cycle/release_control_test_network.py
index c8540b28..a3954c87 100644
--- a/src/pandapipes/test/api/release_cycle/release_control_test_network.py
+++ b/src/pandapipes/test/api/release_cycle/release_control_test_network.py
@@ -73,6 +73,9 @@ def release_control_test_network(max_iter_hyd, save=False):
     pp.create_ext_grid(net, junction=0, p_bar=3, t_k=300, name="External Grid 0", in_service=True,
                        index=None, type="pt")
     pp.create_ext_grid(net, 1, p_bar=5, t_k=350, name="External Grid 1", type="pt")
+    pp.create_ext_grid(net, junction=4, p_bar=6, t_k=290, name="External Grid 2", in_service=True,
+                       index=None, type="pt")
+    pp.create_ext_grid(net, 5, p_bar=5, t_k=290, name="External Grid 3", type="pt")
 
     # sinks
     pp.create_sink(net, junction=2, mdot_kg_per_s=0.2, scaling=1., name="Sink 0", index=None,
@@ -189,6 +192,9 @@ def release_control_test_network_water(max_iter_hyd, save=False):
     pp.create_ext_grid(net, junction=0, p_bar=3, t_k=300, name="External Grid 0", in_service=True,
                        index=None, type="pt")
     pp.create_ext_grid(net, 1, p_bar=5, t_k=350, name="External Grid 1", type="pt")
+    pp.create_ext_grid(net, junction=4, p_bar=6, t_k=290, name="External Grid 2", in_service=True,
+                       index=None, type="pt")
+    pp.create_ext_grid(net, 5, p_bar=5, t_k=290, name="External Grid 3", type="pt")
 
     # sinks
     pp.create_sink(net, junction=2, mdot_kg_per_s=0.2, scaling=1., name="Sink 0", index=None,
diff --git a/src/pandapipes/test/api/test_components/test_circ_pump_mass.py b/src/pandapipes/test/api/test_components/test_circ_pump_mass.py
index 59988e2d..158cf54d 100644
--- a/src/pandapipes/test/api/test_components/test_circ_pump_mass.py
+++ b/src/pandapipes/test/api/test_components/test_circ_pump_mass.py
@@ -28,7 +28,7 @@ def test_circulation_pump_constant_mass(use_numba):
     pandapipes.create_pipe_from_parameters(net, j1, j2, k_mm=1., length_km=0.43380, diameter_m=0.1022)
     pandapipes.create_pipe_from_parameters(net, j3, j4, k_mm=1., length_km=0.26370, diameter_m=0.1022)
     pandapipes.create_circ_pump_const_mass_flow(net, j4, j1, 5, 5, 300, type='pt')
-    pandapipes.create_heat_exchanger(net, j2, j3, 0.1, qext_w=200000)
+    pandapipes.create_heat_exchanger(net, j2, j3, qext_w=200000)
     pandapipes.create_sink(net, j1, 2)
     pandapipes.create_source(net, j4, 2)
 
@@ -50,8 +50,9 @@ def test_circulation_pump_constant_mass(use_numba):
     p_diff = np.abs(1 - res_junction.p_bar.values / data['p'].dropna().values)
     t_diff = np.abs(1 - res_junction.t_k.values / data['t'].dropna().values)
     v_diff = np.abs(1 - res_pipe / data['v'].dropna().values)
-    mdot_diff = np.abs(1 - res_pump['mdot_flow_kg_per_s'].values / data['mdot'].dropna().values)
-    deltap_diff = np.abs(1 - res_pump['deltap_bar'].values / data['deltap'].dropna().values)
+    mdot_diff = np.abs(1 - res_pump['mdot_from_kg_per_s'].values / data['mdot'].dropna().values)
+    deltap_diff = np.abs(
+        1 - (res_pump['p_to_bar'].values - res_pump['p_from_bar'].values) / data['deltap'].dropna().values)
 
     assert np.all(p_diff < 0.01)
     assert np.all(t_diff < 0.01)
diff --git a/src/pandapipes/test/api/test_components/test_circ_pump_pressure.py b/src/pandapipes/test/api/test_components/test_circ_pump_pressure.py
index 925eea90..95fcc505 100644
--- a/src/pandapipes/test/api/test_components/test_circ_pump_pressure.py
+++ b/src/pandapipes/test/api/test_components/test_circ_pump_pressure.py
@@ -29,7 +29,7 @@ def test_circulation_pump_constant_pressure(use_numba):
     pandapipes.create_pipe_from_parameters(net, j3, j4, k_mm=1., length_km=0.26370,
                                            diameter_m=0.1022)
     pandapipes.create_circ_pump_const_pressure(net, j4, j1, 5, 2, 300, type='pt')
-    pandapipes.create_heat_exchanger(net, j2, j3, 0.1, qext_w=200000)
+    pandapipes.create_heat_exchanger(net, j2, j3, qext_w=200000)
     pandapipes.create_sink(net, j1, 2)
     pandapipes.create_source(net, j4, 2)
 
@@ -51,8 +51,9 @@ def test_circulation_pump_constant_pressure(use_numba):
     p_diff = np.abs(1 - res_junction.p_bar.values / data['p'].dropna().values)
     t_diff = np.abs(1 - res_junction.t_k.values / data['t'].dropna().values)
     v_diff = np.abs(1 - res_pipe / data['v'].dropna().values)
-    mdot_diff = np.abs(1 - res_pump['mdot_flow_kg_per_s'].values / data['mdot'].dropna().values)
-    deltap_diff = np.abs(1 - res_pump['deltap_bar'].values / data['deltap'].dropna().values)
+    mdot_diff = np.abs(1 - res_pump['mdot_from_kg_per_s'].values / data['mdot'].dropna().values)
+    deltap_diff = np.abs(
+        1 - (res_pump['p_to_bar'].values - res_pump['p_from_bar'].values) / data['deltap'].dropna().values)
 
     assert np.all(p_diff < 0.01)
     assert np.all(t_diff < 0.01)
diff --git a/src/pandapipes/test/api/test_components/test_ext_grid.py b/src/pandapipes/test/api/test_components/test_ext_grid.py
index 7a3536f7..af3fe4aa 100644
--- a/src/pandapipes/test/api/test_components/test_ext_grid.py
+++ b/src/pandapipes/test/api/test_components/test_ext_grid.py
@@ -227,7 +227,7 @@ def test_t_type_tee_2zu_2ab(use_numba):
                                            u_w_per_m2k=5)
     pandapipes.create_fluid_from_lib(net, "water", overwrite=True)
     max_iter_hyd = 4 if use_numba else 4
-    max_iter_therm = 5 if use_numba else 5
+    max_iter_therm = 8 if use_numba else 8
     pandapipes.pipeflow(net, stop_condition="tol", max_iter_hyd=max_iter_hyd,
                         max_iter_therm=max_iter_therm, friction_model="nikuradse",
                         transient=False, nonlinear_method="automatic", tol_p=1e-4,
@@ -301,7 +301,7 @@ def test_t_type_tee_2zu_2ab2(use_numba):
                                            u_w_per_m2k=5)
     pandapipes.create_fluid_from_lib(net, "water", overwrite=True)
     max_iter_hyd = 3 if use_numba else 3
-    max_iter_therm = 6 if use_numba else 6
+    max_iter_therm = 8 if use_numba else 8
     pandapipes.pipeflow(net, stop_condition="tol", max_iter_hyd=max_iter_hyd,
                         max_iter_therm=max_iter_therm, friction_model="nikuradse",
                         transient=False, nonlinear_method="automatic", tol_p=1e-4,
@@ -377,7 +377,7 @@ def test_t_type_tee_2zu_2ab3(use_numba):
     pandapipes.create_fluid_from_lib(net, "water", overwrite=True)
 
     max_iter_hyd = 3 if use_numba else 3
-    max_iter_therm = 6 if use_numba else 6
+    max_iter_therm = 8 if use_numba else 8
     pandapipes.pipeflow(net, stop_condition="tol", max_iter_hyd=max_iter_hyd,
                         max_iter_therm=max_iter_therm, friction_model="nikuradse",
                         transient=False, nonlinear_method="automatic", tol_p=1e-4,
diff --git a/src/pandapipes/test/api/test_components/test_flow_control.py b/src/pandapipes/test/api/test_components/test_flow_control.py
index bfa3c0f2..b476db0d 100644
--- a/src/pandapipes/test/api/test_components/test_flow_control.py
+++ b/src/pandapipes/test/api/test_components/test_flow_control.py
@@ -18,15 +18,15 @@ def test_flow_control_simple_heat(use_numba):
     pandapipes.create_heat_exchanger(net, j3, j4, 0.1, 50000, 1)
     pandapipes.create_heat_exchanger(net, j6, j7, 0.1, 50000, 1)
 
-    pandapipes.create_flow_control(net, j2, j3, 2, 0.1)
-    pandapipes.create_flow_control(net, j5, j6, 2, 0.1, control_active=False)
+    pandapipes.create_flow_control(net, j2, j3, 2)
+    pandapipes.create_flow_control(net, j5, j6, 2, control_active=False)
 
     pandapipes.create_ext_grid(net, j1, p_bar=5, t_k=360, type="pt")
 
     pandapipes.create_sink(net, j8, 3)
 
     max_iter_hyd = 3 if use_numba else 3
-    max_iter_therm = 4 if use_numba else 4
+    max_iter_therm = 6 if use_numba else 6
     pandapipes.pipeflow(net, max_iter_hyd=max_iter_hyd, max_iter_therm=max_iter_therm,
                         mode='sequential', use_numba=use_numba)
 
diff --git a/src/pandapipes/test/api/test_components/test_heat_consumer.py b/src/pandapipes/test/api/test_components/test_heat_consumer.py
index c5caaf30..747281aa 100644
--- a/src/pandapipes/test/api/test_components/test_heat_consumer.py
+++ b/src/pandapipes/test/api/test_components/test_heat_consumer.py
@@ -30,13 +30,13 @@ def test_heat_consumer_equivalence(simple_heat_net):
     net2 = copy.deepcopy(simple_heat_net)
     juncs = net.junction.index
 
-    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=MDOT[0], qext_w=QEXT[0])
-    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=MDOT[1], qext_w=QEXT[1])
+    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], controlled_mdot_kg_per_s=MDOT[0], qext_w=QEXT[0])
+    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], controlled_mdot_kg_per_s=MDOT[1], qext_w=QEXT[1])
     pandapipes.pipeflow(net, mode='sequential')
 
     j_mid = pandapipes.create_junctions(net2, 2, pn_bar=5, tfluid_k=283.15)
     pandapipes.create_flow_controls(net2, juncs[[1, 2]], j_mid, MDOT, diameter_m=0.1022)
-    pandapipes.create_heat_exchangers(net2, j_mid, juncs[[4, 3]], 0.1022, QEXT)
+    pandapipes.create_heat_exchangers(net2, j_mid, juncs[[4, 3]], qext_w=QEXT)
     pandapipes.pipeflow(net2, mode='sequential')
 
     assert np.allclose(net.res_junction.values, net2.res_junction.iloc[:-2, :].values)
@@ -47,13 +47,12 @@ def test_heat_consumer_equivalence_bulk(simple_heat_net):
     net2 = copy.deepcopy(simple_heat_net)
     juncs = net.junction.index
 
-    pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, controlled_mdot_kg_per_s=MDOT,
-                                     qext_w=QEXT)
+    pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], controlled_mdot_kg_per_s=MDOT, qext_w=QEXT)
     pandapipes.pipeflow(net, mode='sequential')
 
     j_mid = pandapipes.create_junctions(net2, 2, pn_bar=5, tfluid_k=283.15)
     pandapipes.create_flow_controls(net2, juncs[[1, 2]], j_mid, MDOT, diameter_m=0.1022)
-    pandapipes.create_heat_exchangers(net2, j_mid, juncs[[4, 3]], 0.1022, QEXT)
+    pandapipes.create_heat_exchangers(net2, j_mid, juncs[[4, 3]], qext_w=QEXT)
     pandapipes.pipeflow(net2, mode='sequential')
 
     assert np.allclose(net.res_junction.values, net2.res_junction.iloc[:-2, :].values)
@@ -76,27 +75,27 @@ def test_heat_consumer_equivalence2():
     net4 = copy.deepcopy(net)
     net5 = copy.deepcopy(net)
 
-    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
-    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=mdot[1], qext_w=qext[1])
-    pandapipes.pipeflow(net, mode="bidirectional")
+    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
+    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], controlled_mdot_kg_per_s=mdot[1], qext_w=qext[1])
+    pandapipes.pipeflow(net, mode="bidirectional", iter=6)
 
-    pandapipes.create_heat_consumer(net2, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
-    pandapipes.create_heat_consumer(net2, juncs[2], juncs[3], 0.1022, treturn_k=263.4459264973806, qext_w=qext[1])
-    pandapipes.pipeflow(net2, mode="bidirectional", iter=25)
+    pandapipes.create_heat_consumer(net2, juncs[1], juncs[4], controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
+    pandapipes.create_heat_consumer(net2, juncs[2], juncs[3], treturn_k=263.4459264973806, qext_w=qext[1])
+    pandapipes.pipeflow(net2, mode="bidirectional", iter=23)
 
-    pandapipes.create_heat_consumer(net3, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
-    pandapipes.create_heat_consumer(net3, juncs[2], juncs[3], 0.1022, deltat_k=17.82611044059695, qext_w=qext[1])
-    pandapipes.pipeflow(net3, mode="bidirectional")
+    pandapipes.create_heat_consumer(net3, juncs[1], juncs[4], controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
+    pandapipes.create_heat_consumer(net3, juncs[2], juncs[3], deltat_k=17.82611044059695, qext_w=qext[1])
+    pandapipes.pipeflow(net3, mode="bidirectional", iter=6)
 
-    pandapipes.create_heat_consumer(net4, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
-    pandapipes.create_heat_consumer(net4, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=mdot[1],
+    pandapipes.create_heat_consumer(net4, juncs[1], juncs[4], controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
+    pandapipes.create_heat_consumer(net4, juncs[2], juncs[3], controlled_mdot_kg_per_s=mdot[1],
                                     treturn_k=263.4459264973806)
-    pandapipes.pipeflow(net4, mode="bidirectional")
+    pandapipes.pipeflow(net4, mode="bidirectional", iter=7)
 
-    pandapipes.create_heat_consumer(net5, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
-    pandapipes.create_heat_consumer(net5, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=mdot[1],
+    pandapipes.create_heat_consumer(net5, juncs[1], juncs[4], controlled_mdot_kg_per_s=mdot[0], qext_w=qext[0])
+    pandapipes.create_heat_consumer(net5, juncs[2], juncs[3], controlled_mdot_kg_per_s=mdot[1],
                                     deltat_k=17.82611044059695)
-    pandapipes.pipeflow(net5, mode="bidirectional")
+    pandapipes.pipeflow(net5, mode="bidirectional", iter=6)
 
     assert np.allclose(net2.res_junction, net.res_junction)
     assert np.allclose(net2.res_pipe, net.res_pipe)
@@ -114,19 +113,19 @@ def test_heat_consumer_creation_not_allowed(simple_heat_net):
 
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters
-        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=MDOT[0], qext_w=None,
+        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], controlled_mdot_kg_per_s=MDOT[0], qext_w=None,
                                         treturn_k=None)
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, controlled_mdot_kg_per_s=MDOT,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], controlled_mdot_kg_per_s=MDOT,
                                          qext_w=[QEXT[0], None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, qext_w=QEXT,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], qext_w=QEXT,
                                          controlled_mdot_kg_per_s=[MDOT[0], None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, qext_w=QEXT,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], qext_w=QEXT,
                                          controlled_mdot_kg_per_s=None)
 
 
@@ -135,41 +134,41 @@ def test_heat_consumer_creation_not_allowed_2(simple_heat_net):
     juncs = net.junction.index
     with pytest.raises(AttributeError):
         # check for more than 2 set parameters
-        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, controlled_mdot_kg_per_s=MDOT[0],
+        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], controlled_mdot_kg_per_s=MDOT[0],
                                         qext_w=QEXT[0], treturn_k=390)
     with pytest.raises(AttributeError):
         # check for deltat_k and treturn_k given
-        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, deltat_k=20, treturn_k=390)
+        pandapipes.create_heat_consumer(net, juncs[1], juncs[4], deltat_k=20, treturn_k=390)
 
     with pytest.raises(AttributeError):
         # check for more than 2 set parameters
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, deltat_k=[30, 40],
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], deltat_k=[30, 40],
                                          treturn_k=[390, 385])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, controlled_mdot_kg_per_s=MDOT,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], controlled_mdot_kg_per_s=MDOT,
                                          deltat_k=[20, None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, qext_w=QEXT, deltat_k=[20, None])
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], qext_w=QEXT, deltat_k=[20, None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, controlled_mdot_kg_per_s=MDOT,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], controlled_mdot_kg_per_s=MDOT,
                                          treturn_k=[390, None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in some consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, qext_w=QEXT, treturn_k=[390, None])
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], qext_w=QEXT, treturn_k=[390, None])
     with pytest.raises(AttributeError):
         # check for less than 2 set parameters in all consumers
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, qext_w=QEXT, treturn_k=None,
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], qext_w=QEXT, treturn_k=None,
                                          controlled_mdot_kg_per_s=None)
     with pytest.raises(AttributeError):
         # check for deltat_k and treturn_k given
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, deltat_k=[30, 40],
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], deltat_k=[30, 40],
                                          treturn_k=[390, 385])
     with pytest.raises(AttributeError):
         # check for deltat_k and treturn_k given as single values
-        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], 0.1022, deltat_k=30, treturn_k=390)
+        pandapipes.create_heat_consumers(net, juncs[[1, 2]], juncs[[4, 3]], deltat_k=30, treturn_k=390)
 
 
 def test_heat_consumer_qext_zero():
@@ -181,8 +180,8 @@ def test_heat_consumer_qext_zero():
                                             text_k=273.15)
     pandapipes.create_circ_pump_const_pressure(net, juncs[-1], juncs[0], 5, 2, 300, type='pt')
 
-    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, treturn_k=263.4459264973806, qext_w=0)
-    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=1, qext_w=7500)
+    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], treturn_k=263.4459264973806, qext_w=0)
+    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], controlled_mdot_kg_per_s=1, qext_w=7500)
 
     pandapipes.pipeflow(net, mode="bidirectional")
 
@@ -196,18 +195,18 @@ def test_heat_consumer_result_extraction():
                                             diameter_m=0.1022, system=["flow"] * 2 + ["return"] * 2, alpha_w_per_m2k=10,
                                             text_k=273.15)
     pandapipes.create_circ_pump_const_pressure(net, juncs[-1], juncs[0], 5, 2, 300, type='pt')
-    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], 0.1022, treturn_k=263.4459264973806, qext_w=7500)
-    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], 0.1022, controlled_mdot_kg_per_s=1, qext_w=7500)
+    pandapipes.create_heat_consumer(net, juncs[1], juncs[4], treturn_k=263.4459264973806, qext_w=7500)
+    pandapipes.create_heat_consumer(net, juncs[2], juncs[3], controlled_mdot_kg_per_s=1, qext_w=7500)
 
     # create not connected pipe to test for active inactive missmatch
     pandapipes.create_junctions(net, 2, pn_bar=5, tfluid_k=286)
     pandapipes.create_pipe_from_parameters(net, 6, 7, k_mm=0.1, length_km=1,
                                            diameter_m=0.1022, alpha_w_per_m2k=10, text_k=273.15)
 
-    pandapipes.pipeflow(net, mode="bidirectional")
+    pandapipes.pipeflow(net, mode="bidirectional", iter=13)
 
     #hydraulics only to check for lookup heat transfer error
-    pandapipes.pipeflow(net)
+    pandapipes.pipeflow(net, iter=3)
 
 
 if __name__ == '__main__':
diff --git a/src/pandapipes/test/api/test_components/test_heat_exchanger.py b/src/pandapipes/test/api/test_components/test_heat_exchanger.py
index fbb27f06..06686c79 100644
--- a/src/pandapipes/test/api/test_components/test_heat_exchanger.py
+++ b/src/pandapipes/test/api/test_components/test_heat_exchanger.py
@@ -19,11 +19,10 @@ def test_heat_exchanger(use_numba):
         :rtype:
         """
     net = pandapipes.create_empty_network("net", add_stdtypes=False)
-    d = 75e-3
 
     pandapipes.create_junction(net, pn_bar=5, tfluid_k=283)
     pandapipes.create_junction(net, pn_bar=5, tfluid_k=283)
-    pandapipes.create_heat_exchanger(net, 0, 1, d, qext_w=20000)
+    pandapipes.create_heat_exchanger(net, 0, 1, qext_w=20000)
     pandapipes.create_ext_grid(net, 0, p_bar=5, t_k=330, type="pt")
     pandapipes.create_sink(net, 1, mdot_kg_per_s=1)
 
diff --git a/src/pandapipes/test/api/test_convert_format.py b/src/pandapipes/test/api/test_convert_format.py
index 59380702..2edd9be1 100644
--- a/src/pandapipes/test/api/test_convert_format.py
+++ b/src/pandapipes/test/api/test_convert_format.py
@@ -21,32 +21,35 @@
 from pandapipes.test.api.release_cycle.release_control_test_network import release_control_test_network_gas, \
     release_control_test_network_water, release_control_test_network
 
-
 @pytest.mark.slow
 @pytest.mark.parametrize("pp_version, use_numba", test_params)
 def test_convert_format(pp_version, use_numba):
     if version.parse(pp_version) >= version.parse(minimal_version_two_nets):
         names = ["_gas", "_water"]
         net_gas = release_control_test_network_gas(max_iter_hyd=5)
-        net_water = release_control_test_network_water(max_iter_hyd=10)
+        net_water = release_control_test_network_water(max_iter_hyd=11)
     else:
         names = [""]
-        net_old = release_control_test_network(max_iter_hyd=10)
+        net_old = release_control_test_network(max_iter_hyd=11)
     for name in names:
         if "_gas" in name:
             net_ref = net_gas
             max_iter_hyd = 5 if use_numba else 5
         elif "_water" in name:
             net_ref = net_water
-            max_iter_hyd = 10 if use_numba else 10
+            max_iter_hyd = 11 if use_numba else 11
         else:
             net_ref = net_old
-            max_iter_hyd = 10 if use_numba else 10
+            max_iter_hyd = 11 if use_numba else 11
         filename = os.path.join(folder, "example_%s%s.json" % (pp_version, name))
         if not os.path.isfile(filename):
             raise ValueError("File for %s grid of version %s does not exist" % (name, pp_version))
         try:
             net = pp.from_json(filename, convert=False)
+            if not "_gas" in name:
+                pp.create_ext_grid(net, junction=4, p_bar=6, t_k=290, name="External Grid 2", index=None)
+                pp.create_ext_grid(net, junction=5, p_bar=5, t_k=290, name="External Grid 3")
+
         except:
             raise UserWarning("Can not load %s network saved in pandapipes version %s"
                               % (name, pp_version))
diff --git a/src/pandapipes/test/api/test_create.py b/src/pandapipes/test/api/test_create.py
index 3db40790..c2dc946a 100644
--- a/src/pandapipes/test/api/test_create.py
+++ b/src/pandapipes/test/api/test_create.py
@@ -94,21 +94,20 @@ def test_create_heat_exchanger(create_empty_net):
     net = copy.deepcopy(create_empty_net)
     pandapipes.create_junction(net, 1, 293, index=8, geodata=(0, 1))
     pandapipes.create_junction(net, 1, 293, index=9, geodata=(2, 2))
-    pandapipes.create_heat_exchanger(net, 8, 9, 0.3, qext_w=200, index=2)
+    pandapipes.create_heat_exchanger(net, 8, 9, qext_w=200, index=2)
 
     assert len(net.junction) == 2
     assert len(net.heat_exchanger) == 1
     assert np.all(net.heat_exchanger.index == [2])
     assert net.heat_exchanger.at[2, "from_junction"] == 8
     assert net.heat_exchanger.at[2, "to_junction"] == 9
-    assert net.heat_exchanger.at[2, "diameter_m"] == 0.3
     assert net.heat_exchanger.at[2, "qext_w"] == 200
     assert net.heat_exchanger.at[2, "loss_coefficient"] == 0
 
     with pytest.raises(UserWarning):
-        pandapipes.create_heat_exchanger(net, 8, 10, 0.3, qext_w=200)
+        pandapipes.create_heat_exchanger(net, 8, 10, qext_w=200)
     with pytest.raises(UserWarning):
-        pandapipes.create_heat_exchanger(net, 8, 9, 0.3, qext_w=200, index=2)
+        pandapipes.create_heat_exchanger(net, 8, 9, qext_w=200, index=2)
 
 
 def test_create_pipe(create_empty_net):
diff --git a/src/pandapipes/test/api/test_network_tables.py b/src/pandapipes/test/api/test_network_tables.py
index 04ff9c8e..ce5bf389 100644
--- a/src/pandapipes/test/api/test_network_tables.py
+++ b/src/pandapipes/test/api/test_network_tables.py
@@ -77,7 +77,7 @@ def test_additional_tables():
 
     pandapipes.add_new_component(net, HeatExchanger)
     hex_input = list(copy.deepcopy(net.heat_exchanger.columns))
-    pandapipes.create_heat_exchanger(net, 0, 1, 0.2, qext_w=20000)
+    pandapipes.create_heat_exchanger(net, 0, 1, qext_w=20000)
     hex_input_create = list(net.heat_exchanger.columns)
     assert hex_input == hex_input_create, "Input does not equal create-table"
 
diff --git a/src/pandapipes/test/api/test_special_networks.py b/src/pandapipes/test/api/test_special_networks.py
index 90e63497..b04769f8 100644
--- a/src/pandapipes/test/api/test_special_networks.py
+++ b/src/pandapipes/test/api/test_special_networks.py
@@ -28,7 +28,7 @@ def test_one_node_net(use_numba):
     create_sink(net, j, 0.01)
     create_source(net, j, 0.02)
 
-    max_iter_hyd = 1 if use_numba else 1
+    max_iter_hyd = 2 if use_numba else 2
     pandapipes.pipeflow(net, use_numba=use_numba, max_iter_hyd=max_iter_hyd)
 
     assert np.isclose(net.res_ext_grid.values + net.res_sink.values - net.res_source.values, 0)
@@ -60,7 +60,7 @@ def test_two_node_net(use_numba):
     create_ext_grid(net, j, 1, 298.15)
     create_sink(net, j, 0.01)
     create_source(net, j, 0.02)
-    max_iter_hyd = 1 if use_numba else 1
+    max_iter_hyd = 2 if use_numba else 2
     pandapipes.pipeflow(net, use_numba=use_numba, max_iter_hyd=max_iter_hyd)
 
     assert np.all(np.isclose(net.res_ext_grid.values + net.res_sink.values - net.res_source.values,
diff --git a/src/pandapipes/test/openmodelica_comparison/pipeflow_openmodelica_comparison.py b/src/pandapipes/test/openmodelica_comparison/pipeflow_openmodelica_comparison.py
index c70eaa1c..9169f786 100644
--- a/src/pandapipes/test/openmodelica_comparison/pipeflow_openmodelica_comparison.py
+++ b/src/pandapipes/test/openmodelica_comparison/pipeflow_openmodelica_comparison.py
@@ -171,14 +171,8 @@ def retrieve_velocity_liquid(net, element="pipe"):
     v_valid = pd.notnull(v_om)
     v_om = v_om.loc[v_valid]
     v_om[v_om == 0] += 0.0001
-
-    if element == "pipe":
-        v_mean_pandapipes = net.res_pipe.v_mean_m_per_s.loc[v_valid].values.astype(
-            np.float64).round(4)
-
-    if element == "valve":
-        v_mean_pandapipes = net.res_valve.v_mean_m_per_s.loc[v_valid].values.astype(
-            np.float64).round(4)
+    v_mean_pandapipes = net['res_%s' %element].v_mean_m_per_s.loc[v_valid].values.astype(
+        np.float64).round(4)
 
     v_mean_pandapipes[v_mean_pandapipes == 0] += 0.0001
 
diff --git a/src/pandapipes/test/openmodelica_comparison/test_water_openmodelica.py b/src/pandapipes/test/openmodelica_comparison/test_water_openmodelica.py
index da8de0c4..982fac1b 100644
--- a/src/pandapipes/test/openmodelica_comparison/test_water_openmodelica.py
+++ b/src/pandapipes/test/openmodelica_comparison/test_water_openmodelica.py
@@ -46,7 +46,7 @@ def test_case_combined_mixed_sj(use_numba, log_results=False):
 @pytest.mark.parametrize("use_numba", [True, False])
 def test_case_combined_versatility_pc(use_numba, log_results=False):
     net = nw.water_combined_versatility()
-    max_iter_hyd = 14 if use_numba else 14
+    max_iter_hyd = 15 if use_numba else 15
     p_diff, v_diff_abs = pipeflow_openmodelica_comparison(
         net, log_results,
         max_iter_hyd=max_iter_hyd,
@@ -58,7 +58,7 @@ def test_case_combined_versatility_pc(use_numba, log_results=False):
 @pytest.mark.parametrize("use_numba", [True, False])
 def test_case_combined_versatility_sj(use_numba, log_results=False):
     net = nw.water_combined_versatility(method="swamee-jain")
-    max_iter_hyd = 14 if use_numba else 14
+    max_iter_hyd = 15 if use_numba else 15
     p_diff, v_diff_abs = pipeflow_openmodelica_comparison(
         net, log_results,
         max_iter_hyd=max_iter_hyd,
diff --git a/src/pandapipes/test/pipeflow_internals/test_inservice.py b/src/pandapipes/test/pipeflow_internals/test_inservice.py
index 2079318a..6827e923 100644
--- a/src/pandapipes/test/pipeflow_internals/test_inservice.py
+++ b/src/pandapipes/test/pipeflow_internals/test_inservice.py
@@ -412,7 +412,7 @@ def test_connectivity_heat4(complex_heat_connectivity_grid, use_numba):
 
     net2 = copy.deepcopy(net)
 
-    max_iter_hyd = 9 if use_numba else 9
+    max_iter_hyd = 10 if use_numba else 10
     max_iter_therm = 6 if use_numba else 6
     pandapipes.pipeflow(net, max_iter_hyd=max_iter_hyd, max_iter_therm=max_iter_therm,
                         mode='sequential', check_connectivity=True, use_numba=use_numba)
@@ -436,7 +436,7 @@ def test_connectivity_heat5(complex_heat_connectivity_grid, use_numba):
     net.ext_grid.loc[2, 'in_service'] = False
     net.ext_grid.loc[1, 'type'] = 'p'
 
-    max_iter_hyd = 9 if use_numba else 9
+    max_iter_hyd = 10 if use_numba else 10
     max_iter_therm = 3 if use_numba else 3
     pandapipes.pipeflow(net, max_iter_hyd=max_iter_hyd, max_iter_therm=max_iter_therm,
                         check_connectivity=True, mode='sequential', use_numba=use_numba)
diff --git a/src/pandapipes/test/pipeflow_internals/test_pipeflow_analytic_comparison.py b/src/pandapipes/test/pipeflow_internals/test_pipeflow_analytic_comparison.py
index 7e7e4149..1212961a 100644
--- a/src/pandapipes/test/pipeflow_internals/test_pipeflow_analytic_comparison.py
+++ b/src/pandapipes/test/pipeflow_internals/test_pipeflow_analytic_comparison.py
@@ -390,7 +390,7 @@ def test_temperature_internal_nodes_2zu_2ab(use_numba):
     pandapipes.create_fluid_from_lib(net, "water", overwrite=True)
 
     max_iter_hyd = 3 if use_numba else 3
-    max_iter_therm = 4 if use_numba else 4
+    max_iter_therm = 5 if use_numba else 5
     pandapipes.pipeflow(
         net,
         stop_condition="tol",
@@ -452,7 +452,7 @@ def test_temperature_internal_nodes_masche_1load(use_numba):
     pandapipes.create_fluid_from_lib(net, "water", overwrite=True)
 
     max_iter_hyd = 3 if use_numba else 3
-    max_iter_therm = 4 if use_numba else 4
+    max_iter_therm = 5 if use_numba else 5
     pandapipes.pipeflow(
         net,
         stop_condition="tol",
diff --git a/src/pandapipes/test/plotting/test_collections.py b/src/pandapipes/test/plotting/test_collections.py
index 5c3cd901..c6091bb8 100644
--- a/src/pandapipes/test/plotting/test_collections.py
+++ b/src/pandapipes/test/plotting/test_collections.py
@@ -47,11 +47,11 @@ def test_collection_lengths():
                                            geodata=[(4, 0), (6, 0)])
     pandapipes.create_pipe_from_parameters(net, j5, j6, 0.2, diameter_m=d, k_mm=0.1,
                                            geodata=[(4, -2), (6, -2), (6, 0)])
-    pandapipes.create_heat_exchanger(net, j6, j7, d, qext_w=20000)
+    pandapipes.create_heat_exchanger(net, j6, j7, qext_w=20000)
     pandapipes.create_pump_from_parameters(net, j7, j8, 'P1')
     pandapipes.create_pressure_control(net, j8, j9, j9, 10.)
-    pandapipes.create_flow_control(net, j9, j10, 0.5, 0.1)
-    pandapipes.create_flow_control(net, j9, j11, 0.5, 0.1, control_active=False)
+    pandapipes.create_flow_control(net, j9, j10, 0.5)
+    pandapipes.create_flow_control(net, j9, j11, 0.5, control_active=False)
 
     pipe_coll_direct = plot.create_pipe_collection(net, use_junction_geodata=True)
     pipe_coll_real = plot.create_pipe_collection(net)
diff --git a/src/pandapipes/test/test_toolbox.py b/src/pandapipes/test/test_toolbox.py
index 27111a3a..91cefad3 100644
--- a/src/pandapipes/test/test_toolbox.py
+++ b/src/pandapipes/test/test_toolbox.py
@@ -274,7 +274,10 @@ def test_pit_extraction():
     for name in names:
         filename = os.path.join(folder, "example_%s%s.json" % (max_ver, name))
         net = pandapipes.from_json(filename)
-        max_iter_hyd = 10 if '_water' in name else 5
+        if not "_gas" in name:
+            pandapipes.create_ext_grid(net, junction=4, p_bar=6, t_k=290, name="External Grid 2", index=None)
+            pandapipes.create_ext_grid(net, junction=5, p_bar=5, t_k=290, name="External Grid 3")
+        max_iter_hyd = 11 if '_water' in name else 5
         pandapipes.pipeflow(net, max_iter_hyd=max_iter_hyd)
 
         node_table, branch_table = pandapipes.get_internal_tables_pandas(net)
diff --git a/src/pandapipes/toolbox.py b/src/pandapipes/toolbox.py
index b3850894..610dfb0c 100644
--- a/src/pandapipes/toolbox.py
+++ b/src/pandapipes/toolbox.py
@@ -18,7 +18,7 @@
 from pandapipes.create import create_empty_network
 from pandapipes.idx_branch import branch_cols
 from pandapipes.idx_node import node_cols, \
-    T as TYPE_T, P as TYPE_P, PC as TYPE_PC, NONE as TYPE_NONE, L as TYPE_L
+    T as TYPE_T, P as TYPE_P, PC as TYPE_PC, L as TYPE_L
 from pandapipes.pandapipes_net import pandapipesNet
 from pandapipes.topology import create_nxgraph
 
@@ -526,7 +526,7 @@ def check_pressure_controllability(net, to_junction, controlled_junction):
 #     logger.info("dropped %d %s elements with %d switches" % (len(trafos), table, num_switches))
 
 
-pit_types = {TYPE_P: "P", TYPE_L: "L", TYPE_NONE: "NONE", TYPE_T: "T", TYPE_PC: "PC", 0: "NONE"}
+pit_types = {TYPE_P: "P", TYPE_L: "L", TYPE_T: "T", TYPE_PC: "PC"}
 int_cols = ["FROM_NODE", "TO_NODE", "ELEMENT_IDX", "EXT_GRID_OCCURENCE", "EXT_GRID_OCCURENCE_T"]
 bool_cols = ["ACTIVE"]